Approximate Bayesian Model Inversion for PDEs with Heterogeneous and State-Dependent Coefficients

TL;DR

This paper introduces two approximate Bayesian inference methods for PDE parameter estimation with space- and state-dependent coefficients, offering accurate and computationally efficient alternatives to traditional MCMC techniques.

Contribution

The paper develops two novel variational Bayesian methods, Laplace-EM and DSVI-EB, for efficient parameter inference in PDE models with complex coefficients, improving computational cost and accuracy.

Findings

Both methods accurately estimate posterior densities and hyperparameters.

Laplace-EM provides higher accuracy but requires Hessian computations.

DSVI-EB is less accurate but only needs gradient information.

Abstract

We present two approximate Bayesian inference methods for parameter estimation in partial differential equation (PDE) models with space-dependent and state-dependent parameters. We demonstrate that these methods provide accurate and cost-effective alternatives to Markov Chain Monte Carlo simulation. We assume a parameterized Gaussian prior on the unknown functions, and approximate the posterior density by a parameterized multivariate Gaussian density. The parameters of the prior and posterior are estimated from sparse observations of the PDE model's states and the unknown functions themselves by maximizing the evidence lower bound (ELBO), a lower bound on the log marginal likelihood of the observations. The first method, Laplace-EM, employs the expectation maximization algorithm to maximize the ELBO, with a Laplace approximation of the posterior on the E-step, and minimization of a Kullback-Leibler divergence on the M-step. The second method, DSVI-EB, employs the doubly stochastic variational inference (DSVI) algorithm, in which the ELBO is maximized via gradient-based stochastic optimization, with nosiy gradients computed via simple Monte Carlo sampling and Gaussian backpropagation. We apply these methods to identifying diffusion coefficients in linear and nonlinear diffusion equations, and we find that both methods provide accurate estimates of posterior densities and the hyperparameters of Gaussian priors. While the Laplace-EM method is more accurate, it requires computing Hessians of the physics model. The DSVI-EB method is found to be less accurate but only requires gradients of the physics model.

Tables

Table 1. Table 1 : Reference and estimated hyperparameters, and simple MC estimate of the ELBO, for the one-dimensional linear-diffusion problem

			Hyperparameters
		$\widehat{ℱ}$	$\sigma$	$\lambda$
Reference			$1.000$	$0.150$
Laplace-EM		$-37.37\left(5\right)$	$0.608$	$0.132$
DSVI	Full rank	$-43.38\left(8\right)$	$0.551$	$0.120$
	Chevron $k=20$	$-49.54\left(9\right)$	$0.653$	$0.144$
	Chevron $k=10$	$-50.47\left(9\right)$	$0.672$	$0.147$
	Chevron $k=5$	$-47.85\left(6\right)$	$0.681$	$0.148$
	Mean field	$-49.34\left(6\right)$	$0.687$	$0.151$

Table 2. Table 2 : Reference and estimated hyperparameters, and simple MC estimate of the ELBO, for the one-dimensional nonlinear-diffusion problem

			Hyperparameters
		$\widehat{ℱ}$	$\sigma$	$\lambda$
Laplace-EM		$-12.68\left(3\right)$	$4.650$	$6.893$
DSVI	Full rank	$-13.97\left(4\right)$	$4.024$	$5.968$
	Chevron $k=10$	$-14.56\left(4\right)$	$4.007$	$6.265$
	Chevron $k=5$	$-15.02\left(4\right)$	$4.343$	$6.778$
	Chevron $k=2$	$-16.23\left(4\right)$	$4.768$	$7.516$
	Mean field	$-16.66\left(4\right)$	$5.353$	$9.297$

Equations

$${\mathbf{u}}_{\mathrm{s}}={\mathbf{H}}_u\mathbf{u}+{\mathbf{ϵ}}_u,{\mathbf{ϵ}}_{us}\sim \mathcal{N}(0,{\sigma }_{us}^2{\mathbf{I}}_{M_s}),$$

$${\mathbf{y}}_{\mathrm{s}}={\mathbf{H}}_y\mathbf{y}+{\mathbf{ϵ}}_y,{\mathbf{ϵ}}_{ys}\sim \mathcal{N}(0,{\sigma }_{ys}^2{\mathbf{I}}_{N_s}),$$

$$\log p({\mathcal{D}}_s\mid \mathbf{y})\equiv -\frac{1}{2{\sigma }_{us}^2}\Vert {\mathbf{u}}_s-{\mathbf{H}}_u\mathbf{u}{\Vert }_2^2-\frac{1}{2{\sigma }_{ys}^2}\Vert {\mathbf{y}}_s-{\mathbf{H}}_y\mathbf{y}{\Vert }_2^2+\text{const.},$$

$$p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta })=\frac{p({\mathcal{D}}_s\mid \mathbf{y})p(\mathbf{y}\mid \mathbf{\theta })}{p({\mathcal{D}}_s\mid \mathbf{\theta })},$$

$$p({\mathcal{D}}_s\mid \mathbf{\theta })=\int p({\mathcal{D}}_s\mid \mathbf{y})p(\mathbf{y}\mid \mathbf{\theta })\mathrm{d}\mathbf{y}.$$

$$\hat{\mathbf{y}}\equiv \underset{\mathbf{y}}{\mathrm{argmax}}\log p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta })=\underset{\mathbf{y}}{\mathrm{argmax}}\log p(\mathbf{y},{\mathcal{D}}_s\mid \mathbf{\theta }),$$

$$p(\mathbf{y}\mid \mathbf{\theta })=\mathcal{N}(\mathbf{y}\mid 0,{\mathbf{C}}_p(\mathbf{\theta })\equiv C(\mathbf{\Xi },\mathbf{\Xi }\mid \mathbf{\theta })),$$

$$\begin{array}{rl}{D}_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta }))& =-\int q(\mathbf{y})\log \frac{p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta })}{q(\mathbf{y})}\mathrm{d}\mathbf{y}\\ & =-\int q(\mathbf{y})\log \frac{p({\mathcal{D}}_s\mid \mathbf{y})p(\mathbf{y}\mid \mathbf{\theta })}{q(\mathbf{y})p({\mathcal{D}}_s\mid \mathbf{\theta })}\mathrm{d}\mathbf{y}\\ & =-\mathcal{F}\left[q(\mathbf{y}),\mathbf{\theta }\right]+\log p({\mathcal{D}}_s\mid \mathbf{\theta }).\end{array}$$

$$\mathcal{F}\left[q(\mathbf{y}),\mathbf{\theta }\right]={\mathbb{E}}_{q(\mathbf{y})}\left[p({\mathcal{D}}_s\mid \mathbf{y})\right]-D_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid \mathbf{\theta })),$$

$$\log p({\mathcal{D}}_s\mid \mathbf{\theta })=\mathcal{F}\left[q(\mathbf{y}),\mathbf{\theta }\right]+D_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta })).$$

$$(\hat{q}(\mathbf{y}),\hat{\mathbf{\theta }})\equiv \underset{\left(q(\mathbf{y}),\mathbf{\theta }\right)}{\mathrm{argmax}}\mathcal{F}\left[q(\mathbf{y}),\mathbf{\theta }\right].$$

$${\hat{q}}^{(j+1)}(\mathbf{y})\text{ set as }\underset{q(\mathbf{y})}{\mathrm{argmax}}\mathcal{F}\left[q(\mathbf{y}),{\mathbf{\theta }}^{(j)}\right]$$

$${\hat{\mathbf{\theta }}}^{(j+1)}\text{ set as }\underset{\theta }{\mathrm{argmax}}\mathcal{F}\left[q^{(j+1)}(\mathbf{y}),\mathbf{\theta }\right].$$

