Efficient sampling of high-dimensional Gaussian fields: the non-stationary / non-sparse case

This paper is devoted to the problem of sampling Gaussian fields in high dimension. Solutions exist for two specific structures of inverse covariance : sparse and circulant. The proposed approach is valid in a more general case and especially as it emerges in inverse problems. It relies on a perturbation-optimization principle: adequate stochastic perturbation of a criterion and optimization of the perturbed criterion. It is shown that the criterion minimizer is a sample of the target density. The motivation in inverse problems is related to general (non-convolutive) linear observation models and their resolution in a Bayesian framework implemented through sampling algorithms when existing samplers are not feasible. It finds a direct application in myopic and/or unsupervised inversion as well as in some non-Gaussian inversion. An illustration focused on hyperparameter estimation for super-resolution problems assesses the effectiveness of the proposed approach.

💡 Research Summary

The paper addresses the challenging problem of sampling high‑dimensional Gaussian random fields when the precision matrix does not possess the special structures (sparsity or circulant) that enable existing fast samplers. Classical approaches work well only when the precision matrix Q is either sparse—allowing block Gibbs updates or Cholesky‑based direct solves—or circulant—allowing independent sampling in the Fourier domain followed by an inverse FFT. In many Bayesian inverse‑problem settings, however, the precision matrix takes the form

  Q = Σₖ Mₖᵀ Rₖ⁻¹ Mₖ,

where each Rₖ is a relatively small, easily invertible covariance (e.g., noise or prior covariances) and each Mₖ is a linear operator such as the forward model H or the identity. This Q is generally dense and non‑circulant, making the standard methods computationally infeasible.

The authors propose a novel “Perturbation‑Optimization” (PO) algorithm that generates an exact sample from N(0, Q⁻¹) without ever forming Q⁻¹ explicitly. The algorithm consists of two steps:

1. Perturbation (P): Independently draw ηₖ ∼ N(0, Rₖ) for k = 1,…,K. Because each Rₖ is simple, sampling is cheap.
2. Optimization (O): Minimize the quadratic criterion

  J(x | η₁,…,η_K) = Σₖ (ηₖ − Mₖ x)ᵀ Rₖ⁻¹ (ηₖ − Mₖ x).

The minimizer satisfies the linear system

  Q x = Σₖ Mₖᵀ Rₖ⁻¹ ηₖ,

so the solution is

  x̂ = Q⁻¹ Σₖ Mₖᵀ Rₖ⁻¹ ηₖ.

A rigorous proof shows that x̂ is a zero‑mean Gaussian vector with covariance Q⁻¹, i.e., it is an exact draw from the target distribution. The authors also extend the result to non‑zero means by perturbing with ζₖ ∼ N(mₖ, Rₖ), leading to a sample from N(Q⁻¹ Σₖ Mₖᵀ Rₖ⁻¹ mₖ, Q⁻¹).

The key computational advantage is that the algorithm never requires an explicit factorization of Q. Instead, one needs to (i) generate the ηₖ, (ii) compute the right‑hand side Σₖ Mₖᵀ Rₖ⁻¹ ηₖ, and (iii) solve the linear system Q x = rhs. The linear solve can be performed with any iterative method (e.g., conjugate gradient) that only needs matrix‑vector products with Q, which are cheap because Q is a sum of terms Mₖᵀ Rₖ⁻¹ Mₖ. In many imaging inverse problems, Mₖ are convolution, down‑sampling, or identity operators, and Rₖ⁻¹ are diagonal or block‑diagonal, so each product can be implemented efficiently.

The paper then embeds the PO sampler into a Bayesian Gibbs framework for linear inverse problems. With a forward model y = H x + n, Gaussian noise n ∼ N(m_n, R_n) and Gaussian prior x ∼ N(m_x, R_x) (both possibly conditioned on hyper‑parameters θ), the posterior p(x | θ, y) is Gaussian with precision Q = Hᵀ R_n⁻¹ H + R_x⁻¹. Since H is typically non‑circulant (e.g., blur followed by decimation) and dense, standard samplers fail. By setting K = 2, M₁ = H, M₂ = I, R₁ = R_n, R₂ = R_x, the PO algorithm yields an exact draw from the posterior at each Gibbs iteration, enabling joint estimation of x and hyper‑parameters θ.

To demonstrate practicality, the authors apply the method to a super‑resolution (SR) problem. Low‑resolution observations are modeled as y = P C x + n, where C is a convolution (blur) and P is a down‑sampling operator. The prior on x is a Gaussian with precision γ_x Dᵀ D (D is a Laplacian). Hyper‑parameters γ_n and γ_x have Jeffreys’ priors, leading to Gamma conditional posteriors. A Gibbs sampler alternates: (i) sampling γ_n and γ_x from their Gamma conditionals, (ii) sampling x from its Gaussian posterior using the PO algorithm (K = 2). The experiments show rapid convergence of the hyper‑parameter chains, accurate reconstruction of the high‑resolution image, and credible intervals that contain the ground truth, illustrating the ability to quantify uncertainty.

The paper also mentions two additional applications: (a) electromagnetic inverse scattering, where the forward model is bilinear in object and induced current, leading to conditional Gaussian structures; (b) fluorescence microscopy with structured illumination, where non‑stationary illumination yields a non‑invariant covariance. In both cases, the PO sampler enables Gibbs‑based inference that would otherwise be impossible.

In conclusion, the Perturbation‑Optimization framework provides a general, exact, and computationally feasible way to sample high‑dimensional Gaussian fields whenever the precision matrix can be expressed as a sum of low‑rank or easily invertible terms. Its reliance on simple Gaussian perturbations and iterative linear solvers makes it attractive for a wide range of Bayesian inverse problems. Future work could explore extensions to non‑linear operators, preconditioning strategies to accelerate the linear solves, and distributed implementations for truly massive datasets.