Efficient Derivative-free Bayesian Inference for Large-Scale Inverse Problems

The focus of this work is on efficient derivative-free Bayesian inference approaches for large scale inverse problems, in which the goal is to estimate probability densities for uncertain parameters, given noisy observations derived from the output of a model that depends on the parameters. Such approaches are highly desirable for numerous models arising in science and engineering applications, often defined through partial differential equations (PDEs). These include, to name a few, global climate model calibration [1, 2], material constitutive relation calibration [3, 4, 5], seismic inversion in geophysics [6, 7, 8, 9, 10], and biomechanics inverse problems [11, 12]. Such problems may feature multiple scales, may include chaotic dynamics, or may involve turbulent phenomena; as a result the forward models are typically very expensive to evaluate. Moreover, the forward solvers are often given as a black box (e.g., off-the-shelf solvers [13] or multiphysics systems requiring coupling of different solvers [14, 15]), and may not be differentiable due to the numerical methods used (e.g., embedded boundary method [16, 17] and adaptive mesh refinement [18, 19]) or because of the inherently discontinuous physics (e.g. in fracture [20] or cloud modeling [21, 22]).

Traditional methods for derivative-free Bayesian inference to estimate the posterior distribution include specific instances of the Markov chain Monte Carlo methodology [23, 24, 25, 26] (MCMC), such as random walk Metropolis or the preconditioned Crank-Nicolson (pCN) algorithm [26], and Sequential Monte Carlo methods [27, 28] (SMC), which are in any case often interwoven with MCMC. These methods typically require $\mathcal{O}(10^{4})$ iterations, or more, to reach statistical convergence for the complex forward models which motivate our work. Given that each forward run can be expensive, conducting $\mathcal{O}(10^{4})$ runs is often computationally unfeasible. We present an approach based on the Kalman filter methodology, which aims to estimate the first two moments of the posterior distribution. We demonstrate that, in numerical tests across a range of examples, the proposed methodologies converge within $\mathcal{O}(10)$ iterations, using $\mathcal{O}(10)$ embarrassingly-parallel model evaluations per step, resulting in orders of magnitude reduction in cost over derivative-free MCMC and SMC methods. We also demonstrate favorable performance in comparison with existing Kalman-based Bayesian inversion techniques.

In Subsection 1.2, we outline the Bayesian approach to inverse problems, describing various approaches to sampling, formulated as dynamical systems on probability measures, and introducing our novel mean field approach. In Subsection 1.3, we discuss pathwise stochastic dynamical systems which realize such dynamics at the level of measures, and discuss filtering algorithms which may be applied to them for the purposes of approximate inversion. Subsection 1.4 highlights the novel contributions in this paper, building on the context established in the two preceding subsections. Subsection 1.5 summarizes notational conventions that we adopt throughout.

Inverse problems can be formulated as recovering unknown parameters $\theta\in\mathbb{R}^{N_{\theta}}$ from noisy observation $y\in\mathbb{R}^{N_{y}}$ related through

Here $\mathcal{G}$ denotes a forward model mapping parameters to output observables, and $\eta$ denotes observational noise; for simplicity we will assume known Gaussian statistics: $\eta\sim\mathcal{N}(0,\Sigma_{\eta})$. In the Bayesian perspective, $\theta$ and $y$ are treated as random variables. Given the prior $\rho_{\rm prior}(\theta)$ on $\theta$, the inverse problem can be formulated as finding the posterior $\rho_{\rm post}(\theta)$ on $\theta$ given $y$ [29, 30, 31]:

and $Z(y)$ is the normalization constant

We focus on the case, where the prior $\rho_{\rm prior}$ is (or is approximated as) Gaussian with mean and covariance $r_{0}$ and $\Sigma_{0}$, respectively. Then the posterior $\rho_{\rm post}(\theta)$ can be written as

Since filtering methods are at the heart of our proposed methodology, we provide here a brief summary of a few key concepts. Filtering methods may be deployed to approximate the posterior distribution given by Eq. 2. The inverse problem is first paired with a dynamical system for the parameter [54, 55, 56, 57], leading to a hidden Markov model, to which filtering methods may be applied. In its most basic form, the hidden Markov model takes the form