$$\begin{gathered} \log p(\mathbf{y}\mid {\mathcal{D}}_s,\mathbf{\theta })=-\log p({\mathcal{D}}_s\mid \mathbf{\theta })+\log p(\hat{\mathbf{y}},{\mathcal{D}}_s\mid \mathbf{\theta }) \\ +\frac{1}{2}(\mathbf{y}-\hat{\mathbf{y}}{)}^{\top }\left[{\nabla \nabla \log p(\mathbf{y},{\mathcal{D}}_s\mid \mathbf{\theta })|}_{\mathbf{y}=\hat{\mathbf{y}}}\right](\mathbf{y}-\hat{\mathbf{y}})+\dots , \end{gathered}$$

$$\begin{array}{rl}{\mathbf{\mu }}_q& \equiv \underset{\mathbf{y}}{\mathrm{argmin}}\left[-\log p(\mathbf{y},{\mathcal{D}}_s\mid \mathbf{\theta })\right]\\ & =\underset{\mathbf{y}}{\mathrm{argmin}}\left\{-\log p({\mathcal{D}}_s\mid \mathbf{y})+\frac{1}{2}\left[{\mathbf{y}}^{\top }{\mathbf{C}}_p^{-1}(\mathbf{\theta })\mathbf{y}+\log \det {\mathbf{C}}_p(\mathbf{\theta })+N\log 2\pi \right]\right\},\end{array}$$

$${\mathbf{\Sigma }}_q\equiv -{\nabla \nabla \log p(\mathbf{y},{\mathcal{D}}_s\mid \mathbf{\theta })|}_{\mathbf{y}=\hat{\mathbf{y}}}=\mathbf{H}+{\mathbf{C}}_p^{-1}(\mathbf{\theta }),$$

$${\mathbf{\theta }}^{(j+1)}=\underset{\mathbf{\theta }}{\mathrm{argmin}}{D}_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid {\mathbf{\theta }}^{(j)})).$$

$$D_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid \mathbf{\theta }))=\frac{1}{2}\left[\mathrm{tr}\left({\mathbf{C}}_p^{-1}{\mathbf{\Sigma }}_q\right)+{\hat{\mathbf{y}}}^{\top }{\mathbf{C}}_p^{-1}\hat{\mathbf{y}}-N+\log \frac{\det {\mathbf{C}}_p}{\det {\mathbf{\Sigma }}_q}\right].$$

$$\frac{\partial }{\partial {\theta }_i}{D}_{\mathrm{KL}}(q(\mathbf{y})\Vert p(\mathbf{y}\mid \mathbf{\theta }))=-\frac{1}{2}{\hat{\mathbf{y}}}^{\top }{\mathbf{C}}_p^{-1}\frac{\partial {\mathbf{C}}_p}{\partial {\theta }_i}{\mathbf{C}}_p^{-1}\hat{\mathbf{y}}+\frac{1}{2}\mathrm{tr}\left[{\mathbf{C}}_p^{-1}\frac{\partial {\mathbf{C}}_p}{\partial {\theta }_i}\left({\mathbf{I}}_N-{\mathbf{C}}_p^{-1}{\mathbf{\Sigma }}_q\right)\right].$$

$$\max {\left\{\left|{\theta }_i^{(j+1)}-{\theta }_i^{(j)}\right|/\left|{\theta }_i^{\mathrm{s}}\right|\right\}}_{i=1}^{N_{\theta }}\le \mathrm{rtol},$$

$$(\hat{\mathbf{\varphi }},\hat{\mathbf{\theta }})=\underset{(\mathbf{\varphi },\mathbf{\theta })}{\mathrm{argmax}}\mathcal{F}\left[q(\mathbf{y}\mid \mathbf{\varphi }),\mathbf{\theta })\right].$$

$$\begin{gathered} \mathcal{F}\left[\mathbf{\varphi }(\mathbf{y}\mid \mathbf{\varphi }),\mathbf{\theta }\right]={\mathbb{E}}_{\mathcal{N}(\mathbf{z}\mid 0,{\mathbf{I}}_N)}[\log p({\mathcal{D}}_s\mid \mathbf{y})] \\ +{\mathbb{E}}_{\mathcal{N}(\mathbf{z}\mid 0,{\mathbf{I}}_N)}[\log p(\mathbf{y}\mid \mathbf{\theta })]+\log \det {\mathbf{R}}_q+\mathcal{H}[\mathcal{N}(\mathbf{z}|0,{\mathbf{I}}_N)], \end{gathered}$$

$$f(\mathbf{z};\mathbf{\varphi },\mathbf{\theta })=\log p({\mathcal{D}}_s\mid \mathbf{y})+\log \det {\mathbf{R}}_q-\frac{1}{2}\left[{\mathbf{y}}^{\top }{\mathbf{C}}_p^{-1}\mathbf{y}+\log \det {\mathbf{C}}_p-N\right],$$

$${\nabla }_{{\mathbf{\mu }}_q}f(\mathbf{z};\mathbf{\varphi },\mathbf{\theta })$$

$${\nabla }_{{\mathbf{R}}_q}f(\mathbf{z};\mathbf{\varphi },\mathbf{\theta })$$

$$\frac{\partial }{\partial {\theta }_i}f(\mathbf{z};\mathbf{\varphi },\mathbf{\theta })$$

$$f_n(\mathbf{\varphi },\mathbf{\theta })=\frac{1}{n}\sum _{k=1}^{n}f\left({\mathbf{z}}^{(k)};\mathbf{\varphi },\mathbf{\theta }\right),{\mathbf{z}}^{(k)}\sim \mathcal{N}(0,{\mathbf{I}}_N),$$

$${\nabla }_{(\cdot )}{f}_n(\mathbf{\varphi },\mathbf{\theta })=\frac{1}{n}\sum _{k=1}^n{\nabla }_{(\cdot )}f({\mathbf{z}}^{(k)};\mathbf{\varphi },\mathbf{\theta }),$$

$${\nabla }_{{\mathbf{\omega }}_q}f(\mathbf{z};\mathbf{\varphi },\mathbf{\theta })=\left[-{\nabla }_{\mathbf{y}}\log p({\mathcal{D}}_s\mid \mathbf{y})+{\mathbf{C}}_p^{-1}y\right]\circ \mathbf{z}\circ \exp ({\mathbf{\omega }}_q)-{\mathbf{I}}_N,$$

$$\frac{\partial }{\partial x}\left[k(x)\frac{\partial }{\partial x}u(x)\right]=0,x\in [0,1],$$

$$u(0)=u_{\mathrm{L}},u(1)=u_{\mathrm{R}},$$

$$C(x,x^{\prime }\mid \mathbf{\theta })={\sigma }^2\exp \left[-{\left(x-x^{\prime }\right)}^2/2{\lambda }^2\right]+{\sigma }_n^2{1}_{x=x^{\prime }},$$

$$\frac{\partial }{\partial x}\left[k(u(x))\frac{\partial }{\partial x}u(x)\right]=0,x\in [0,1],$$

Full text

1 Introduction

Partial differential equation (PDE) models of many physical systems involve space-dependent parameters and constitutive relationships that are usually only partially observed. Model inversion aims to estimate these unknown functions of space and the system’s state from sparse measurements of the state, associated quantities of interest, and the unknown functions themselves. Bayesian inference provides a probabilistic framework for model inversion [1], in which data is assimilated by computing the posterior density of the parameters in terms of the likelihood of the observations given the PDE model and the prior density of the parameters codifying modeling assumptions. Unlike deterministic parameter estimation methods [2, 3], the Bayesian framework provides a probabilistic characterization of the estimated parameter that can be employed for quantifying uncertainty and evaluating modeling assumptions. For linear problems with Gaussian likelihoods and Gaussian priors, Bayesian inference can be done exactly (known in the context of state estimation for dynamical systems as the Kalman filter [4]). Unfortunately, physics models define nonlinear maps between the state and the model parameters, preventing carrying out exact inference even in the case of Gaussian likelihoods and Gaussian priors. The Markov Chain Monte Carlo (MCMC) method is robust for general nonlinear problems but is computationally expensive [5]. Despite recent advances in Hamiltonian Monte Carlo and ensemble and parallel MCMC [6, 7, 8], the number of forward simulations and likelihood evaluations required by MCMC sampling poses a challenge for model inversion of PDE models with high-dimensional parameters. Here, we propose two cost-effective alternatives to MCMC for estimating unknown parameters and constitutive relationships in PDE models.

