A fast algorithm for parabolic PDE-based inverse problems based on Laplace transforms and flexible Krylov solvers

Abstract

We consider the problem of estimating parameters in large-scale weakly nonlinear inverse problems for which the underlying governing equations is a linear, time-dependent, parabolic partial differential equation. A major challenge in solving these inverse problems using Newton-type methods is the computational cost associated with solving the forward problem and with repeated construction of the Jacobian, which represents the sensitivity of the measurements to the unknown parameters. Forming the Jacobian can be prohibitively expensive because it requires repeated solutions of the forward and adjoint time-dependent parabolic partial differential equations corresponding to multiple sources and receivers. We propose an efficient method based on a Laplace transform-based exponential time integrator combined with a flexible Krylov subspace approach to solve the resulting shifted systems of equations efficiently. Our proposed solver speeds up the computation of the forward and adjoint problems, thus yielding significant speedup in total inversion time. We consider an application from Transient Hydraulic Tomography (THT), which is an imaging technique to estimate hydraulic parameters related to the subsurface from pressure measurements obtained by a series of pumping tests. The algorithms discussed are applied to a synthetic example taken from THT to demonstrate the resulting computational gains of this proposed method.

Multipreconditioned Gmres for Shifted Systems

“…Matrix Functions. Matrix function evaluations are relevant in many applications; for example, the evaluation of exp(−tA)b is important in the time-integration of large-scale dynamical systems [2]. In the field of statistics and uncertainty quantification, several computations involving a symmetric positive definite covariance matrix A can be expressed in terms of matrix functions.…”

“…The next set of basis vectors V (2) , with orthonormal columns, are generated by orthogonalizing the columns of AZ (1) against V (1) = r 0 /β, and performing a thin QR factorization. The orthogonalization coefficients generated in this process are stored in matrices H (j,1) , j = 1, 2.…”

“…The orthogonalization coefficients generated in this process are stored in matrices H (j,1) , j = 1, 2. The space of search directions is increased by applying each preconditioner to this new matrix; i.e., Z (2) = [P −1 1 V (2) . .…”

“…. P −1 np V (2) ] ∈ R n×n 2 p . The general procedure at step k is to block orthogonalize AZ (k) against the matrices V (1) , .…”

F-Operators for the Construction of Closed Form Solutions to Linear Homogenous PDEs with Variable Coefficients

“…Local fractional Fourier transform operator via Mittag-Leffler function defined on the fractal set is used to solve PDEs in [2]. An efficient computational scheme based on a Laplace transform-based exponential time integrator combined with a flexible Krylov subspace approach is proposed in [3] to solve linear, time-dependent, parabolic PDEs. Fourier wavelets are used to construct solutions to partial differential equations in [4].…”