Data-informed posterior approximation for Bayesian linear inverse problems

Technical Summary: Data-informed posterior approximation for Bayesian linear inverse problems

Problem Statement
The paper addresses the computational challenges inherent in large-scale Bayesian linear inverse problems, formulated as $y = Gx + \eta$, where $x \in \mathbb{R}^n$ is the unknown parameter, $y \in \mathbb{R}^m$ is the observation, and $\eta$ is Gaussian noise. The primary difficulty arises when the parameter dimension $n$ is extremely large ($n \gg m$), making the direct computation of the posterior distribution's mean and covariance intractable. Specifically, the inversion of the matrix $G^\top \Gamma^{-1} G + \lambda \Sigma^{-1}$ (where $\Sigma$ is the prior covariance and $\lambda$ is a hyperparameter) is computationally prohibitive. Existing dimension reduction techniques, such as those based on the Karhunen–Loève expansion or likelihood-informed subspaces (LIS), often rely on explicit factorizations of the prior covariance or operate entirely within the high-dimensional parameter space, failing to fully exploit the interaction between the forward operator and the observation data. Furthermore, the selection of the hyperparameter $\lambda$ typically requires repeated large-scale computations in the parameter space.

Methodology
The authors propose a framework that shifts the computational focus from the parameter space to the data space. The core of the methodology involves the following components:

1. Data-Space Geometric Framework:
   The paper defines a quotient Hilbert space $(\mathbb{R}^m/\mathcal{N}(M), \langle \cdot, \cdot \rangle_M)$, where $M = G\Sigma G^\top$ and $\mathcal{N}(M)$ is its null space. This space is endowed with an inner product $\langle [u], [v] \rangle_M = u^\top M v$. The authors establish an isometric embedding $\iota$ from this data space into the parameter space $(\mathbb{R}^n, \langle \cdot, \cdot \rangle_{\Sigma^{-1}})$. They prove that the prior-to-posterior update is confined entirely to the image of this embedding, $\text{Im}(\iota) = \mathcal{R}(\Sigma G^\top)$, which has a dimension equal to $\text{rank}(G)$. This reveals that the informative directions of the posterior are intrinsically low-dimensional and independent of the ambient parameter dimension $n$.

2. Quotient-Space Golub–Kahan Bidiagonalization (Q-GKB):
   To construct data-informed Krylov subspaces without explicitly forming $M$ or computing orthogonal projections onto its range, the authors develop a novel iterative algorithm called Q-GKB. This method applies the Golub–Kahan bidiagonalization process directly to the operators defined on the quotient space. The algorithm generates two sequences of vectors, $\{u_i\}$ and $\{v_i\}$, which form orthonormal bases for paired Krylov subspaces in the data space. Crucially, Q-GKB is a matrix-free method that relies only on matrix-vector products with $M$ and $\Gamma^{-1}$, avoiding the storage and factorization of large dense matrices.

3. Integrated Empirical Bayesian Inference:
   The framework integrates empirical Bayesian inference to simultaneously estimate the hyperparameter $\lambda$ and approximate the posterior. By projecting the marginal likelihood optimization problem onto the low-dimensional Krylov subspaces generated by Q-GKB, the method computes an approximate negative log-marginal likelihood $L^{(k)}(\lambda)$ using only small-scale matrices derived from the iteration. This allows for the efficient determination of $\lambda_k$ at each step without returning to the high-dimensional parameter space.

4. Posterior Approximation:
   The approximate posterior mean and covariance are constructed within the data-informed Krylov subspace. The approximate mean coincides with the iteratively regularized solution, and the approximate covariance is computed via the inversion of small-scale matrices. The paper provides theoretical error bounds for these approximations using the Förstner distance and Kullback–Leibler (KL) divergence, showing that the error is governed by the decay of the generalized eigenvalues of the matrix pair $(M, \Gamma)$.

Key Contributions

• Geometric Characterization: The paper provides a rigorous mathematical characterization of the low-dimensional nature of posterior updates through the quotient-space geometry of the data space, establishing an isometric correspondence between data and parameter spaces.
• Q-GKB Algorithm: A new quotient-space formulation of the Golub–Kahan bidiagonalization process is developed, enabling the construction of data-informed Krylov subspaces in a fully matrix-free manner.
• Simultaneous Estimation: The integration of empirical Bayesian inference into the Q-GKB iteration enables simultaneous hyperparameter estimation and posterior approximation, eliminating the need for repeated high-dimensional computations.
• Theoretical Guarantees: The authors derive explicit error bounds for the approximate posterior mean and covariance, linking the convergence rate to the spectral properties of the data space.

Numerical Results
The method is validated on three representative problems:

1. One-dimensional Fredholm integral equation: Demonstrates rapid convergence of the hyperparameter and solution reconstruction, with the effective data-space dimension being significantly smaller than the parameter dimension.
2. Two-dimensional image deblurring: Shows the method's ability to handle moderately ill-posed problems. The Q-GKB approach converges faster than the direct LIS method for the posterior mean and achieves comparable accuracy for the covariance, while being computationally feasible on standard hardware where LIS is not.
3. X-ray computed tomography (CT): A large-scale application ($n=65,536$) where the prior covariance is handled via a matrix-free FFT-based approach. The results show that the method captures the posterior mean and uncertainty (variance) effectively, with the posterior mean converging rapidly and the covariance stabilizing as the subspace dimension increases.

Significance and Claims
The paper claims that the data-space formulation offers a fundamentally different computational paradigm for Bayesian inverse problems. By shifting the focus to the data space, the method circumvents the computational bottlenecks associated with high-dimensional parameter spaces and explicit prior covariance factorizations. The authors assert that their framework provides a rigorous geometric understanding of data-informed subspaces and offers a scalable, matrix-free algorithm for simultaneous hyperparameter estimation and posterior approximation. The numerical experiments support the theoretical findings, demonstrating that the method is effective, stable, and capable of handling large-scale problems where direct access to the posterior distribution is computationally prohibitive. The paper concludes that this perspective opens new paths for posterior inference and suggests future directions such as extending the framework to infinite-dimensional and nonlinear inverse problems.