Gaussian process (GP) regression, known as kriging in spatial geophysics, is commonly used to construct probabilistic models of heterogeneous parameters; therefore, GPs serve as a reasonable choice of prior for unknown parameters. In the context of Bayesian inference with GP priors, GP regression is equivalent to exact Bayesian inference for assimilating direct measurements of unknown parameters. Similarly, the marginal likelihood of parameter measurements can be computed in closed form, therefore allowing for model selection to be carried out by empirical Bayesian inference, also known as type-II maximum likelihood estimation [9].

Assimilating measurements of the state of PDE models is on the other hand less straightforward. Recently, a framework has been proposed to combine GP priors on the state and a discretization of the governing PDEs to assimilate state observations [10, 11]. In this framework, state estimation and type-II maximum likelihood can be carried out in closed form when the governing equations are linear on the state; for the nonlinear case, inference is carried out approximately via linearization of the governing equations.

Parameter estimation for PDE models presents another layer of challenge as governing equations commonly induce nonlinear relations between parameters and states. A common example is the Laplace equation with space-dependent unknown diffusion coefficient, which is linear on the state, but induces a nonlinear relation between the state and the diffusion coefficient. For the general case of parameter estimation with nonlinearity introduced by the physics model, approximate Bayesian inference methods are necessary. The standard approximate inference tool is MCMC sampling of the Bayesian posterior. Given unbounded computational resources, MCMC will provide arbitrarily accurate results, but in practice MCMC often requires an intractable amount of forward simulations of the physics model. Algorithms such as Hamiltonian Monte Carlo (HMC) and the Metropolis-adjusted Langevin algorithm (MALA) employ first-order information in the form sensitivities of the physics model to improve the mixing and convergence of the Markov chain random walk, but nevertheless the total number of forward and sensitivity simulations remains a challenge. As an alternative, approaches such as the Laplace approximation and variational inference aim to approximate the exact posterior with an analytical, parameterized density.

In this manuscript we propose employing approximate Bayesian inference with GP priors to approximate the posterior of PDE parameters and to estimate the hyperparameters of their GP prior. We propose two optimization-based methods: The first, Laplace-EM, is based on the Laplace approximation [9, 12, 13] and the expectation maximization (EM) algorithm [14, 12]. The second, doubly stochastic variational inference for empirical Bayes inference (DSVI-EB) is based on the DSVI algorithm [15, 16]. The proposed methods employ first and second-order information, i.e., gradient and Hessian of physics models, evaluated via the discrete adjoint method. Both presented methods enjoy advantageous computational properties over MCMC and other approximate Bayesian inference algorithms such as expectation propagation [17] and the Laplace approximation-based method of [13] that renders each of them attractive for model inversion depending on the nature of the inversion problem. In particular, the Laplace-EM method is accurate for approximating the unimodal posteriors of the numerical examples of this manuscript, but requires computing Hessians. On the other hand, DSVI-EB is less accurate but only requires computing gradients, and can be trivially parallelized. We note that Gaussian mixtures can be employed in variational inference to approximate multimodal posteriors [18], but in the present work we limit our focus to unimodal posteriors. Furthermore, both methods are applicable to non-factorizing likelihoods, do not require computing moments of the likelihood, and do not require third- or higher order derivatives of the physics model. Finally, variational inference and the Laplace approximation have been employed for model inversion [19, 20, 21, 22, 13], but to the best of our knowledge, have not been used in the context of the empirical Bayesian framework to estimate GP prior hyperparameters, with the exception of the work of [13]. That work, based on the Laplace approximation, requires computing third-order derivatives of the physics model, which may be costly to compute, whereas the presented methods do not require third-order derivatives.

The manuscript is structured as follows: In Section 2 we formulate the empirical Bayesian inference problem for physics models and GP priors. In Section 3 we propose the approximate Bayesian inference and summarize the expectation maximization (EM) algorithm. The Laplace-EM algorithm is introduced in Section 4, and the DSVI-EB algorithm is described in Section 5. The computational complexity of the algorithms is discussed in Section 6. The application of the proposed methods are presented in Section 7. Finally, conclusions are given in Section 8.

2 Problem formulation

We consider physical systems modeled by a stationary PDEs over the simulation domain $\Omega\subset\mathbb{R}^{d}$, $d\in[1,3]$. We denote by $u\colon\Omega\to U\subset\mathbb{R}$ the system’s state, and by $y\colon\Omega\times U\to\mathbb{R}$ the system’s parameter, an unknown scalar function of space and the system’s state. Our goal is to estimate the unknown function $y(x,u)$ from sparse, noisy measurements of $u(\mathbf{x})$ and $y(\mathbf{x},u)$. The PDE and boundary conditions are discretized for numerical computations, resulting in the set of $M$ algebraic equations $\mathbf{L}(\mathbf{u},\mathbf{y})=0$, where $\mathbf{u}\in\mathbb{R}^{M}$ denotes the vector of $M$ state degrees of freedom, and $\mathbf{y}\in\mathbb{R}^{N}$ denotes the discretized parameter vector, corresponding to the value of $y(\mathbf{x},u)$ at the $N$ discrete locations $\{\xi_{i}\in\Omega\times U\}^{N}_{i=1}$. Furthermore, we denote by $\bm{\Xi}$ the matrix of coordinates $\bm{\Xi}\equiv(\xi_{1},\dots,\xi_{N})$.

We assume that the sparse observations of the discrete state and parameters, $\mathbf{u}_{\mathrm{s}}$ and $\mathbf{y}_{\mathrm{s}}$, respectively, are collected with iid normal observation errors, that is,

[TABLE]

where $\mathbf{u}_{\mathrm{s}}\in\mathbb{R}^{M_{s}}$, $M_{s}\ll M$ are the state observations, $\mathbf{y}_{\mathrm{s}}\in\mathbb{R}^{N_{s}}$, $N_{s}\ll N$ are the parameter observations, $\mathbf{H}_{u}\in\mathbb{R}^{M_{s}\times M}$ is the state observation operator, $\mathbf{H}_{y}\in\mathbb{R}^{N_{s}\times N}$ is the parameter observation operator, and $\bm{\epsilon}_{us}$ and $\bm{\epsilon}_{ys}$ are observation errors satisfying $\mathbb{E}[\bm{\epsilon}_{us}\bm{\epsilon}^{\top}_{ys}]=0$. Then, the likelihood of the observations $\mathcal{D}_{s}\equiv\{\mathbf{u}_{s},\mathbf{y}_{s}\}$ given $\mathbf{y}$ is defined as

[TABLE]

where $\mathbf{u}$ satisfies the physics constraint $\mathbf{L}(\mathbf{u},\mathbf{y})=0$ given $\mathbf{y}$, and the constant is independent of $\mathbf{y}$.

In probabilistic terms, our goal is to estimate the posterior density of $\mathbf{y}$ given the data $\mathcal{D}_{s}$. By Bayes’ theorem, this posterior is given by

[TABLE]