here $\theta_{n}$ is the unknown state vector, $y_{n+1}$ is the output of the observation model, and $\eta_{n+1}\sim\mathcal{N}(0,\Sigma_{\eta})$ is the observation error at the $n^{th}$ iteration. Any filtering method can be applied to estimate $\theta_{n}$ given observation data $\{y_{\ell}^{\dagger}\}_{\ell=1}$. The Kalman filter [58] can be applied to this setting provided the forward operator $\mathcal{G}$ is linear and the initial state $\theta_{0}$ and the observation errors are Gaussian. The Kalman filter has been extended to nonlinear and non-Gaussian settings in manifold ways, including but not limited to, the extended Kalman filter (EKF, or sometimes ExKF) [59, 60], the ensemble Kalman filters (EnKF)  [61, 62, 63], and the unscented Kalman (UKF) filter [64, 57]. We refer to the extended, ensemble and unscented Kalman filters as approximate Kalman filters to highlight the fact that, outside the linear setting where the Kalman filter [58] is exact, they are all uncontrolled approximations designed on the principle of matching first and second moments.

More precisely, the EnKF uses Monte Carlo sampling to estimate desired means and covariances empirically. Its update step is of the form (7) and can be either deterministic or stochastic. The ensemble adjustment/transform filters are particle approximations of square root filters, a deterministic approach to matching first and second moment information [65]. The unscented Kalman filter uses quadrature, and is also a deterministic method; it may also be viewed as approximating a square root filter. The stochastic EnKF on the other hand compares the data $y$ to model generated data and its update step is intrinsically stochastic, that is, the vector $b$ in (7) itself is random.

All of the filtering methods to estimate $\theta_{n}$ given $\{y_{\ell}^{\dagger}\}_{\ell=1}^{n}$ that we have described so far may be employed in the setting where $y_{\ell}^{\dagger}\equiv y$; repeated exposure of the parameter to the data helps the system to learn the parameter from the data. In order to maintain statistical consistency, an $N$-fold insertion of the same data $y$ requires an appropriate modification of the data likelihood function and the resulting Bayesian inference step becomes

Initializing with $\rho_{0}(\theta)=\rho_{\rm prior}(\theta)$, after $N$ iterations, $\rho_{N}(\theta)$ is equal to the posterior density. The filtering distribution for (12) recovers this exactly if $y_{\ell}^{\dagger}\equiv y$ and if the variance of $\eta$ is rescaled by $N;$ use of ensemble Kalman methods in this setting leads to approximate Bayesian inference, which is intuitively accurate when the posterior is close to Gaussian. We note that the resulting methodology can be viewed as a homotopy method, such as SMC [27] and transport variants [47], which seek to deform the prior into the posterior in one unit time with a finite number of inner steps $N$ – foundational papers introducing ensemble Kalman methods in this context are [54, 55, 66]. Adaptive time-stepping strategies in this context are explored in [68, 69, 70]. Throughout this paper, we will denote the resulting methods as iterative extended Kalman filter (IEKF), iterative ensemble Kalman filter (IEnKF), iterative unscented Kalman filter (IUKF), iterative ensemble adjustment Kalman filter (IEAKF) and iterative ensemble transport Kalman filter (IETKF).

We emphasize that multiple insertions of the same data $y$ without the adjustment (13) of the data likelihood function, and/or over arbitrary numbers of steps, leads to the class of optimization-based Kalman inversion methods: EKI [56], Tikhonov-regularized EKI, termed TEKI [71] and unscented Kalman inversion, UKI [72]; see also [73] for recent adaptive methodologies which are variants on TEKI. These variants of the Kalman filter lead to efficient derivative-free optimization approaches to approximating the maximum likelihood estimator or maximum a posteriori estimator in the asymptotic limit as $n\to\infty$. The purpose of our paper is to develop similar ideas, based on iteration to infinity in $n$, but to tackle the problem of sampling from the posterior $\rho_{\rm post}(\theta)$ rather than the optimization problem. To achieve these we introduce a novel mean-field stochastic dynamical system, generalizing (12) and apply ensemble Kalman methods to it. This leads to Bayesian analogues of EKI and the UKI. To avoid proliferation of nomenclature, we will also refer to these as EKI and UKI relying on the context to determine whether the optimization or Bayesian approach is being adopted; in this paper our focus is entirely on the Bayesian context. We will also use ensemble adjustment and transform filters, denoted as EAKF and ETKF, noting that these two may be applied in either the optimization (using (12)) or Bayesian (using the novel mean-field stochastic dynamical system introduced here) context, but that here we only study the Bayesian problem. The main conclusions of our work are two-fold, concerning the application of Kalman methods to solve the Bayesian inverse problem: that with carefully chosen underlying mean-field dynamical system, such that the prediction and analysis steps approximate Eq. 9 and replicate  Eq. 10, iterating to infinity leads to more efficient and robust methods than the homotopy methods which transport prior to posterior in a finite number of steps; and that deterministic implementations of ensemble Kalman methods, and variants, are superior to stochastic methods.

