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

TL;DR

This paper introduces a low-complexity sampling method for Gaussian posteriors in high-dimensional linear inverse problems, leveraging subspace splitting and orthogonality, applicable even when covariance matrices are infeasible to compute.

Contribution

It proposes a novel subspace splitting approach that enables efficient sampling from Gaussian posteriors without forming the covariance matrix, especially useful in high-dimensional or underdetermined problems.

Findings

Efficient sampling demonstrated in high-dimensional settings.

Applicable to non-Gaussian and hierarchical models.

Avoids covariance matrix computation in large-scale problems.

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.