where $p(\mathbf{y}\mid\bm{\theta})$ is the parameterized prior density of $\mathbf{y}$, with hyperparameters $\bm{\theta}$, and $p(\mathcal{D}_{s}\mid\bm{\theta})$ is the marginal likelihood or evidence of the data, given by

[TABLE]

If one is not interested in the uncertainty in estimating $\mathbf{y}$ given the data, one can compute in lieu of the full posterior ((4)) the maximum a posteriori (MAP) point estimate of $\mathbf{y}$, defined as the mode of the posterior, that is,

[TABLE]

where $p(\mathbf{y},\mathcal{D}_{s}\mid\bm{\theta})=p(\mathcal{D}_{s}\mid\mathbf{y})p(\mathbf{y}\mid\bm{\theta})$ is the joint density of the data and the parameters given $\bm{\theta}$. Here we used the fact that the marginal likelihood $p(\mathcal{D}_{s}\mid\bm{\theta})$ is independent of $\mathbf{y}$.

We employ a zero-mean Gaussian process prior, that is,

[TABLE]

where $\mathcal{N}(\cdot\mid\bm{\mu},\bm{\Sigma})$ denotes the multivariate normal density with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$, and $C(\cdot,\cdot\mid\bm{\theta})$ is a parameterized covariance kernel.

The posterior density depends on the prior hyperparameters, which can be chosen based on prior expert knowledge, or estimated from data. In the empirical Bayes approach, also known as type-II maximum likelihood or marginal likelihood estimation, point estimates of the hyperparameters are obtained by maximizing the marginal likelihood with respect to $\bm{\theta}$, i.e., $\hat{\bm{\theta}}\equiv\operatorname*{arg\,max}_{\bm{\theta}}p(\mathcal{D}_{s}\mid\bm{\theta})$. In the fully Bayes approach, we instead pose a hyperprior on the hyperparameters, which is updated with data by the Bayes’ theorem. In this work we will pursue the empirical Bayes approach.

Due to the nonlinear map from $\mathbf{y}$ to $\mathbf{u}$ defined by the physics constraint, the Bayesian inference problem of evaluating the posterior and marginal likelihood cannot be done in closed form. Exact inference therefore requires sampling the posterior via MCMC, which is intractable for sufficiently large $N$ and $M$. As a consequence, estimating hyperparameters via marginal likelihood estimation is also intractable. As an alternative to exact inference, in this work we propose various approximate inference algorithms.

3 Approximate inference and Expectation Maximization

The goal is to approximate the exact posterior $p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta})$ by a density $q(\mathbf{y})$. The Kullback-Leibler (KL) divergence $D_{\mathrm{KL}}(q(\mathbf{y})\;\|\;p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta}))$ provides a means to rewriting the marginal likelihood, (5), in terms of $q(\mathbf{y})$. Namely, substituting (4) into the definition of the KL divergence gives

[TABLE]

where $\mathcal{F}\left[q(\mathbf{y}),\bm{\theta}\right]$ is given by

[TABLE]

and $\mathbb{E}_{q(\mathbf{y})}\left[\cdot\right]\equiv\int(\cdot)q(\mathbf{y})\,\mathrm{d}\mathbf{y}$ denotes expectation with respect to the density $q(\mathbf{y})$. Reorganizing (8) we have the following alternative expression for (5):

[TABLE]

Given that the KL divergence is always non-negative, we have the inequality $\log p(\mathcal{D}_{s}\mid\bm{\theta})\geq\mathcal{F}\left[q(\mathbf{y}),\bm{\theta}\right]$; therefore, the operator $\mathcal{F}$ is often called the evidence lower bound (ELBO). The inequality becomes an equality when $q(\mathbf{y})=p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta})$, that is, when the variational density is equal to the exact posterior. In the empirical Bayes setting, this suggest the strategy of selecting both $q$ and $\bm{\theta}$ by maximizing the ELBO [14], i.e.,

[TABLE]

Instead of maximizing over $q(\mathbf{y})$ and $\bm{\theta}$ simultaneously, we can do it iteratively by alternating two maximization steps, resulting in the iterative scheme

[TABLE]

thus recovering the expectation maximization (EM) algorithm [14].

It remains to specify how the maximization problems (10) and (11) will be solved, particularly with respect to how to optimize over the space of possible densities $q(\mathbf{y})$ approximating the true posterior. The approximate inference algorithms presented in this manuscript are based on two families of approximations of the posterior. The Laplace-EM algorithm (Section 4) uses a local approximation around the MAP for a given $\bm{\theta}$, and optimizes the ELBO using the EM algorithm, (11). The DSVI-EB algorithm (Section 5) uses a parameterized density $q(\mathbf{y}\mid\bm{\phi})$ with variational parameters $\bm{\phi}$ to be selected jointly with $\bm{\theta}$ via stochastic optimization.

4 Laplace-EM algorithm

The Laplace approximation is an approach for approximating unimodal posteriors densities. It consists of fitting a multivariate Gaussian density around the MAP for a given choice of hyperparameters $\bm{\theta}$. The $j$th E-step of the EM algorithm, (11), consists of finding the posterior for a given set of hyperparameters, $\bm{\theta}^{(j)}$. This suggests we can replace the E-step by a Laplace approximation to the posterior, giving raise to the Laplace-EM algorithm.

We proceed to briefly describe the Laplace approximation. Expanding up to second order the log posterior (see (4)) around the MAP ((6)) yields

[TABLE]

where $\left.\nabla\nabla\log p(\mathbf{y},\mathcal{D}_{s}\mid\bm{\theta})\right|_{\mathbf{y}=\hat{\mathbf{y}}}$ denotes the Hessian of the log joint density $\log p(\mathbf{y},\mathcal{D}_{s}\mid\bm{\theta})$ around the MAP. This quadratic expression suggests approximating the posterior by the multivariate Gaussian density $q(\mathbf{y})\equiv\mathcal{N}(\mathbf{y}\mid\bm{\mu}_{q},\bm{\Sigma}_{q})$ with mean $\bm{\mu}_{q}$ given by the MAP and covariance $\bm{\Sigma}_{q}$ given by the Hessian of the log joint density. In other words, we have

[TABLE]

and

[TABLE]

where $\mathbf{H}\equiv-\left.\nabla\nabla\log p(\mathcal{D}_{s}\mid\mathbf{y})\right|_{\mathbf{y}=\bm{\mu}_{q}}$ denotes the Hessian of the likelihood around $\bm{\mu}_{q}$. Note that $\bm{\Sigma}_{q}$ and $\bm{\mu}_{q}$ depend on $\bm{\theta}$ indirectly through the dependence of the MAP on $\bm{\theta}$.

In this work we solve the minimization problem (12) via gradient-based optimization. The necessary gradient of $\log p(\mathcal{D}_{s}\mid\mathbf{y})$ with respect to $\mathbf{y}$ is computed via the discrete adjoint method described in C. The Hessian of $\log p(\mathcal{D}_{s}\mid\mathbf{y})$ with respect to $\mathbf{y}$, necessary to evaluate (13), is also computed via the discrete adjoint method.

We propose the Laplace-EM algorithm, where the Laplace approximation provides an approximation to the E-step. For the M-step, we keep the Laplace approximation fixed and maximize the ELBO with respect to the hyperparameters of the GP prior, $\bm{\theta}$. From (9) we see that for fixed $q(\mathbf{y})$, $\bm{\theta}$ appears only through the KL divergence $D_{\mathrm{KL}}(q(\mathbf{y})\;\|\;p(\mathbf{y}\mid\bm{\theta}))$; therefore, it suffices to minimize this KL divergence at the M-step. The Laplace-EM M-step reads

[TABLE]

For the GP prior, this KL divergence is given in closed form by

[TABLE]

If using a gradient method, the gradient of the KL divergence is given by

[TABLE]

