Optimal low-rank approximations of Bayesian linear inverse problems

In the Bayesian approach to inverse problems, data are often informative, relative to the prior, only on a low-dimensional subspace of the parameter space. Significant computational savings can be achieved by using this subspace to characterize and approximate the posterior distribution of the parameters. We first investigate approximation of the posterior covariance matrix as a low-rank update of the prior covariance matrix. We prove optimality of a particular update, based on the leading eigendirections of the matrix pencil defined by the Hessian of the negative log-likelihood and the prior precision, for a broad class of loss functions. This class includes the F\"{o}rstner metric for symmetric positive definite matrices, as well as the Kullback-Leibler divergence and the Hellinger distance between the associated distributions. We also propose two fast approximations of the posterior mean and prove their optimality with respect to a weighted Bayes risk under squared-error loss. These approximations are deployed in an offline-online manner, where a more costly but data-independent offline calculation is followed by fast online evaluations. As a result, these approximations are particularly useful when repeated posterior mean evaluations are required for multiple data sets. We demonstrate our theoretical results with several numerical examples, including high-dimensional X-ray tomography and an inverse heat conduction problem. In both of these examples, the intrinsic low-dimensional structure of the inference problem can be exploited while producing results that are essentially indistinguishable from solutions computed in the full space.

Similar Publications

We describe a simple and effective technique, the Eigenvector Method for Umbrella Sampling (EMUS), for accurately estimating small probabilities and expectations with respect to a given target probability density. In EMUS, we apply the principle of stratified survey sampling to Markov chain Monte Carlo (MCMC) simulation: We divide the support of the target distribution into regions called strata, we use MCMC to sample (in parallel) from probability distributions supported in each of the strata, and we weight the data from each stratum to assemble estimates of general averages with respect to the target distribution. We demonstrate by theoretical results and computational examples that EMUS can be dramatically more efficient than direct Markov chain Monte Carlo when the target distribution is multimodal or when the goal is to compute tail probabilities. Read More

We develop a finite element method for elliptic partial differential equations on so called composite surfaces that are built up out of a finite number of surfaces with boundaries that fit together nicely in the sense that the intersection between any two surfaces in the composite surface is either empty, a point, or a curve segment, called an interface curve. Note that several surfaces can intersect along the same interface curve. On the composite surface we consider a broken finite element space which consists of a continuous finite element space at each subsurface without continuity requirements across the interface curves. Read More

We investigate a finite element formulation of the exponentiated Hencky-logarithmic model whose strain energy function is given by \[ W_\mathrm{eH}(\boldsymbol{F}) = \dfrac{\mu}{k}\, e^{\displaystyle k \left\lVert\mbox{dev}_n \log\boldsymbol{U}\right\rVert^2} + \dfrac{\kappa}{2 \hat{k}}\, e^{\displaystyle \hat{k} [\mbox{tr} (\log\boldsymbol{U})]^2 }\,, \] where $\mu>0$ is the (infinitesimal) shear modulus, $\kappa>0$ is the (infinitesimal) bulk modulus, $k$ and $\hat{k}$ are additional dimensionless material parameters, $\boldsymbol{U}=\sqrt{\boldsymbol{F}^T\boldsymbol{F}}$ and $\boldsymbol{V}=\sqrt{\boldsymbol{F}\boldsymbol{F}^T}$ are the right and left stretch tensor corresponding to the deformation gradient $\boldsymbol{F}$, $\log$ denotes the principal matrix logarithm on the set of positive definite symmetric matrices, $\mbox{dev}_n \boldsymbol{X} = \boldsymbol{X}-\frac{\mbox{tr} \boldsymbol{X}}{n}\boldsymbol{1}$ and $\lVert \boldsymbol{X} \rVert = \sqrt{\mbox{tr}\boldsymbol{X}^T\boldsymbol{X}}$ are the deviatoric part and the Frobenius matrix norm of an $n\times n$-matrix $\boldsymbol{X}$, respectively, and $\mbox{tr}$ denotes the trace operator. To do so, the equivalent different forms of the constitutive equation are recast in terms of the principal logarithmic stretches by use of the spectral decomposition together with the undergoing properties. We show the capability of our approach with a number of relevant examples, including the challenging "eversion of elastic tubes" problem. Read More

The Poisson-Boltzmann equation (PBE) is a nonlinear elliptic parametrized partial differential equation that arises in biomolecular modeling and is a fundamental tool for structural biology. It is used to calculate electrostatic potentials around an ensemble of fixed charges immersed in an ionic solution. Efficient numerical computation of the PBE yields a high number of degrees of freedom in the resultant algebraic system of equations, ranging from several hundred thousands to millions. Read More

In this paper, we propose an efficient numerical scheme for the approximate solution of the time fractional diffusion-wave equation with reaction term based on cubic trigonometric basis functions. The time fractional derivative is approximated by the usual finite difference formulation and the derivative in space is discretized using cubic trigonometric B-spline functions. A stability analysis of the scheme is conducted to confirm that the scheme does not amplify errors. Read More

In this paper, we present the compact representation for matrices belonging to the the Broyden class of quasi-Newton updates, where each update may be either rank-one or rank-two. This work extends previous results solely for the restricted Broyden class of rank-two updates. In this article, it is not assumed the same Broyden update is used each iteration; rather, different members of the Broyden class may be used each iteration. Read More

Mathematical methods of image inpainting involve the discretization of given continuous models. We present a method that avoids the standard pointwise discretization by modeling known variational approaches, in particular total variation (TV), using a finite dimensional spline space. Besides the analysis of the resulting model, we present a numerical implementation based on the alternating method of multipliers. Read More

Sparse data approximation has become a popular research topic in signal processing. However, in most cases only a single measurement vector (SMV) is considered. In applications, the multiple measurement vector (MMV) case is more usual, i. Read More

In this paper we compare numerically two different coarse space definitions for two-level domain decomposition preconditioners for the Helmholtz equation, both in two and three dimensions. While we solve the pure Helmholtz problem without absorption, the preconditioners are built from problems with absorption. In the first method, the coarse space is based on the discretization of the problem with absorption on a coarse mesh, with diameter constrained by the wavenumber. Read More

The construction of fast iterative solvers for the indefinite time-harmonic Maxwell's system at high-frequency is a problem of great current interest. Some of the difficulties that arise are similar to those encountered in the case of the high-frequency Helmholtz equation. Here we investigate how two-level domain-decomposition preconditioners recently proposed for the Helmholtz equation work in the Maxwell case, both from the theoretical and numerical points of view. Read More