The methods we propose are exact in the setting of linear $\mathcal{G}$ and Gaussian prior density $\rho_{\rm prior}$; but, for nonlinear $\mathcal{G}$, the Kalman-based filters we employ generally do not converge to the exact posterior distribution, due to the Gaussian ansatz used when deriving the method; negative theoretical results and numerical evidence are reported in [74, 75]. Nonetheless, practical experience demonstrates that the methodology can be effective for problems with distributions close to Gaussian, a situation which arises in many applications.

Finally, we note that we also include comparisons with the ensemble Kalman sampler [75, 76, 77], which we refer to as the EKS, an ensemble based Bayesian inversion method derived from discretizing a mean-field stochastic differential equation and which is also based on iteration to infinity, that is, on the invariance principle of the posterior distribution; and we include comparison with the CBS approach [52] mentioned above, another methodology which also iterates a mean-field dynamical system to infinity to approximate the posterior.

The key idea underlying this work is the development of an efficient derivative-free Bayesian inference approach based on applying Kalman-based filtering methods to a hidden Markov model arising from a novel mean-field dynamical system. Stemming from this, our main contributions are as follows: ¹¹1In making these statements, we acknowledge that for linear Gaussian problems it is possible to solve the Bayesian inverse problem exactly in one step, or multiple steps, using the Kalman filter in transport/coupling mode, when initialized correctly and with a large enough ensemble. However, the transport/coupling methods are not robust to perturbations from initialization, non-Gaussianity and so forth, whereas the methods we introduce are. Our results substantiate this claim.

1. 1.

   In the setting of linear Gaussian inverse problems, we prove that the filtering distribution of the mean field model converges exponentially fast to the posterior distribution.

2. 2.

   We generalize the inversion methods EKI, UKI, EAKI and ETKI from the optimization to the Bayesian context by applying the relevant variants on Kalman methodologies to the novel mean-field dynamical system (Bayesian) rather than to (12) (optimization).

3. 3.

   We study and compare application of both deterministic and stochastic Kalman methods to the novel mean-field dynamical system, demonstrating that the deterministic methods (UKI, EAKI and ETKI) outperform the stochastic method (EKI); this may be attributed to smooth, noise-free approximations resulting from deterministic approaches.

4. 4.

   We demonstrate that the application of Kalman methods to the novel mean-field dynamical system outperforms the application of Kalman filters to transport/coupling models – the IEnKF, IUKF, IEAKF and IETKF approaches; this may be attributed to the exponential convergence underlying the filter for the novel mean-field dynamical system.

5. 5.

   We also demonstrate that the application of Kalman methods to the novel mean-field dynamical system outperforms the EKS, when Euler-Maruyama discretization is used, because the continuous-time formulation requires very small time-steps, and CBS which suffers from stochasticity, similarly to the EKI.

6. 6.

   We propose several strategies, including low-rank approximation and a bi-fidelity approach, to reduce the computational and memory cost.

7. 7.

   We demonstrate, on both linear and nonlinear model problems (including inference for subsurface geophysical properties in porous medium flow), that application of deterministic Kalman methods to approximate the filtering distribution of the novel mean-field dynamical system delivers mean and covariance which are close to the truth or to those obtained with the pCN MCMC method. The latter uses $\mathcal{O}(10^{4})$ model evaluations or more whilst for our method only $\mathcal{O}(10)$ iterations are required with $\mathcal{O}(10)$ ensemble members, leading to only $\mathcal{O}(10^{2})$ model evaluations, two orders of magnitude savings.

8. 8.

   The method is applied to perform Bayesian parameter inference of subgrid-scale parameters arising in an idealized global climate model, a problem currently far beyond the reach of state-of-the-art MCMC methods such as pCN and variants.