The Laplace-EM algorithm is summarized in Algorithm 1. In practice, the EM iterations are halted once either a maximum number of iterations are completed, or once the relative change in hyperparameters is below a certain threshold, that is, when

[TABLE]

where $N_{\theta}$ is the number of prior hyperparameters, the $\theta^{s}_{i}$, $i\in[1,N_{\theta}]$ are prescribed hyperparameter scales (that provide a sense of the magnitude of the hyperparameters), and $\mathrm{rtol}$ is the prescribed tolerance.

We note that the Laplace approximation is a commonly used tool for unimodal non-Gaussian inference [9, 12, 13]. Directly maximizing with respect to $\bm{\theta}$, the Laplace approximation to the marginal likelihood requires evaluating third-order derivatives of the log-likelihood function (3) with respect to $\mathbf{y}$. This is due to the implicit dependence of the MAP on $\bm{\theta}$, which requires evaluating third-order derivatives of the physics constraint. The use of the expectation maximization algorithm allows us to side-step the need of third-order derivatives. Other methods for non-Gaussian inference such as expectation propagation [17] require multiple evaluations of the moments of the likelihood function, and are therefore not considered in this work.

5 Doubly stochastic variational inference

In variational inference (VI) [23], we restrict our choice of $q(\mathbf{y})$ to a parameterized family $q(\mathbf{y}\mid\bm{\phi})$. In this context we refer to $q$ as the variational density and $\bm{\phi}$ as the variational parameters. Following (10), we will estimate the variational parameters and the GP prior hyperparameters simultaneously by maximizing the corresponding ELBO, that is,

[TABLE]

In this section we present our proposed implementation of variational inference for empirical Bayes. The main challenges of VI are (i) approximating the expectations on the expression for the ELBO, (9), and (ii) optimizing such approximations. To address these challenges we employ the doubly stochastic variational inference (DSVI) framework [15, 16], in which a noisy simple Monte Carlo estimate of the ELBO (1st source of stochasticity) is minimized via a gradient-based stochastic optimization algorithm (2nd source of stochasticity). In particular, we employ stochastic gradient ascent with the adaptive step-size sequence proposed by [24]. The gradients of the ELBO estimate with respect to variational parameters and prior hyperparameters are computed via Gaussian backpropagation [25, 16, 26, 24].

5.1 Gaussian backpropagation

To maximize the ELBO via gradient-based stochastic optimization, we construct unbiased estimates of the ELBO and its gradients with respect to $\bm{\phi}$ and $\bm{\theta}$. Computing the gradient $\nabla_{\bm{\phi}}\mathcal{F}$ is not trivial as it involves expectations over $q$, which depends on $\bm{\phi}$.

We restrict ourselves to the multivariate Gaussian variational family $q(\mathbf{y}\mid\bm{\phi})=\mathcal{N}(\mathbf{y}\mid\bm{\mu}_{q},\bm{\Sigma}_{q})$, with variational mean $\bm{\mu}_{q}\in\mathbb{R}^{N}$ and covariance $\bm{\Sigma}_{q}=\mathbf{R}_{q}\mathbf{R}_{q}^{\top}\in\mathbb{R}^{N\times N}$, where $\mathbf{R}_{q}$ is a lower triangular factor matrix. For this choice we have $\bm{\phi}=\{\bm{\mu}_{q},\mathbf{R}_{q}\}$. Similar to the Laplace approximation, this choice is justified for unimodal posteriors. We then introduce the change of variables $\mathbf{y}=\bm{\mu}_{q}+\mathbf{R}_{q}\mathbf{z}$, with $\mathbf{z}\sim\mathcal{N}(0,\mathbf{I}_{N})$. Substituting this change of variables into (9), we can rewrite the involved expectations in terms of expectations over $\mathcal{N}(\mathbf{z}\mid 0,\mathbf{I}_{N})$, resulting in

[TABLE]

where $\mathbf{y}=\bm{\mu}_{q}+\mathbf{R}_{q}\mathbf{z}$, and $\mathcal{H}[\mathcal{N}(\mathbf{z}\mid 0,\mathbf{I}_{N})]=N(1+\log 2\pi)/2$ is the differential entropy of the standard multivariate normal. We then define the following unbiased estimate of the ELBO,

[TABLE]

with gradients

[TABLE]

where again $\mathbf{y}=\bm{\mu}_{q}+\mathbf{R}_{q}\mathbf{z}$. The details of the derivations of (20)–(22) are presented in A. It can be verified that $\mathbb{E}_{\mathcal{N}(\mathbf{z}\mid 0,\mathbf{I}_{N})}[f(\mathbf{z};\bm{\phi},\bm{\theta})]=\mathcal{F}\left[q(\mathbf{y}\mid\bm{\phi}),\bm{\theta}\right]$, so that the estimates are unbiased.

The variance of the estimate (19) and its gradients can be reduced by using the simple Monte Carlo (MC) or batch estimate and the corresponding gradients

[TABLE]

where $n$ is the size of the batch. The variance of the batch estimate (23) is lower by a factor of $n$, but requires computing the gradients of the log-likelihood $n$ times.

Our DSVI algorithm is summarized in Algorithm 2. The stochastic gradient ascent algorithm with adaptive step-size sequence, proposed by [24], is reproduced in B for completeness.

5.2 Parameterization of the factor matrix $\mathbf{R}_{q}$

It remains to discuss the parameterization of the factor $\mathbf{R}_{q}$. In this manuscript, we consider three alternatives: a full parameterization, the so-called mean field parameterization, and a constrained Chevron parameterization. The sparsity patterns of these parameterizations are shown in Figure 1.

In the full rank parameterization [24], we take $\mathbf{R}_{q}$ to be the non-unique Cholesky factor, that is, a $N\times N$ lower triangular matrix with unconstrained entries (Figure 1(a)). In this case, we have $\bm{\phi}\in\mathbb{R}^{N+N(N+1)/2}$. The number of variational parameters is therefore $O(N^{2})$, which may render their optimization difficult. In order to address this challenge, we can employ the mean field and Chevron parameterizations, which result in a total number of variational parameters that is linear on $N$.

In the mean field parameterization, we take $\mathbf{R}_{q}$ to be a strictly positive diagonal matrix, i.e., $\mathbf{R}_{q}=\operatorname{diag}[\exp(\bm{\omega}_{q})]$, with $\bm{\omega}_{q}\in\mathbb{R}^{N}$ (Figure 1(b)) and $\exp(\cdot)$ understood as element-wise. This parameterization assumes that the variational density covariance is diagonal, and the exponential ensures that the non-zero entries of $\mathbf{R}_{q}$ are strictly positive. In this case, we have $\bm{\phi}\equiv\{\bm{\mu}_{q},\bm{\omega}_{q}\}\in\mathbb{R}^{2N}$. The gradient of $f(\mathbf{z};\bm{\phi},\bm{\theta})$ with respect to $\bm{\omega}_{q}$ is given by

[TABLE]

where $\circ$ denotes element-wise product, and exponentiation is taken as element-wise. The details of the derivation are presented in A. This parameterization assumes that the posterior components of $\mathbf{y}$ are essentially uncorrelated and therefore cannot resolve the correlations of the true posterior, which are expected to be non-trivial for highly correlated prior covariance structures and for a small number of observations. As a consequence, the mean field parameterization tends to underestimate the posterior variance [23].

Finally, the constrained Chevron parameterization is similar to the full parameterization, but we set entries below the diagonal and for column number larger than the Chevron parameter $k<N$ to zero (Figure 1(c)) [27]. In this case, we have $\bm{\phi}\in\mathbb{R}^{N+(2N-k)(k+1)/2}$. The number of variational parameters for this parameterization is $O(N(k+1))$, a reduction with respect to the full parameterization, while maintaining some degree of expressivity for capturing correlations of the true posterior.

