Fast Krylov Subspace Methods For Geostatistical Inverse Problems

Fast Krylov Subspace Methods for Geostatistical Inverse Problems

Inverse problems are part of a mathematical framework used to estimate parameters that characterize a physical system but are difficult to measure directly. An example application is hydraulic tomography which is a method of imaging the subsurface. Water is pumped at designated pumping wells and the measured pressure response is recorded at corresponding observation wells. These noisy measurements of pressure are then used to obtain tomographic reconstructions of important hydrogeological parameters such as hydraulic conductivity and storage. Inverse problems of this type are particularly challenging both mathematically as they are ill-posed, and computationally as they generally require repeated solutions of large-scale partial differential equations. In this dissertation, we discuss in detail a method of imaging based on oscillatory pumping tests. We discuss methods to extract the signal from the noise and examine the duration of the transient. For this class of inverse problems, solving shifted systems of linear equations is a major computational bottleneck. In the dissertation, we describe an iterative algorithm for solving shifted linear systems. Krylov subspace methods are particularly appealing because of their shift-invariant property. By exploiting this property, only a single Krylov basis is computed and the solution for multiple shifts can be obtained at a cost that is nearly equal to the cost of solving a single system. We then show how the time dependent inverse problems can be accelerated using these Krylov subspace solvers using a Laplace-transform approach.

SIAM Journal on Scientific Computing

Fast Algorithms for Geostatistical Inverse Problems and Uncertainity Quantification

Inverse problems are ubiquitous in science and engineering, especially in the field of geosciences. Our work seeks to develop and implement computationally efficient statistical methods for optimizing the use of limited (and possibly noisy) environmental data to accurately determine (or image) heterogeneous subsurface geological properties, and quantify the corresponding predictive uncertainty to provide a sound basis for management or policy decision making. This work tackles inverse problems using the Geostatistical approach that stochastically models unknowns as random fields. However, a direct implementation of the Geostatistical approach is challenging due to high computational costs in identifying small scale features. The costs occur because solving inverse problems requires several expensive simulations of partial differential equations as well as representing high dimensional random fields on irregular grids and complicated domains. Our approach uses Hierarchical matrices to efficiently represent dense covariance matrices and solves the resulting intermediary system of equations to compute the maximum a posteriori (MAP) estimate using preconditioned iterative methods. The resulting cost is reduced from $O(N^2)$ to $O(N\log N)$, where N is the number of unknowns to be determined. Uncertainty quantification in the geostatistical approach can be performed by computing the posterior covariance at the MAP estimate. We derive efficient representation of the posterior covariance matrix at the MAP point as the sum of the prior covariance matrix and a low-rank update that contains information from the dominant generalized eigenmodes of the data misfit part of the Hessian and the inverse covariance matrix. The rank of the low-rank update is typically independent of the dimension of the unknown parameter. We provide an efficient randomized algorithm for computing the dominant eigenmodes of the generalized eigenvalue problem (and as a result, the low-rank decomposition) that avoids forming square-roots of the covariance matrix or its inverse. As a result, we have a method that scales almost linearly with the dimension of unknown parameter space and the data dimension. Further, we show how to efficiently compute some measures of uncertainty that are based on scalar invariants of the posterior covariance matrix. The resulting uncertainty measures can be used in the context of optimal experimental design. Finally, we consider a specific application, namely oscillatory hydraulic tomography (OHT). In OHT, the reconstruction of hydrogeological parameters, such as hydraulic conductivity and specific storage, using limited discrete measurements of pressure (head) obtained from sequential oscillatory pumping tests, leads to a nonlinear inverse problem. The adjoint approach for computing the Jacobian requires repeated solution of the forward (and adjoint) problem for multiple frequencies which requires solutions of a shifted system of linear equations. We develop flexible preconditioned Krylov subspace solvers specifically designed for shifted systems. We analyze the convergence of the solver and when an iterative solver is used for inverting the preconditioner matrices.