Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Subspace Splitting Fast Sampling from Gaussian Posterior Distributions of Linear Inverse Problems

Daniela Calvetti, Erkki Somersalo

TL;DR

This work tackles efficient sampling from Gaussian posteriors in high-dimensional linear inverse problems, where the data dimension is often far smaller than the unknown parameter dimension. It introduces a splitting-based Randomize-Then-Optimize (RTO) approach that leverages a low-dimensional data-space adjoint problem to generate independent posterior samples without forming the full covariance, and it supports matrix-free forward models via Krylov/Lanczos techniques. The method generalizes to non-identity priors through whitening and linear transformations, enabling sampling in transformed coordinates and backward-translation to the original variables, with extensions to conditionally Gaussian hierarchies. Demonstrations on tomography (fanbeam, cross-borehole), electrical impedance tomography, and magnetoencephalography show substantial computational gains and practical applicability for uncertainty quantification and proposal generation in Bayesian workflows.

Abstract

It is well-known that the posterior density of linear inverse problems with Gaussian prior and Gaussian likelihood is also Gaussian, hence completely described by its covariance and expectation. Sampling from a Gaussian posterior may be important in the analysis of various non-Gaussian inverse problems in which a estimates from a Gaussian posterior distribution constitute an intermediate stage in a Bayesian workflow. Sampling from a Gaussian distribution is straightforward if the Cholesky factorization of the covariance matrix or its inverse is available, however when the unknown is high dimensional, the computation of the posterior covariance maybe unfeasible. If the linear inverse problem is underdetermined, it is possible to exploit the orthogonality of the fundamental subspaces associated with the coefficient matrix together with the idea behind the Randomize-Then-Optimize approach to design a low complexity posterior sampler that does not require the posterior covariance to be formed. The performance of the proposed sampler is illustrated with a few computed examples, including non-Gaussian problems with non-linear forward model, and hierarchical models comprising a conditionally Gaussian submodel.

Subspace Splitting Fast Sampling from Gaussian Posterior Distributions of Linear Inverse Problems

TL;DR

This work tackles efficient sampling from Gaussian posteriors in high-dimensional linear inverse problems, where the data dimension is often far smaller than the unknown parameter dimension. It introduces a splitting-based Randomize-Then-Optimize (RTO) approach that leverages a low-dimensional data-space adjoint problem to generate independent posterior samples without forming the full covariance, and it supports matrix-free forward models via Krylov/Lanczos techniques. The method generalizes to non-identity priors through whitening and linear transformations, enabling sampling in transformed coordinates and backward-translation to the original variables, with extensions to conditionally Gaussian hierarchies. Demonstrations on tomography (fanbeam, cross-borehole), electrical impedance tomography, and magnetoencephalography show substantial computational gains and practical applicability for uncertainty quantification and proposal generation in Bayesian workflows.

Abstract

It is well-known that the posterior density of linear inverse problems with Gaussian prior and Gaussian likelihood is also Gaussian, hence completely described by its covariance and expectation. Sampling from a Gaussian posterior may be important in the analysis of various non-Gaussian inverse problems in which a estimates from a Gaussian posterior distribution constitute an intermediate stage in a Bayesian workflow. Sampling from a Gaussian distribution is straightforward if the Cholesky factorization of the covariance matrix or its inverse is available, however when the unknown is high dimensional, the computation of the posterior covariance maybe unfeasible. If the linear inverse problem is underdetermined, it is possible to exploit the orthogonality of the fundamental subspaces associated with the coefficient matrix together with the idea behind the Randomize-Then-Optimize approach to design a low complexity posterior sampler that does not require the posterior covariance to be formed. The performance of the proposed sampler is illustrated with a few computed examples, including non-Gaussian problems with non-linear forward model, and hierarchical models comprising a conditionally Gaussian submodel.

Paper Structure

This paper contains 7 sections, 5 theorems, 108 equations, 6 figures, 1 table.

Table of Contents

  1. Introduction
  2. Splitting-based Randomize-Then-Optimze sampling of a Gaussian posterior
  3. Splitting-based sampling and adjoint problem
  4. Matrix-free sampling
  5. Priors and linear transformations
  6. Computed examples
  7. Discussion

Key Result

Theorem 2.1

Assume that ${\mathsf A}\in{\mathbb R}^{m\times n}$, and $b\in{\mathbb R}^m$ is an indirect noisy observation (linsys) of $x\in{\mathbb R}^n$, where and $x$ and $\varepsilon$ are mutually independent. We define where $b^j = b+\eta^j$ and $x_0^j = x_0 + \nu^j$, the random perturbations being drawn from the densities Then the points $x^j$ are independent samples from the posterior density $\pi_{x

Figures (6)

  • Figure 1: The configuration of the fanbeam tomography example. In the figure, three equally spaced illumination views are shown.
  • Figure 2: The generative model for the simulated data is shown in the left panel. The middle and the right panels show the estimated posterior mean and marginal standard deviations computed from a sample of size $10\,000$.
  • Figure 3: Left: The generative structure with one seismic source in the left borehole, and a line of receivers in the right hole. The other three panels show the sample-based estimated posterior mean, and the lower and upper quartiles of 25% and 75%, computed pixel-by-pixel.
  • Figure 4: On the top left, the generative model $\gamma =\log(\sigma/\sigma_0+1)$ for the data is shown, and the bottom left panel shows the corresponding estimate with the linearized model. The four panels in the middle and on the right show four random draws from the prior distribution.
  • Figure 5: The source space consisting of 15 002 dipoles, their locations marked by gray dots. The selected active parcel of dipoles are located on the anterior part of the right circular insular sulcus, are indicated in red, while the active dipoles generating the data with stars. The blue dots indicate the dipole positions on the anterior horizontal part of the lateral fissure, corresponding to the parcel of dipoles having the second largest posterior expectation for the variance.
  • ...and 1 more figures

Theorems & Definitions (5)

  • Theorem 2.1
  • Lemma 2.2
  • Lemma 2.3
  • Theorem 2.4
  • Theorem 2.5