The DSVI-EB method follows the automatic differentiation variational inference (ADVI) algorithm [24], in which gradients of the joint probability $p(\mathcal{D}_{s}\mid\mathbf{y})$ with respect to variational parameters are computed using Gaussian backpropagation and reverse-mode automatic differentiation. ADVI is formulated for the full Bayes case and implements the full and mean-field parameterizations of $\mathbf{R}_{q}$. In comparison, our work is formulated for the empirical Bayes case, employs the adjoint method to compute gradients of physics solvers, and implements the constrained Chevron parameterization in addition to the full and mean-field parameterizations.

Two schemes are common in the literature for computing the gradients of the noisy ELBO estimate: the reinforce algorithm [28], also known as the likelihood ratio method or the log-derivative trick, and Gaussian backpropagation [26], also known as the reparameterization trick [25, 16]. The reinforce algorithm employs gradients of the variational density with respect to its parameters, which is convenient as it only requires zero-order information of the physics model. Unfortunately the reinforce estimates of the ELBO gradients are well-known to be of high variance and must be paired with a variance reduction technique [15]. Gaussian backpropagation, on the other hand, results in lower-variance gradient estimates at the cost of requiring first-order information of the physics model.

An alternative formulation of VI is presented in [18], where the authors employ mixtures of diagonal multivariate Gaussian densities as the variational posterior, and approximate the ELBO using a second order Taylor expansion around the mean of each mixture component. The mean, diagonal covariance and mixture weights are estimated by maximizing the ELBO via coordinate ascent. This entirely deterministic approach is formulated for the inference problem and does not consider optimization over prior hyperparameters. A similar approach is also presented in [29] in the context of empirical Bayes.

6 Computational cost

In this section, we discuss the computational effort of the Laplace-EM and DSVI-EB algorithms. We compute the gradient and the Hessian of the likelihood via the discrete adjoint method (see C for details). Note that the Laplace-EM method requires both gradients and Hessians, while the DSVI-EB method only requires gradients. For the physics constraint $\mathbf{L}(\mathbf{u},\mathbf{y})=0$, the computation of the gradient requires the solution of one (linear) backward sensitivity problem of size $M\times M$. Similarly, the computation of the Hessian requires one backward sensitivity problem and $N$ forward sensitivity problems, each of size $M\times M$. For the following discussion we assume that the cost of each forward and backward sensitivity problem is of order $O(M^{\gamma})$, $\gamma>1$.

For the Laplace-EM algorithm we discuss the cost per each EM cycle. Each E-step requires one Cholesky factorization of $\mathbf{C}_{p}(\bm{\theta}^{(j)})$, of cost $N^{3}$, the solution of (12) via gradient-based optimization, and the computation of the Hessian. Therefore, the cost of each E-step is $O(\max(M^{\gamma},N^{3}))$. Each iteration of the M-step requires one Cholesky factorization of $\mathbf{C}_{p}(\bm{\theta})$. Therefore, the total cost of each EM cycle is again $O(\max(M^{\gamma},N^{3}))$.

For the DSVI-EB algorithm, each iteration requires evaluating $n$ gradients and one Cholesky factorization of $\mathbf{C}_{p}$. If $n$ is chosen independent of $M$, we have that the total cost per iteration is also $O(\max(M^{\gamma},N^{3}))$.

Finally, in general we expect the number of iterations for each $E$- and $M$-step, and the number of EM cycles and DSVI iterations, to increase with increasing $N$. The analysis of how said numbers scale with $N$ is beyond the scope of this manuscript.

7 Numerical experiments

In this section, we present the application of the proposed approximate inference algorithms to the identification of the diffusion coefficient in diffusion equations. In particular, we are interested in identifying the diffusion coefficient $k(\mathbf{x},u)$ of the homogeneous diffusion equation $\nabla\cdot(k(\mathbf{x},u)\nabla u)=0$ in $\Omega\subset\mathbb{R}^{d}$, from both measurements of the diffusion coefficient and of the state $u$. For the linear case ($k\equiv k(\mathbf{x})$), the diffusion equation models phenomena such as stationary heat transfer and Darcy flow. For the nonlinear case ($k\equiv k(u)$), one recovers the so-called Richards equation for horizontal flows in unsaturated porous media.

7.1 Linear diffusion problem

We consider the one-dimensional diffusion equation with Dirichlet boundary conditions

[TABLE]

where $u\colon[0,1]\to\mathbb{R}$ is the state and $k\colon[0,1]\to\mathbb{R}^{+}$ is the diffusion coefficient. The state is discretized into $M$ degrees of freedom $u_{i}$ organized into the vector $\mathbf{u}\in\mathbb{R}^{M}$. The diffusion coefficient is discretized into $N$ degrees of freedom $k_{i}=\exp y_{i}$ corresponding to $N$ spatial coordinates $\{x_{i}\}^{N}_{i=1}$, with the $y_{i}$ organized into the vector $\mathbf{y}\in\mathbb{R}^{N}$. The discretized problem (26) and (27) is of algebraic form $\mathbf{L}(\mathbf{u},\mathbf{y})\equiv\mathbf{S}(\mathbf{y})\mathbf{u}-\mathbf{b}(\mathbf{y})=0$, where $\mathbf{S}\colon\mathbb{R}^{N}\to\mathbb{R}^{M\times M}$ and $\mathbf{b}\colon\mathbb{R}^{N}\to\mathbb{R}^{M}$. In (26), $\mathbf{y}$ can only be identified from measurements of $u$ up to an additive constant [1], and measurements of $\mathbf{y}$ are required to estimate it uniquely.

We apply the presented model inversion algorithms to estimating a synthetic diffusion coefficient from $10$ measurements of the state $u$ and one measurement of the log-diffusion coefficient $y$. The reference values of $\mathbf{y}$ and $\mathbf{u}$ and the corresponding observations are shown in Figure 2. The reference $\mathbf{y}$ is taken as a realization of the zero-mean GP with squared exponential covariance

[TABLE]

where $\bm{\theta}\equiv(\sigma,\lambda)$, and $\sigma_{n}$ is set to $1\text{\times}{10}^{-2}$. We refer to $\sigma$ and $\lambda$ as the standard deviation and correlation length, respectively, of the prior covariance. The reference values of $\bm{\theta}$ are presented in Table 1. The $\mathbf{y}$ and $\mathbf{u}$ observations are taken at randomly selected degrees of freedom, with observation error standard deviations $\sigma_{us}=\sigma_{ys}=$1\text{\times}{10}^{-3}$$. Finally, the boundary conditions are set to $u_{\mathrm{L}}=1.0$ and $u_{\mathrm{R}}=0.0$, and the numbers $M$ and $N$ of degrees of freedom is set to $50$.

7.1.1 Empirical Bayesian inference

Figure 3 shows the estimated diffusion coefficient, together with the $95\%$ confidence intervals centered around the posterior mean, computed using the Laplace-EM algorithm and DSVI-EB with Chevron parameterization and $k=20$. It can be seen that both methods provide accurate estimates of the reference field, with the reference field falling inside the confidence interval of the estimates (with localized exceptions for DSVI-EB with Chevron parameterization in the vicinity of the $x=1.0$ boundary, as shown in Figure 3(b)).

