Efficient Krylov subspace methods for uncertainty quantification in large Bayesian linear inverse problems

Abstract

Summary Uncertainty quantification for linear inverse problems remains a challenging task, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. This work exploits Krylov subspace methods to develop and analyze new techniques for large‐scale uncertainty quantification in inverse problems. In this work, we assume that generalized Golub‐Kahan‐based methods have been used to compute an estimate of the solution, and we describe efficient methods to explore the posterior distribution. In particular, we use the generalized Golub‐Kahan bidiagonalization to derive an approximation of the posterior covariance matrix, and we provide theoretical results that quantify the accuracy of the approximate posterior covariance matrix and of the resulting posterior distribution. Then, we describe efficient methods that use the approximation to compute measures of uncertainty, including the Kullback‐Liebler divergence. We present two methods that use the preconditioned Lanczos algorithm to efficiently generate samples from the posterior distribution. Numerical examples from dynamic photoacoustic tomography demonstrate the effectiveness of the described approaches.

Computational tools for inversion and uncertainty estimation in respirometry

“…Since the posterior density function is Gaussian, variance estimates for the solution can be obtained by computing the diagonal entries of Γ post . Furthermore, samples from the posterior can be obtained using efficient Krylov subspace methods [ 13 , 14 ].…”

Computationally efficient methods for large-scale atmospheric inverse modeling

“…However, neither computing nor storing the covariance matrix is feasible, making further uncertainty estimation challenging. Instead, we follow the approach described in (Chung et al, 2018;Saibaba et al, 2020) for the fixed mean case, where an approximation to the posterior covariance matrix is obtained using the computed vectors generated using the generalized Golub-Kahan bidiagonalization process. An advantage of this approach is that, by storing partial information while computing the MAP estimate, we can approximately compute the uncertainty associated with the MAP estimate (e.g., posterior variance) with minimal additional cost, and no further accesses to the forward and adjoint models.…”

“…More precisely, in addition to storing the information required for storing Q, we only need to store nk +k additional entries corresponding to the matrices Z K k and ∆ K k . Furthermore, the error in the low-rank approximation can be analyzed using similar techniques as in (Saibaba et al, 2020). Similar to the approach described in Section 3.2, the posterior variance, which corresponds to the diagonal entries of Γ post , can be approximated using the diagonal entries of (B6).…”

“…The purpose of this study is to integrate several state-of-the-art computational and mathematical tools with AIM-tools that have been developed for and have had considerable success in other scientific fields (e.g., passive seismic tomography, medical imaging). Specifically, we investigate the use of generalized hybrid (genHyBR) projection methods for surface flux estimation and extend their use for inverse problems where the mean of the fluxes is not known a priori (sometimes referred to as geostatistical inverse modeling) (Chung and Saibaba, 2017;Saibaba et al, 2020). We address these challenges in our present work.…”

“…Discussion started: 20 January 2022 c Author(s) 2022. CC BY 4.0 License.Note that in previous works(Chung et al, 2018;Saibaba et al, 2020), we found that additional reorthogonalization of the generalized Golub-Kahan basis vectors yielded more accurate results, so we perform them in the numerical experiments.Algorithm 2 Compute the diagonals of the low-rank term YZ . Call as [d] =LowRankDiag(Y, Z) Require: Matrices Y, Z ∈ R n×k defining the outer product YZ 1: for i = 1, .…”