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

When modeling global satellite data to recover a planetary magnetic or gravitational potential field and evaluate it elsewhere, the method of choice remains their analysis in terms of spherical harmonics. When only regional data are available, or when data quality varies strongly with geographic location, the inversion problem becomes severely ill-posed. In those cases, adopting explicitly local methods is to be preferred over adapting global ones (e. Read More


It is well known that the equation $x'(t)=Ax(t)+f(t)$, where $A$ is a square matrix, has a unique bounded solution $x$ for any bounded continuous free term $f$, provided the coefficient $A$ has no eigenvalues on the imaginary axis. This solution can be represented in the form \begin{equation*} x(t)=\int_{-\infty}^{\infty}\mathcal G(t-s)x(s)\,ds. \end{equation*} The kernel $\mathcal G$ is called Green's function. Read More


This article reviews the application of advanced Monte Carlo techniques in the context of Multilevel Monte Carlo (MLMC). MLMC is a strategy employed to compute expectations which can be biased in some sense, for instance, by using the discretization of a associated probability law. The MLMC approach works with a hierarchy of biased approximations which become progressively more accurate and more expensive. Read More


We outline the construction of compatible B-splines on 3D surfaces that satisfy the continuity requirements for electromagnetic scattering analysis with the boundary element method (method of moments). Our approach makes use of Non-Uniform Rational B-splines to represent model geometry and compatible B-splines to approximate the surface current, and adopts the isogeometric concept in which the basis for analysis is taken directly from CAD (geometry) data. The approach allows for high-order approximations and crucially provides a direct link with CAD data structures that allows for efficient design workflows. Read More


The goal of this paper is to achieve a computational model and corresponding efficient algorithm for obtaining a sparse representation of the fitting surface to the given scattered data. The basic idea of the model is to utilize the principal shift invariant (PSI) space and the l1 norm minimization. In order to obtain different sparsity of the approximation solution, the problem is represented as a multilevel LASSO (MLASSO) model with different regularization parameters. Read More


This analysis proposes an analytical-numerical approach for providing solutions of a class of nonlinear fractional Klein-Gordon equation subjected to appropriate initial conditions in Caputo sense by using the Fractional Reduced Differential Transform Method (FRDTM). This technique provides the solutions very accurately and efficiently in convergent series formula with easily computable coefficients. The behavior of the approximate series solution for different values of fractional-order "a" is shown graphically. Read More


We introduce a novel nonlinear imaging method for the acoustic wave equation based on model order reduction. The objective is to image the discontinuities of the acoustic velocity, a coefficient of the scalar wave equation from the discretely sampled time domain data measured at an array of transducers that can act as both sources and receivers. We treat the wave equation along with transducer functionals as a dynamical system. Read More


In engineering, it is a common desire to couple existing simulation tools together into one big system by passing information from subsystems as parameters into the subsystems under influence. As executed at fixed time points, this data exchange gives the global method a strong explicit component. Globally, such an explicit cosimulation schemes exchange time step can be seen as a step of an one-step method which is explicit in some solution components. Read More


The interval subset sum problem (ISSP) is a generalization of the well-known subset sum problem. Given a set of intervals $\left\{[a_{i,1},a_{i,2}]\right\}_{i=1}^n$ and a target integer $T,$ the ISSP is to find a set of integers, at most one from each interval, such that their sum best approximates the target $T$ but cannot exceed it. In this paper, we first study the computational complexity of the ISSP. Read More


The numerical approximation of high-frequency wave propagation in inhomogeneous media is a challenging problem. In particular, computing high-frequency solutions by direct simulations requires several points per wavelength for stability and usually requires many points per wavelength for a satisfactory accuracy. In this paper, we propose a new method for the acoustic wave equation in inhomogeneous media in the time domain to achieve superior accuracy and stability without using a large number of unknowns. Read More