The estimated prior hyperparameters are presented in Table 1, together with a simple MC estimate of the ELBO, $\hat{\mathcal{F}}$, computed using $1\text{\times}{10}^{4}$ realizations of the estimated posterior. It can be seen that hyperparameter estimates are very similar for the Laplace-EM method and the DSVI method. In particular, estimates of the correlation length are close to reference values, while the standard deviation is underestimated across all methods. Estimates are also similar for the different factor matrix parameterizations on the DSVI method, with the exception of the full rank parameterization that resulted in more pronounced underestimation of both the standard deviation and correlation length. In terms of the ELBO, the Laplace-EM method results in the highest value, followed by the DSVI method with full rank factor parameterization. This indicates that full rank representations of the covariance matrix of the estimated posterior density result in better estimates of the true posterior, whereas reduced representations such as Chevron and mean field are less accurate. In practice, it can be seen in Figure 3(b) that the reduced representation is less capable of resolving the uncertainty of the $y$ estimate in the vicinity of the $x=1.0$ boundary. Nevertheless, reduced representations are not worse than fuller representations for estimating prior hyperparameters.

7.1.2 Comparison against MCMC

We proceed to evaluate the accuracy of the proposed inference algorithms at approximating the posterior density $p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta})$. For this purpose, we employ MCMC simulation as the benchmark as it is known to converge to the exact posterior density. In order to restrict the focus to the approximation of the posterior density, we set the prior hyperparameters to fixed values equal to the reference values, $\bm{\theta}_{\mathrm{ref}}$, and employ the Laplace-EM 111Note that in this context the Laplace-EM algorithm is reduced to the Laplace approximation for given $\bm{\theta}$ (e.g. a single E-step in the EM algorithm), but we will refer to the associated results as Laplace-EM results for the sake of convenience. and DSVI algorithms to estimate the posterior $p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta}_{\mathrm{ref}})$.We compare the estimated posterior mean and standard deviation against the sample mean and standard deviation computed from $1\text{\times}{10}^{4}$ MCMC realizations of the posterior generated using the No-U-Turn Sampler (NUTS) [8].

Figure 4 presents a point-wise comparison against MCMC of the estimated mean and standard deviation computed using both Laplace-EM and DSVI with the full rank parameterization and the Chevron parameterization with $k=20$ and $5$. It can be seen that all estimates of the mean are very accurate, which indicates that the presented approximate inference algorithms provide accurate estimates of the mean of multimodal posterior densities. This result is expected for the Laplace-EM algorithm where the estimated posterior mean is set to the MAP, but for the DSVI algorithm this is less of a given.

For the standard deviation, the Laplace-EM method provides the most accurate estimates, followed by DSVI with the full rank parameterization. This reinforces the conclusion drawn previously that full rank representations of the covariance lead to better estimates of the true posterior. Furthermore, it can be seen that the Chevron representation accurately resolves the bulk of point-wise standard deviation values (clustered at the bottom left of each plot in Figure 4(b)) but leads to noticeable underestimation of the larger point-wise values (i.e. the top half of Figure 4(b)). The underestimation of the standard deviation is more pronounced for decreasing $k$, and is the most pronounced for the mean field parameterization (not shown), which as remarked previously tends to underestimate the variance of the posterior [23].

Finally, in Figure 5 we present the posterior mean and variance for Laplace-EM and DSVI with Chevron parameterization and $k=20$, obtained for fixed $\bm{\theta}=\bm{\theta}_{\mathrm{ref}}$. Comparing Figure 3 against Figure 5 reveals that even though the empirical Bayes estimation procedure results in a standard deviation estimate lower than the reference value (see Table 1), the posterior density with empirical Bayes hyperparameter estimates, $p(\mathbf{y}\mid\mathcal{D}_{s},\hat{\bm{\theta}})$, is a good approximate to the posterior density with reference hyperparameters, $p(\mathbf{y}\mid\mathcal{D}_{s},\bm{\theta}_{\mathrm{ref}})$.

7.2 Nonlinear diffusion problem

We consider the one-dimensional nonlinear diffusion equation with Dirichlet boundary conditions

[TABLE]

where $u\colon[0,1]\to(-\infty,0]$ is the state and $k\colon(\infty,0]\to\mathbb{R}^{+}$ is the diffusion coefficient. Similarly to Section 7.1, the state is discretized into $M$ degrees of freedom $u_{i}$ organized into the vector $\mathbf{u}\in\mathbb{R}^{M}$. The diffusion coefficient function is discretized into $N$ degrees of freedom $k_{i}=\exp y_{i}$ corresponding to $N$ values of $u$ over $[u_{\min},0]$ (where $u_{\mathrm{min}}<u_{\mathrm{L}}$), organized into the vector $\mathbf{y}\in\mathbb{R}^{N}$. The discretized problem (29) and (30) is of the algebraic form $\mathbf{L}(\mathbf{u},\mathbf{y})\equiv\mathbf{S}(\mathbf{u},\mathbf{y})\mathbf{u}-\mathbf{b}(\mathbf{u},\mathbf{y})=0$, where $\mathbf{S}\colon\mathbb{R}^{M}\times\mathbb{R}^{N}\to\mathbb{R}^{M\times M}$ and $\mathbf{b}\colon\mathbb{R}^{M}\times\mathbb{R}^{N}\to\mathbb{R}^{M}$. Inspection of (29) reveals that $y(u)\equiv\log k(u)$ can only be identified over the range $[u_{\mathrm{L}},u_{\mathrm{R}}]$ and up to an additive constant222This can be verified by introducing the Kirchhoff transformation $f(u)=\int^{u}_{u_{\min}}k(u)\,\mathrm{d}u$, with which (29) can be written as a linear equation on $f$..To disambiguate the estimate, we provide measurements of $y(u)$ at $u=u_{\min}$ and $u=0$.

We apply the presented model inversion algorithms to estimating a function $k(u)$ from $5$ measurements of the state $u$ and 2 measurements of $y(u)\equiv\log k(u)$ at $u=u_{\min}$ and $u=0$. The reference diffusion coefficient is $k(u)=\exp u$ ($y(u)=u$). The $\mathbf{u}$ observations are taken at randomly selected degrees of freedom, and are shown in Figure 6. Observation error standard deviations $\sigma_{us}$ and $\sigma_{ys}$ are set to $1\text{\times}{10}^{-2}$. Boundary conditions are set to $u_{\mathrm{L}}=-2.0$ and $u_{\mathrm{R}}=-0.5$, and $u_{\min}$ is set to $-2.5$. Finally, the numbers $M$ and $N$ are set to $50$ and $21$, respectively.

Figure 7 shows the estimated diffusion coefficient using the proposed model inversion methods, together with the 95% confidence intervals centered around the posterior mean. Presented are the results for the Laplace-EM method and the DSVI-EB method with Chevron parameterization and $k=5$. It can be seen that the estimated posterior mean and confidence intervals for both methods are nearly identical. As prior covariance $C(u,u\mid\bm{\theta})$, we employ the squared exponential model (28) with $\sigma_{n}$ set to $1\text{\times}{10}^{-2}$. As in the linear case, both methods accurately estimate the reference function $y(u)$, and the reference function falls inside the 95% confidence intervals provided by the estimated posterior covariance.

Table 2 presents the estimated hyperparameters of the prior for the Laplace-EM and the DSVI-EB method, together with simple MC estimates of the ELBO computed using $1\text{\times}{10}^{4}$ realizations of the corresponding estimated posterior densities. It can be seen that estimated hyperparameters are different for the different methods (note that here we don’t have reference values for the hyperparameters of the prior, as the reference $k(u)$ is not drawn from a GP model). Nevertheless, it can be seen that both methods result in similar estimates of $y$. In agreement with the linear case, the Laplace-EM and the DSVI-EB method with full rank parameterization result in the largest values of ELBO. Additionally, it can be seen that the ELBO decreases with increasing sparsity of the posterior covariance factor parameterization (i.e. with decreasing Chevron factor $k$), being the lowest for the mean field parameterization. This illustrates the compromise between the sparsity of the covariance factor and its expressive capacity for approximating the true posterior, that is, that less sparse covariance factors produce more accurate approximate posteriors.