The remainder of the paper is organized as follows. In Section 2, the mean field dynamical system, various algorithms which approximate its filtering distribution, and a complete analysis in the linear setting, are all presented. These correspond to our contributions (i) and (ii). In Section 3, strategies to speed up the algorithm and improve the robustness for real-world problems are presented. These correspond to our contribution (vi). Numerical experiments are provided in Section 4; these serve to empirically confirm the theory and demonstrate the effectiveness of the framework for Bayesian inference. These correspond to our contributions (iii), (iv), (v), (vii) and (viii). We make concluding remarks in Section 5.

The code is accessible online:

https://github.com/Zhengyu-Huang/InverseProblems.jl

$A\succ B$ and $A\succeq B$ denote $A-B$ positive-definite or positive-semidefinite, for symmetric matrices $A,B.$ $\|\cdot\|,\langle\cdot,\cdot\rangle$ denote Euclidean norm and inner-product. We use $\mathbb{Z}^{+}=\{0,1,2,\cdots\}$ to denote the set of natural numbers; $\mathcal{N}(\cdot,\cdot)$ to denote Gaussian distributions; and $\varrho(\cdot)$ to denote the spectral radius. As encountered in Subsection 1.2 we make use of the similar symbol $\rho$ for densities; these should be easily distinguished from spectral radius by context and by a different font.

Following the discussion from Section 1.2.2, we propose an implementation of (8) to solve inverse problems by pairing the parameter-to-data map with a dynamical system for the parameter, and then employ techniques from filtering to estimate the parameter given the data.

We introduce the prediction step

Here $\theta_{n+1}$ is the unknown state vector and $\omega_{n+1}\sim\mathcal{N}(0,\Sigma_{\omega,n+1})$ is the independent, zero-mean Gaussian evolution error, which will be chosen such that (14) mimics (9) for Gaussian densities. The analysis step (10) follows exactly from the observation model (first introduced in [71])

here we have defined the augmented forward map

with $\nu_{n+1}\sim\mathcal{N}(0,\Sigma_{\nu,n+1})$ the independent, zero-mean Gaussian observation error and $x_{n+1}$ the output of the observation model at time $n+1$. We define artificial observation $x_{n+1}^{\dagger}$ using the following particular instance of the data, constructed from the one observation $y$ and the prior mean $r_{0}$ and assumed to hold for all $n\geq 1:$

We will apply filtering methods to condition $\theta_{n}$ on $Y_{n}:=\{x_{1}^{\dagger},x_{2}^{\dagger},\cdots,x_{n}^{\dagger}\}$, the observation set at time $n$. As we will see later, the choice of $\{x_{\ell}^{\dagger}\}_{l=1}$ leads to the correct posterior.

Let $C_{n}$ denote the covariance of the conditional random variable $\theta_{n}|Y_{n}$. Then the error covariance matrices $\{\Sigma_{\omega,n+1}\}$ and $\{\Sigma_{\nu,n+1}\}$ in the extended dynamical system (14)-(15) are chosen at the $n^{th}$ iteration, as follows:

Here ${\color[rgb]{0,0,0}{0<\Delta\tau<1}}$, and in our numerical studies we choose ${\color[rgb]{0,0,0}{\Delta\tau=1/2}}$, although other choices are possible. Since the artificial evolution error covariance $\Sigma_{\omega,n+1}$ in (14) is updated based on $C_{n}$, the conditional covariance of $\theta_{n}|Y_{n}$, it follows that (14) is a mean-field dynamical system: it depends on its own law, specifically on the law of $\theta_{n}|Y_{n}$. Details underpinning the choices of the error covariance matrices $\Big{\{}\Sigma_{\omega,n}\Big{\}}$ and $\Big{\{}\Sigma_{\nu,n}\Big{\}}$ are given in Subsections 2.2 and 2.3: the matrices are chosen so that, for linear Gaussian problems, the prediction and analysis steps follow Eqs. 9 and 10, and the converged mean and covariance of the resulting filtering distribution for $\theta_{n}|Y_{n}$ under the prediction step (14) and the observation model (15) match the posterior mean and covariance.

Denote by $\rho_{n}$, the conditional density of $\theta_{n}|Y_{n}$. We first introduce a class of conceptual Kalman inversion algorithms which approximate $\rho_{n}$ by considering only first and second order statistics (mean and covariance), and update $\rho_{n}$ sequentially using the standard prediction and analysis steps [45, 46]: $\rho_{n}\mapsto\hat{\rho}_{n+1}$, and then $\hat{\rho}_{n+1}\mapsto\rho_{n+1}$, where $\hat{\rho}_{n+1}$ is the distribution of $\theta_{n+1}|Y_{n}$. The second analysis step is performed by invoking a Gaussian hypothesis. In subsequent subsections, we then apply different methods to approximate the resulting maps on measures, leading to unscented, stochastic ensemble Kalman and adjustment/transform square root Kalman filters.