8 Conclusions and discussion

We have presented two approximate empirical Bayesian methods, Laplace-EM and DSVI-EB, for estimating unknown parameters and constitutive relations in PDE models. Compared to other methods for approximate Bayesian inference, the proposed methods do not require third-order derivatives of the physics model, do not involve computing moments of non-Gaussian likelihoods, and are applicable to non-factorizing likelihoods. Furthermore, the calculation of the batch estimate of the ELBO and its gradients employed in the DSVI-EB method is trivially parallelizable, leading to savings in computational time. The numerical experiments presented show that both methods accurately approximate the posterior density and the hyperparameters of the GP prior. In particular, we find that the Laplace-EM method is more accurately approximate the posterior density, at the cost of computing Hessians of the physics model, which increase the computational cost of each EM cycle. The DSVI-EB method, on the other hand, is less accurate but does not require Hessians. Consistent with the literature, we find that the accuracy of the DSVI-EB method at approximating the posterior decreases with increasing sparsity of the covariance factor parameterization employed.

For a very large number of degrees of freedom of the discretization of the unknown functions, the computational cost of the proposed methods is dominated by the associated cubic complexity. Future work will aim to address the challenge of cubic complexity by employing sparse GP inference.

Appendix A Gaussian backpropagation rules

In this section we present the derivation of the gradients (20)–(22).

For the gradient with respect to the variational mean, (20), we have by the chain rule, in index notation,

[TABLE]

where we have used $\partial y_{j}/\partial\mu_{q,i}=\delta_{ji}$. It follows that $\nabla_{\bm{\mu}_{q}}\log p(\mathcal{D}_{s}\mid\mathbf{y})=\nabla_{\mathbf{y}}\log p(\mathcal{D}_{s}\mid\mathbf{y})$. Similarly, we have

[TABLE]

and thus we recover (20).

For the gradient with respect to the Cholesky factor $\mathbf{R}_{q}$, we note that $\partial y_{k}/\partial R_{q,ij}=\delta_{ki}z_{j}$. By the chain rule, we have

[TABLE]

so that $\nabla_{\mathbf{R}_{q}}\log p(\mathcal{D}_{s}\mid\mathbf{y})=\left[\nabla_{\mathbf{y}}\log p(\mathcal{D}_{s}\mid\mathbf{y})\right]\mathbf{z}^{\top}$. Similarly,

[TABLE]

Finally, we have $\nabla_{\mathbf{R}_{q}}\log\det\mathbf{R}_{q}=(\mathbf{R}^{-1}_{q})^{\top}$, from which we recover (21).

The gradients with respect to the prior hyperparameters, (22) can be derived from [9], Eqs. (A.14) and (A.15).

For the gradient with respect to $\bm{\omega}_{q}$ of the mean-field parameterization $\mathbf{R}_{q}=\operatorname{diag}[\exp{\bm{\omega}_{q}}]$, we employ the relation

[TABLE]

By the chain rule, we have

[TABLE]

summation over $k$ not implied. It follows that $\nabla_{\bm{\omega}_{q}}\log p(\mathcal{D}_{s}\mid\mathbf{y})=\left[\nabla_{\mathbf{y}}\log p(\mathcal{D}_{s}\mid\mathbf{y})\right]\circ\mathbf{z}\circ\exp{\bm{\omega}_{q}}$. Similarly,

[TABLE]

summation over $k$ not implied. Finally,

[TABLE]

from which we recover (25).

Appendix B Stochastic gradient ascent with adaptive step-size sequence

Here we reproduce for completeness the stochastic gradient ascent algorithm with adaptive step-size sequence proposed in [24]. The presentation is expanded to the empirical Bayes context for the update of prior hyperparameters. At each iteration, the variational parameters and prior hyperparameters are updated using the rules

[TABLE]

where $f^{(j)}_{n}\equiv f_{n}(\bm{\phi}^{(j)},\bm{\theta}^{(j)})$, and the vectors of step-sizes $\bm{\rho}^{(j)}_{\bm{\phi}}$ and $\bm{\rho}^{(j)}_{\bm{\theta}}$ are given by

[TABLE]

and the sequence

[TABLE]

where $\sqrt{\cdot}$ and $(\cdot)^{2}$ are understood as element-wise. The parameters $\tau$, $\alpha$, and $\epsilon$ are set to $1.0$, $0.1$, and $1\text{\times}{10}^{-16}$, respectively, while the parameter $\eta>0$ is chosen on a case-by-case basis.

Appendix C Discrete adjoint method for Darcy flow

In this section we describe the computation of the gradient and Hessian of the log-likelihood, $\nabla_{\mathbf{y}}\log p(\mathcal{D}_{s}\mid\mathbf{y})$ via the discrete adjoint method [30, 31]. For this purpose we introduce the function

[TABLE]

so that $\nabla_{\mathbf{y}}\log p(\mathcal{D}_{s}\mid\mathbf{y})=\nabla h(u(\mathbf{y}),\mathbf{y})$ by virtue of (3) (as the constant in (3) is independent of $\mathbf{y}$). In the following we will employ the following notation: Let $a$ be a scalar function, $\mathbf{b}$ and $\mathbf{c}$ be vector functions, and $\gamma$ be a scalar variable; then, $\partial a/\partial\mathbf{b}$ denotes the row vector with entries $\partial a/\partial b_{i}$, $\partial\mathbf{b}/\partial\gamma$ denotes the column vector with entries $\partial b_{i}/\partial\gamma$, and $\partial\mathbf{b}/\partial\mathbf{c}$ be the matrix with $ij$th entry $\partial b_{i}/\partial c_{j}$.

Differentiation $h(\mathbf{u},\mathbf{y})$ with respect to $y_{i}$ gives

[TABLE]

Similarly, differentiating the physics constraint $\mathbf{L}(\mathbf{u},\mathbf{y})=0$ with respect to $\mathbf{y}$ gives

[TABLE]

which implies $\partial\mathbf{u}/\partial y_{j}=-(\partial\mathbf{L}/\partial\mathbf{u})^{-1}(\partial\mathbf{L}/\partial y_{j})$. Substituting this relation into (36) gives the following expression for the $j$th component of the gradient:

[TABLE]

where the adjoint variables $\bm{\lambda}$ satisfies the adjoint equation

[TABLE]

It can be seen that computing the gradient $\nabla h(u(\mathbf{y}),\mathbf{y})$ requires a single linear backward sensitivity problem, (39), of size $M\times M$.

For the Hessian, we differentiate (36) with respect to $y_{i}$, obtaining

[TABLE]

where $D^{2}_{i,j}h$ is given by

[TABLE]

and $\partial^{2}h/\partial\mathbf{u}^{2}$ denotes the Hessian of $h$ with respect to $\mathbf{u}$, i.e. the matrix with $ij$th entry $\partial^{2}h/\partial u_{i}\partial u_{j}$. Similarly, differentiating (37) with respect to $y_{i}$ gives

[TABLE]

where $D^{2}_{i,j}\mathbf{L}$ is given element-wise in a manner similar to (41). (42) implies $\partial^{2}\mathbf{u}/\partial y_{i}\partial y_{j}=-(\partial\mathbf{L}/\partial\mathbf{u})^{-1}D^{2}_{i,j}\mathbf{L}$. Substituting into (40) gives the following expression for the $ij$th component of the Hessian:

[TABLE]

Computing the Hessian therefore requires the solution of $N$ linear forward sensitivity problems, (37), for each $\partial\mathbf{u}/\partial y_{i}$, and one backward sensitivity solution for the adjoint variables, each problem of size $M\times M$.

Acknowledgments

This work was supported by the Applied Mathematics Program within the U.S. Department of Energy Office of Advanced Scientific Computing Research. Pacific Northwest National Laboratory is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.