In the prediction step, assume that $\rho_{n}\approx\mathcal{N}(m_{n},C_{n})$, then under Eq. 14, $\hat{\rho}_{n+1}=\mathcal{N}(\widehat{m}_{n+1},\widehat{C}_{n+1})$ is also Gaussian and satisfies

In the analysis step, we assume that the joint distribution of $\{\theta_{n+1},x_{n+1}\}|Y_{n}$ can be approximated by a Gaussian distribution

where

These expectations are computed by assuming $\theta_{n+1}|Y_{n}\sim\hat{\rho}_{n+1}$ and noting that the distribution of $(\theta_{n+1},x_{n+1})$ is then defined by (14)-(15). This corresponds to projecting ²²2We use the term “projecting” as finding the Gaussian $p$ which matches the first and second moments of a given measure $\pi$ corresponds to finding the closest Gaussian $p$ to $\pi$ with respect to variation in the second argument of the (nonsymmetric) Kullback-Leibler divergence [78][Theorem 4.5]. the joint distribution onto the Gaussian which matches its mean and covariance. Conditioning the Gaussian in Eq. 20 to find $\theta_{n+1}|\{Y_{n},x_{n+1}^{\dagger}\}=\theta_{n+1}|Y_{n+1}$, gives the following expressions for the mean $m_{n+1}$ and covariance $C_{n+1}$ of the approximation to $\rho_{n+1}$ :

Equations 19, 20, 21, 22a and 22b establish a class of conceptual algorithms for application of Gaussian approximation to solve the inverse problems. To make implementable algorithms a high level choice needs to be made: whether to work strictly within the class of Gaussians, that is to impose $\rho_{n}\equiv\mathcal{N}(m_{n},C_{n})$, or whether to allow non-Gaussian $\rho_{n}$ but to insist that the second order statistics of the resulting measures agree with Equations 19, 20, 21, 22a and 22b. In what follows the UKI takes the first perspective; all other methods take the second perspective. For the UKI the method views Equations 19, 20, 21, 22a and 22b as providing a nonlinear map $(m_{n},C_{n})\mapsto(m_{n+1},C_{n+1})$; this map is then approximated using quadrature. For the remaining methods a mean-field dynamical system is used, which is non-Gaussian but matches the aforementioned Gaussian statistics; this mean-field model is then approximated by a finite particle system [79]. The dynamical system is of mean-field type because of the expectations required to calculate Equations 20 and 21 and (18). The continuous time limit of the evolution for the mean and covariance is presented in A; this is obtained by letting $\Delta\tau\to 0.$

In this subsection, we study the algorithm in the context of linear inverse problems, for which $\mathcal{G}(\theta)=G\theta$ for some matrix $G.$ Furthermore we assume that $\rho_{\rm prior}$ is Gaussian $\mathcal{N}(r_{0},\Sigma_{0})$ and recall that the observational noise is $\mathcal{N}(0,\Sigma_{\eta}).$ Thanks to the linear Gaussian structure the posterior is also Gaussian with mean and precisions given by

Furthermore, the equations (21) reduce to

We note that

The update equations (22a)(22b) become

with $\widehat{C}_{n+1}=C_{n}+\Sigma_{\omega,n+1}$. We have the following theorem about the convergence of the algorithm:

To make practical methods for solving nonlinear inverse problems (1) out of the foregoing, the expectations (integrals) appearing in the prediction step (19) as well as in the analysis step via Eq. 21 need to be approximated appropriately. While Eq. 19 can be implemented via a simple rescaling of the covariance matrix or ensemble, respectively, (we use both) the analysis step can be implemented using any nonlinear Kalman filter (we use a variety). In the present work, we focus on both the unscented and ensemble Kalman filters, which lead to the Bayesian implementations of unscented Kalman inversion (UKI), stochastic ensemble Kalman inversion (EKI), ensemble adjustment Kalman inversion (EAKI), and ensemble transform Kalman inversion (ETKI). We now detail these methods. ³³3Recall the discussion in Subsection 1.3 about distinction between optimization and Bayesian implementations of all these methods.