Gaussian Process Regression and Conditional Polynomial Chaos for Parameter Estimation

TL;DR

This paper introduces a novel Bayesian parameter estimation method combining Gaussian Process regression and conditional Polynomial Chaos to efficiently estimate unknown parameters in PDE models, reducing computational costs and improving accuracy.

Contribution

The paper presents a new data-driven surrogate modeling approach that integrates Gaussian Process regression with conditional Polynomial Chaos for PDE parameter estimation.

Findings

Reduces dimensionality of parameter space.

Replaces expensive PDE solutions with surrogate models.

Improves accuracy in estimating hydraulic conductivity.

Abstract

We present a new approach for constructing a data-driven surrogate model and using it for Bayesian parameter estimation in partial differential equation (PDE) models. We first use parameter observations and Gaussian Process regression to condition the Karhunen-Lo\'{e}ve (KL) expansion of the unknown space-dependent parameters and then build the conditional generalized Polynomial Chaos (gPC) surrogate model of the PDE states. Next, we estimate the unknown parameters by computing coefficients in the KL expansion minimizing the square difference between the gPC predictions and measurements of the states using the Markov Chain Monte Carlo method. Our approach addresses two major challenges in the Bayesian parameter estimation. First, it reduces dimensionality of the parameter space and replaces expensive direct solutions of PDEs with the conditional gPC surrogates. Second, the estimated parameter field exactly matches the parameter measurements. In addition, we show that the conditional gPC surrogate can be used to estimate the states variance, which, in turn, can be used to guide data acquisition. We demonstrate that our approach improves its accuracy with application to one- and two-dimensional Darcy equation with (unknown) space-dependent hydraulic conductivity. We also discuss the effect of hydraulic conductivity and head locations on the accuracy of the hydraulic conductivity estimations.

Equations

$$\begin{array}{rl}-\nabla \cdot (\hat{\kappa }(\mathbf{x})\nabla \hat{u}(\mathbf{x}))& =0,\mathbf{x}\in D;\\ \hat{u}(\mathbf{x})& =f(\mathbf{x}),\mathbf{x}\in \partial D_D;\\ \vec{n}\cdot \hat{\kappa }(\mathbf{x})\nabla \hat{u}(\mathbf{x})& =g(\mathbf{x}),\mathbf{x}\in \partial D_L.\end{array}$$

$$\begin{array}{rl}\nabla \cdot (\kappa (\mathbf{x},\omega )\nabla u(\mathbf{x},\omega ))& =0,\mathbf{x}\in D;\\ u(\mathbf{x},\omega )& =f(\mathbf{x}),\mathbf{x}\in \partial D_D;\\ \vec{n}\cdot \kappa (\mathbf{x},\omega )\nabla u(\mathbf{x},\omega )& =g(\mathbf{x}),\mathbf{x}\in \partial D_L.\end{array}$$

$${\omega }^{\ast }={\text{argmin}}_{\omega \in \Omega }\left[\sum _{j=1}^{N_k}|\hat{u}({\mathbf{x}}^j)-u({\mathbf{x}}^j,\omega ){|}^2+\sum _{j=1}^{N_m}|\hat{k}({\mathbf{x}}^j)-k({\mathbf{x}}^j,\omega ){|}^2\right].$$

$${\omega }^{\ast }={\text{argmin}}_{\omega \in \Omega }\sum _{j=1}^{N_k}|\hat{u}({\mathbf{x}}^j)-u({\mathbf{x}}^j,\omega ){|}^2,$$

$$Y(\mathbf{x},\omega )=\overline{Y}(\mathbf{x})+\sum _{n=1}^{\infty }\sqrt{{\lambda }_n}{\xi }_n(\omega ){ϵ}_n(\mathbf{x}),\text{in }{L}^2(\Omega ),$$

$$\lambda ϵ(\mathbf{x})={\int }_D{C}_Y(\mathbf{x},\mathbf{y})ϵ(\mathbf{y})d\mathbf{y}.$$

$$Y_N(\mathbf{x},\omega ):=\overline{Y}(\mathbf{x})+\sum _{n=1}^N\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\xi }_n(\omega ).$$

$${\Sigma }_{i,j}=C_Y({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})=\sum _{n=1}^{\infty }{\lambda }_n{ϵ}_n({\mathbf{x}}^{(i)}){ϵ}_n({\mathbf{x}}^{(j)}),i,j=1,\dots ,N_m,$$

$${ϵ}_n({\mathbf{x}}^{\ast })=({ϵ}_n({\mathbf{x}}^{(1)}),\dots ,{ϵ}_n({\mathbf{x}}^{(N_m)})).$$

$$C_{{\xi }_n,Y}({\mathbf{x}}^{(i)})=\sqrt{{\lambda }_n}{ϵ}_n({\mathbf{x}}^{(i)}),$$

$${\tilde{\mu }}_n$$

$${\tilde{m}}_{n,k}$$

$$\tilde{\mathbf{\eta }}=\tilde{\mathbf{\mu }}+\tilde{M}\mathbf{\eta },$$

$$\begin{array}{rl}{\tilde{Y}}_{N_G}(\mathbf{x},\omega )& =\overline{Y}(\mathbf{x})+\sum _{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\tilde{\eta }}_n\\ & =\overline{Y}(\mathbf{x})+\sum _{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\tilde{\mu }}_n+\sum _{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x})\sum _{l=1}^{N_G}{m}_{n,l}{\eta }_l.\end{array}$$

$$M{\Lambda }^{1/2}R={\Lambda }^{1/2}R-{\Lambda }^{1/2}R{\Sigma }^{-1}{R}^T{\Lambda }^{1/2}{\Lambda }^{1/2}R=0.$$

$$\begin{array}{rl}{C}_{\tilde{Y}}(\mathbf{x},\mathbf{y})& =\mathbb{E}[({\tilde{Y}}_{N_G}(\mathbf{x},\omega )-\mathbb{E}({\tilde{Y}}_{N_G}(\mathbf{x},\omega )))({\tilde{Y}}_{N_G}(\mathbf{y},\omega )-\mathbb{E}({\tilde{Y}}_{N_G}(\mathbf{y},\omega )))]\\ & =\mathbb{E}[\sum _{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x})\sum _{l=1}^{N_G}{\tilde{m}}_{n,l}{\eta }_l\sum _{s=1}^{N_G}\sqrt{{\lambda }_s}{ϵ}_s(\mathbf{y})\sum _{t=1}^{N_G}{\tilde{m}}_{s,t}{\eta }_t]\\ & =\sum _{n,s,l=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\tilde{m}}_{n,l}\sqrt{{\lambda }_s}{ϵ}_s(\mathbf{y}){\tilde{m}}_{s,l}\\ & =\mathbf{ϵ}{\Lambda }^{1/2}\tilde{M}{\tilde{M}}^T{\Lambda }^{1/2}{\mathbf{ϵ}}^T\\ & =\mathbf{ϵ}{\Lambda }^{1/2}\tilde{M}{\Lambda }^{1/2}{\mathbf{ϵ}}^T\end{array}$$

$${\tilde{Y}}_{N_G}(\mathbf{x},\omega )\sim \overline{Y}(\mathbf{x})+\sum _{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\tilde{\mu }}_n+\sum _{i=1}^{N_G-N_m}\sqrt{{\tilde{\lambda }}_i}{\widetilde{ϵ}}_i(\mathbf{x}){\xi }_i,$$

$$\begin{array}{rl}\tilde{\kappa }(\mathbf{x},\omega )& \sim \tilde{\kappa }(\mathbf{x},\mathbf{\xi })\\ & =e^{\overline{Y}(\mathbf{x})+{\sum }_{n=1}^{N_G}\sqrt{{\lambda }_n}{ϵ}_n(\mathbf{x}){\tilde{\mu }}_n}{e}^{{\sum }_{i=1}^{N_G-N_m}\sqrt{{\tilde{\lambda }}_i}{\widetilde{ϵ}}_i(\mathbf{x}){\xi }_i},\end{array}$$

$$\begin{array}{rl}-\nabla \cdot (\tilde{\kappa }(\mathbf{x},\mathbf{\xi })\nabla \tilde{u}(\mathbf{x},\mathbf{\xi }))& =0,\mathbf{x}\in D;\\ \tilde{u}(\mathbf{x},\mathbf{\xi })& =f(\mathbf{x}),\mathbf{x}\in \partial D_D;\\ \vec{n}\cdot \tilde{\kappa }(\mathbf{x},\mathbf{\xi })\nabla \tilde{u}(\mathbf{x},\mathbf{\xi })& =g(\mathbf{x}),\mathbf{x}\in \partial D_L.\end{array}$$

$$\tilde{u}(\mathbf{x},\mathbf{\xi })\approx \sum _{|\mathbf{i}|=0}^P{c}_{\mathbf{i}}(\mathbf{x}){\Phi }_{\mathbf{i}}(\mathbf{\xi }),$$

$$c_{\mathbf{i}}(\mathbf{x})=\mathbb{E}[\tilde{u}(\mathbf{x},\mathbf{\xi }){\Phi }_{\mathbf{i}}(\mathbf{\xi })]=\int \tilde{u}(\mathbf{x},\mathbf{\xi }){\Phi }_{\mathbf{i}}(\mathbf{\xi })\rho (\mathbf{\xi })d\mathbf{\xi }$$

$${\Phi }_{\mathbf{i}}(\mathbf{\xi })={\varphi }_{i_1}({\xi }_1)\dots {\varphi }_{i_N}({\xi }_N),0\le |\mathbf{i}|\le P,$$

$$\mathbb{E}[{\varphi }_m({\xi }_k){\varphi }_n({\xi }_k)]={\delta }_{m,n},0\le n,m\le N.$$

$$c_{\mathbf{i}}(\mathbf{x})\approx {\tilde{c}}_{\mathbf{i}}(\mathbf{x})=\sum _{m=1}^M\tilde{u}(\mathbf{x},{\mathbf{\xi }}^{(m)}){\Phi }_{\mathbf{i}}({\mathbf{\xi }}^{(m)})w^{(m)},$$

$$\tilde{u}(\mathbf{x},\mathbf{\xi })\approx {\tilde{u}}_C(\mathbf{x},\mathbf{\xi })=\sum _{|\mathbf{i}|=0}^P{\tilde{c}}_{\mathbf{i}}(\mathbf{x}){\Phi }_{\mathbf{i}}(\mathbf{\xi })$$

$$\mathbb{E}[\tilde{u}(\mathbf{x},\mathbf{\xi })]\approx \mathbb{E}[{\tilde{u}}_C(\mathbf{x},\mathbf{\xi })]={\tilde{c}}_0(\mathbf{x})$$

$$\mathbb{E}[\tilde{u}(\mathbf{x},\mathbf{\xi })-E[\tilde{u}(\mathbf{x},\mathbf{\xi })]{]}^2\approx \mathbb{E}[{\tilde{u}}_C(\mathbf{x},\mathbf{\xi })-E[{\tilde{u}}_C(\mathbf{x},\mathbf{\xi })]{]}^2=\sum _{|\mathbf{i}|=1}^P{\tilde{c}}_{\mathbf{i}}(\mathbf{x}{)}^2.$$

$${\mathbf{\xi }}^{\ast }={\text{argmin}}_{\mathbf{\xi }\in {\mathbb{R}}^{N_G-N_m}}\sum _{j=1}^{N_k}|\hat{u}({\mathbf{x}}^j)-{\tilde{u}}_C({\mathbf{x}}^j,\mathbf{\xi }){|}^2+\lambda \Vert \mathbf{\xi }-{\mathbf{\xi }}^o{\Vert }^2,$$

$$\hat{u}({\mathbf{x}}^j)={\tilde{u}}_C({\mathbf{x}}^j,\mathbf{\xi })+{\delta }_j,j=1,\dots ,N_k,$$

$$L(\mathbf{\xi })=\prod _{j=1}^{N_k}{p}_{\delta }\left({\delta }_j\right)=\prod _{j=1}^{N_k}{p}_{\delta }(\hat{u}({\mathbf{x}}^j)-{\tilde{u}}_C({\mathbf{x}}^j,\mathbf{\xi })).$$

$$p(\mathbf{\xi }|\hat{u}({\mathbf{x}}^j),j=1,\dots ,N_k)\propto p_{\mathbf{\xi }}\prod _{j=1}^{N_k}{p}_{\delta }(\hat{u}({\mathbf{x}}^j)-{\tilde{u}}_C({\mathbf{x}}^j,\mathbf{\xi })).$$

Full text

1 Introduction

Here, we focus on inverse problems, namely, parameter estimation in partial differential equations (PDEs) models. We are motivated by applications with high-dimensional parameter space. These problems are usually ill posed and deterministic parameter estimation methods are computationaly expensive and require regularization to obtain unique solutions [1, 3, 21, 12]. As an alternative, Bayesian methods have been proposed for obtaining the most probable combination of the parameters [21, 27]. Among Bayesian methods, the Markov chain Monte Carlo method provides a powerful tool for generating the posterior distributions of the parameters [9] and finding the most probable point in the parameter space for the “optimal” choice of parameters. However, MCMC requires a large number of forward solutions of the PDE model to generate posterior distributions that could be computationally expensive. Therefore, accurate and computational efficient surrogates for the PDE model are often sought [4].

Various surrogates have been adopted to accelerate the Bayesian inference for inverse problems [8, 23, 6, 2]. For example, the Gaussian process regression (GPR) or Kriging model was used in [37] to approximate the relation between parameters and the forward solutions. In [20, 19, 18], the generalized polynomial chaos (gPC) expansion was employed to approximate the forward solution of the PDE models as a function of the random parameters in the Karhunen-Loéve (KL) expansion of the unknown coefficient.

In this work, we propose a novel gPC-based surrogate model, which we call the conditional gPC surrogate model, and use it to approximate the solution of a PDE model in the MCMC approach for parameter estimation. The easy-to-evaluate conditional gPC surrogate significantly accelerates the MCMC sampling and the parameter estimation and quarantees that the estimated parameters exactly match the parameter measurements. The conditional gPC is based on the data-driven conditional KL representation of unknown space-dependent parameters [25, 14, 15] that reduces the dimensionality of the parameter space and allows parameter estimation with relatively few measurements. There are two levels of the dimension reduction in the conditional KL and gPC methods. First, the exact infinite KL representation of a random parameter is approximated with a finite KL expansion [16, 10, 26, 32]. At the second level, available observations of the parameter are used to condition the truncated KL expansion and further reduce the dimensionality. The reduction of dimensionality by conditioning on the parameter measurements partially alleviates the curse of dimensionality and allows a high-order gPC construction with relatively few samples (collocation points).

We use the conditional gPC surrogate model in combination with MCMC to estimate a partially observed space-dependent diffusion coefficient in the steady-state diffusion equation given partial observations of the state variables. We model the unknown parameters and states as random fields and represent the parameters with the conditional KL expansion. The partial measurements of the (random) state variable are treated as measurements of one realization of this variable and a minimization problem is formed to estimate the parameter by computing the optimal coefficients (values of the random variables) in the conditional KL expansion. Furthermore, we evaluate several sampling strategies of the states. We find that the gPC surrogate is more sensitive with respect to the parameters at the locations where it has higher variance. Therefore, measuring the state variable where its variance is largest allows more accurate parameter estimation as compared to uniformly and randomly distributed state measurements. Since in practice only few measurements of state are available, the “optimal” choice of measurement locations is important. Moreover, this strategy can be extended to the other regression approaches under the probabilistic framework.

The remainder of this paper is organized as follows. In Section  2, we formulate the spatial dependent coefficient estimation problem via the parameter minimization problem. The conditional KL representation with reduced dimensionality is stated in Section  3. In Section  4, we construct the conditional gPC surrogate and analyze its properties. In Section  5, we reformulate our minimization problem under the Bayesian framework and briefly introduce the MCMC method for parameter estimation. We study the accuracy of our approach through numerical examples, including parameter estimation in one- and two-dimensional elliptic PDEs in Section  6. Conclusions are given in Section  7.

2 Problem Formulation

Consider the (deterministic) steady-state diffusion equation with space-dependent partially known coefficient $\hat{\kappa}(\mathbf{x})$ and appropriate boundary conditions:

[TABLE]

Among other applications, Eq (1) describes flow in geological porous media, where $\hat{\kappa}(\mathbf{x})$ is the hydraulic conductivity and $\hat{u}(\mathbf{x})$ is the hydraulic head. In this application, for financial and technical reasons, $\hat{\kappa}(\mathbf{x})$ only can be measured in a few locations. Numerical treatments of Eq (1) are based on the idea that $\hat{\kappa}(\mathbf{x})$ can be resolved by a finite vector $\hat{\bm{\kappa}}:=\left(\hat{\kappa}(\mathbf{x}_{1}),\dots,\hat{\kappa}(\mathbf{x}_{n})\right)$ consisting of $\hat{\kappa}$ values at the collection of points $\{\mathbf{x}_{i}\}_{i=1}^{n}\subset D$.

We employ the probabilistic approach where we treat (the partially known) conductivity as a random field with prior distribution learned from the $\hat{k}$ measurements. Specifically, we employ the probability space $(\Omega,\mathcal{F},P)$ and assume that $\hat{\kappa}(\mathbf{x})$ is a realization of a spatially heterogenous random field, i.e., there exists a $P$-measurable map $\kappa(\cdot,\omega):\Omega\rightarrow L^{\infty}(D)$ and $\hat{\kappa}(\mathbf{x})=\kappa(\mathbf{x},\omega^{*})$.

This probabilistic approach renders Eq (1) stochastic:

[TABLE]

We propose to estimate $\hat{\kappa}(\mathbf{x})=\kappa(\mathbf{x},\omega^{*})$ given partial measurements $\{\hat{\kappa}(\mathbf{x}^{(i)})\}_{i=1}^{N_{m}}$ and $\{\hat{u}(\mathbf{x}^{(j)})\}_{j=1}^{N_{k}}$ by solving the following minimization problem,

[TABLE]

Solving this “full” parameter estimation problem with Bayesian inference is challenging specifically, it is challenging to minimize both terms in Eq (3). Instead, it is common to solve a simpler minimization problem:

[TABLE]

where $\hat{\kappa}$ measurements are only used to obtain a prior distribution of $k(\mathbf{x},\omega)$ [20, 19, 18]. Therefore, $\kappa(\mathbf{x},\omega^{*})$ found from the optimization problem (4) does not guarantee to much $\hat{\kappa}$ observations, while one found from (3) does. Following common practice, we assume that the conductivity $\kappa(\mathbf{x},\omega)$ in Eq (2) has a lognormal distribution, i.e., $Y(\mathbf{x},\omega)=\ln\kappa(\mathbf{x},\omega)$ has normal distribution [7]. In the work of [20, 19, 18], the “unconditional” KL expansion of $Y(\mathbf{x},\omega)$ and gPC representation of $u(\mathbf{x},\omega)$ in Eq (4) were used. Here, we use a “conditional” KL representation of $Y(\mathbf{x},\omega)$ and conditional gPC stochastic collocation method to construct a surrogate for $u(\mathbf{x},\omega)$ and solve the optimization problem (3). The conditional surrogate has two main advantages: it reduces the dimensionality of $\Omega$; and it exactly satisfies the second term in (3), i.e., $\hat{k}(\mathbf{x}^{j})\equiv\tilde{k}(\mathbf{x}^{j},\omega)$, where $\tilde{k}(\mathbf{x}^{j},\omega)$ is the KL expansion of $\kappa$ conditioned on $\{\hat{\kappa}(\mathbf{x}^{(i)})\}_{i=1}^{N_{m}}$. The latter reduces the optimization problem (3) to the simpler optimization problem (4). Both these advantages of the conditional gPC surrogate reduce the computational cost and improve accuracy of parameter estimation. The conditional KL and gPC methods are described in the following two sections. In Section 5, we present a method for solving the minimization problem (4) and demonstrate that its dimensionality is significantly smaller than those based on the unconditional KL representation of $Y(\mathbf{x},\omega)$.

3 Karhunen-Loève representation of $Y(\mathbf{x},\omega)$

KL representation of random coefficients combined with the gPC collocation method is the most common way to solve the stochastic PDE (2). The main challenge in this approach is that its computational cost exponentially increases with the number of terms in the truncated KL expansion. Until recently, gPC methods had been only applied to stochastic PDEs with random coefficients modeled as second-order stationary random fields, i.e., fields with constant variances and covariance functions depending only on the distance between two points. The number of terms in the truncated KL expansion of the stationary covariance function depends on the correlation length (it increases with the decreasing correlation length). Here, we condition the KL expansion of $Y(\mathbf{x},\omega)$ on the measurement $\{\ln\hat{\kappa}(\mathbf{x}^{(i)})\}_{i=1}^{N_{m}}$. The “conditional” KL representation of $Y(\mathbf{x},\omega)$ has zero variance at $\{\mathbf{x}^{(i)}\}_{i=1}^{N_{m}}$ (for simplicity, we assume that the measurements are exact) and increases away from the measurement locations. In this section, we demonstrate that the dimensionality of conditional KL expansion is smaller than that of the corresponding unconditional KL expansion. In the following two sections, we introduce the unconditional and conditional KL expansions of $Y(\mathbf{x},\omega)$.

3.1 “Unconditional” Karhunen-Loève representation of $Y(\mathbf{x},\omega)$

For a mean-square continuous stochastic process $Y(\mathbf{x},\omega)=\mathbb{E}[Y(\mathbf{x},\omega)]+Y^{\prime}(\mathbf{x},\omega)$ with the mean (expectation with respect to $\omega$) $\overline{Y}(\mathbf{x})=\mathbb{E}[Y(\mathbf{x},\omega)]$, zero-mean fluctuations $Y^{\prime}(\mathbf{x},\omega)$, and “unconditional” covariance function $C_{Y}(\mathbf{x},\mathbf{y}):=\mathbb{E}[Y^{\prime}(\mathbf{x},\omega)Y^{\prime}(\mathbf{y},\omega)]$, $\mathbf{x},\mathbf{y}\in D$, the KL expansion is given by Mercer’s theorem [22] as

[TABLE]

where $\lambda_{n}$ are positive eigenvalues and ${\epsilon}_{n}(\mathbf{x})$ are mutually orthogonal eigenfunctions, i.e., $\int_{D}{\epsilon}_{n}(\mathbf{x}){\epsilon}_{l}(\mathbf{x})d\mathbf{x}=\delta_{n,l}$ and $\delta_{n,l}$ is the Kronecker delta function. The eigenfunctions are found as the solution of the Fredholm integral equation of the second kind:

[TABLE]

The random variables $\xi_{n}$ are mutually uncorrelated and have zero mean and unit variance. Moreover, if $Y(\mathbf{x},\omega)$ is a Gaussian process, then $\xi_{n}$ are independent standard normal random variables. By convention, the eigenvalues $\lambda_{n}$ in the KL expansion are arranged in the decreasing order and the truncated $N$-term KL expansion of $Y(\mathbf{x},\omega)$ is defined as:

[TABLE]

3.2 “Conditional” Karhunen-Loève representation of $Y(\mathbf{x},\omega)$

The observations $\{Y(\mathbf{x}^{(i)})=\ln\kappa(\mathbf{x}^{(i)})\}_{i=1}^{N_{m}}$ are used to compute the (unconditional) stationary covariance describing an unconditional random field $Y$ through the variagram analysis or by maximizing the likelihood function[5]. In addition, the observations can be used to model a conditional random field using GPR as presented in [25]. In the following, we will use $\mathbf{x}^{*}$ to denote the set of $N_{m}$ observation locations $\{\mathbf{x}^{(i)}\}_{i=1}^{N_{m}}$ of the random field and $Y(\mathbf{x}^{*})$ to denote the column vector of $N_{m}$ observations $\{Y(\mathbf{x}^{(i)})\}_{i=1}^{N_{m}}$. Next, we define the covariance matrix of the $Y$ observations:

[TABLE]

and $R$ as the matrix with $n$th row given by the values of the $n$th eigenfunction at the observed locations $\mathbf{x}^{*}$,

[TABLE]

Next, we approximate $Y(\mathbf{x},{\omega})$ with the truncated KL expansion with $N_{G}$ terms $Y(\mathbf{x},{\omega})\approx Y_{N_{G}}(\mathbf{x},{\omega})=\overline{Y}(\mathbf{x})+\sum_{n=1}^{N_{G}}\sqrt{\lambda_{n}}{\epsilon}_{n}(\mathbf{x})\xi_{n}$. Then, $R$ is a $N_{G}\times N_{m}$ matrix and $\Sigma=R^{T}\Lambda R$, where $\Lambda$ is the diagonal matrix with $(\lambda_{1},\dots,\lambda_{N_{G}})^{T}$ being the positive diagonal vector.

In [25], it was demonstrated that $Y_{N_{G}}(\mathbf{x},{\omega})$ conditioned on $Y(\mathbf{x}^{*})$, $\widetilde{Y}_{N_{G}}(\mathbf{x},{\omega}):=[Y_{N_{G}}(\mathbf{x},{\omega})|Y(\mathbf{x}^{*})]$, can be written in term of $\xi_{n}$ conditioned on $Y(\mathbf{x}^{*})$, $\widetilde{\xi_{n}}:=[\xi_{n}|Y(\mathbf{x}^{*})]$, $n=1,...,N_{G}$. Noting that the covariance between $\xi_{n}$ and observation $Y(\mathbf{x}^{(i)})$ is:

[TABLE]

the conditional mean $\tilde{\mu}$ and covariance $\widetilde{M}=(\widetilde{m}_{n,k})\in\mathbb{R}^{N_{G}\times N_{G}}$ of $\{\widetilde{\xi}_{n}\}_{n=1}^{N_{G}}$ is obtained in the form [28]:

[TABLE]

Then, $\{\widetilde{\xi}_{n}\}_{n=1}^{N_{G}}$ has the same distribution as

[TABLE]

where $\tilde{\bm{\mu}}=\Lambda^{1/2}R\Sigma^{-1}({Y}(\mathbf{x}^{*})-\overline{Y(\mathbf{x}^{*})})$, $\widetilde{M}=I-\Lambda^{1/2}R\Sigma^{-1}R^{T}\Lambda^{1/2}$, and $\{\eta_{k}\}_{k=1}^{N_{G}}$ is a sequence of $i.i.d.$ $\mathcal{N}(0,1)$ random variables. As a result, $\tilde{Y}_{N_{G}}(\mathbf{x},\omega)$ can be written in terms of $\{\eta_{k}\}_{k=1}^{N_{G}}$ as

[TABLE]

Note that $\Sigma$ should be a full rank matrix and $rank(\Sigma)=N_{m}$. If $\Sigma$ is not full rank, then it is possible to select a subset of the $N_{m}$ measurements for which $\Sigma$ would be full rank.

Below, we prove a lemma and a theorem that the rank of $\widetilde{M}$ and the dimension of $\widetilde{Y}_{N_{G}}(\mathbf{x},{\omega})$ could be reduced to $N_{G}-N_{M}$.

Lemma 1**.**

If $rank(\Sigma)=N_{m}$, then $rank(\widetilde{M})=N_{G}-N_{m}$.

Proof.

By definition, it is easy to verify that $M^{2}=M$ and $M^{T}=M$. Since $M-M^{2}=M(I-M)=0$, the dimension of the null space of $M$ equals to the dimension of the range space of $I-M$ which is the same as the rank of $I-M$, i.e., $dim(ker(M))=rank(I-M)$. Since $I-M=\Lambda^{1/2}R\Sigma^{-1}R^{T}\Lambda^{1/2}$, then $rank(I-M)\leq rank(\Sigma)=N_{m}$.

Also,

[TABLE]

Then, $dim(ker(M))\geq rank(\Lambda^{1/2}R)$. Since $\Sigma=R^{T}\Lambda R=(\Lambda^{1/2}R)^{T}(\Lambda^{1/2}R)$ and $rank(\Sigma)=N_{m}$, then $rank(\Lambda^{1/2}R)\geq N_{m}$. Thus, we have $dim(ker(M))\geq N_{m}$. In combination with the previous result, we get $dim(ker(M))=N_{m}$ and $rank(M)=N_{G}-dim(ker(M))=N_{G}-N_{m}$. ∎

Theorem 1**.**

Assume ${\epsilon}_{n}(\mathbf{x})$ are orthonormal functions, $\lambda_{n}$ are positive values, $\xi_{n}$ are i.i.d standard normal random variables for $n=1,\dots,N_{G}$ and $rank(\Sigma)=N_{m}$ with $0<N_{m}<N_{G}$, then the conditional random field $\tilde{Y}_{N_{G}}(\mathbf{x},{\omega})$ can be represented by an expansion of $N_{G}-N_{m}$ i.i.d. standard normal random variables.

Proof.

From Eq (11), the covariance matrix of $\widetilde{Y}_{N_{G}}(\mathbf{x},{\omega})$, $C_{\tilde{Y}}$, can be written as

[TABLE]

Using Lemma 1, we write $rank(Cov_{\tilde{Y}})=rank(\widetilde{M})=N_{G}-N_{m}$. Consequently, there are only $N_{G}-N_{m}$ nonzero eigenvalues in the KL expansion of $\tilde{Y}_{N_{G}}(\mathbf{x},{\omega})$ that can be expressed as

[TABLE]

where $\tilde{\lambda}_{i}$ and $\tilde{{\epsilon}}_{i}(\mathbf{x})$, $i=1,\dots,N_{G}-N_{m}$ are eigenvalues and eigenfunctions of $Cov_{\tilde{Y}}$ respectively, and $\xi_{i}$, $i=1,\dots,N_{G}-N_{m}$ are i.i.d. standard Gaussian random variables. ∎

Remark 1**.**

The dimensionality of the random field $\tilde{Y}_{N_{G}}(\mathbf{x},{\omega})$ is $N_{G}-N_{m}$ and is smaller than that of the original unconditioned random field, i.e., $N_{G}$. Conditioning on observation reduces the complexity of the forward problem (1). In addition to the dimensionality reduction, conditioning also reduces the variance of $\tilde{Y}_{N_{G}}(\mathbf{x},{\omega})$ [25].

4 Conditional Stochastic Collocation Method

With the conditioned random conductivity $\tilde{\kappa}(\mathbf{x},{\omega})=\exp[\tilde{Y}_{N_{G}}(\mathbf{x},{\omega})]$, Eq (2) has $N_{G}-N_{m}$ random dimensions and can be efficiently solved with the gPC numerical methods when $N_{G}-N_{m}$ is relatively small. Given that the computational cost of gPC-based methods exponentially increases with the number of random dimensions, the conditioning significantly reduces the computational cost of gPC. Here, we assume that $N_{G}-N_{m}$ is small (i.e., the number of observations is comparable to $N_{G}$), and use the conditional stochastic collocation method [33] as a surrogate for Eq (2).

4.1 Generalized polynomial chaos

Based on Eq (13), the conditional random field $\tilde{\kappa}(\bm{x},\omega)$ can be represented by KL expansion with $N_{G}-N_{m}$ terms as

[TABLE]

and the “conditional” $u$, $\tilde{u}$, satisfies

[TABLE]

The gPC method for Eq (15) is based on an orthogonal polynomial approximation of $\tilde{u}(\mathbf{x},\bm{\xi})$. Let $\mathbf{i}=(i_{1},\dots,i_{N})\in\mathbf{N}_{0}^{n}$ be a multi-index with $|\mathbf{i}|=i_{1}+\dots+i_{N}$ and $P\geq 0$ be an integer. Then, the $P$th-degree gPC expansion of function $\tilde{u}(\mathbf{x},\bm{\xi})$ is defined as

[TABLE]

where

[TABLE]

are the coefficients of the expansion, $\Phi_{\mathbf{i}}(\bm{\xi})$ are the basis functions

[TABLE]

and $\bm{\xi}=(\xi_{1},\dots,\xi_{N_{G}-N_{m}})$ is an $(N_{G}-N_{m})$-dimensional random vector. Here, $\phi_{j}(\xi_{k})$ is the $j$th-degree one-dimensional orthogonal polynomial in $\xi_{k}$ direction satisfying

[TABLE]

For identical independent distributed (i.i.d.) Gaussian random variables $\{\xi_{i}\}_{i=1}^{N}$, $\Phi_{\mathbf{i}}(\bm{\xi})$ are Hermite polynomials. The number of $\Phi_{\mathbf{i}}(\bm{\xi})$ is (\begin{array}[]{c}N_{G}-N_{m}+P\\ P\end{array}) [34, 31]. The classical approximation theory guarantees that the gPC approximation (16) converges to $\tilde{u}(\mathbf{x},\bm{\xi})$ in $L^{2}$-norm as the degree $P$ increases when $\tilde{u}(\mathbf{x},\bm{\xi})$ is square integrable with respect to the probability measure.

4.2 Stochastic collocation method

In the stochastic collocation method, the expansion coefficients ${c}_{\mathbf{i}}(\mathbf{x})$ are approximated as

[TABLE]

where $\{\bm{\xi}^{(m)}\}_{m=1}^{M}$ is a set of quadrature points and $w^{(m)},m=1,\dots,M$ are the corresponding weights. For each collocation point $\bm{\xi}^{(i)}$, $u({\mathbf{x},\bm{\xi}^{(i)}})$ is obtained by solving (15) with $\tilde{\kappa}(\mathbf{x},\bm{\xi})$ replaced by $\tilde{\kappa}(\mathbf{x},\bm{\xi}^{(i)})$.

Then,

[TABLE]

and the mean and variance of $\tilde{u}(\mathbf{x},\bm{\xi})$ can be approximated as

[TABLE]

and

[TABLE]

There are various quadrature rules, including tensor product quadrature rules for low-dimensional $\bm{\xi}$ and the sparse grid methods for moderately dimensional $\bm{\xi}$ [33, 24, 17, 13]. If the solution allows a low-dimensional representation, the compressed sensing can be used to decrease the number of collocation points [35, 36, 11]. The convergence analysis of stochastic collocation methods is given in [33, 30].

Remark 2**.**

In order to solve optimization problems, including (4),(21), and (24) (or to make the optimization problems well posed), additional measurements on solution (state) are required. Function $u(x,\omega)$ with random inputs $\omega$ and variance $\sigma^{2}_{u}(x)$ is more sensitive with respect to $\omega$ at position $x$ where the variance is larger. Therefore, the optimal location of $u(x,\omega)$ measurements for solving optimization problems should collocate with the local maxima of $\sigma^{2}_{u}(x)$. In practice, only a few measurements are affordable; therefore, determining the optional measurement locations for the optimization problem is very important.

5 Bayesian estimation of $\tilde{\kappa}(\mathbf{x},\bm{\xi})$

We assume that in addition to $\{\hat{\kappa}(\mathbf{x}^{j})\}_{j=1}^{N_{m}}$ measurements of $\hat{\kappa}$, we have $\{\hat{u}(\mathbf{x}^{j})\}_{j=1}^{N_{k}}$ measurements of $\hat{u}$. With the conditional KL model (14), we approximate $\hat{\kappa}(\bm{x})$ as $\hat{\kappa}(\bm{x})\approx\tilde{\kappa}(\bm{x},\bm{\xi}^{*})$, where $\bm{\xi}^{*}$ is the $(N_{G}-N_{m})$ dimensional deterministic vector. Then, the estimation of $\hat{\kappa}(\bm{x})$ reduces to the estimation of $\bm{\xi}^{*}$. With the surrogate model $\tilde{u}_{C}(\mathbf{x},\bm{\xi})$, we estimate $\bm{\xi}^{*}$ by solving the following minimization problem:

[TABLE]

where the regularization term $\|\bm{\xi}-\bm{\xi}^{o}\|^{2}$ was proposed in [27] and $\bm{\xi}^{o}\in\mathbb{R}^{N_{G}-N_{m}}$ and $\lambda$ are the regularization parameters. Solving this optimization problem using standard minimization methods could be computationally expensive for large $N_{G}-N_{m}$. Therefore, we employ a Bayesian statistical approach to solve the optimization problem and provide the correspondence between the regularization parameters and the prior statistics of $\bm{\xi}$.

In Eq (21), we assume that the errors $\{\hat{u}(\mathbf{x}^{j})-\tilde{u}_{C}(\mathbf{x}^{j},\bm{\xi})\}_{j=1}^{N_{k}}$ are i.i.d. random variables, that is,

[TABLE]

where $\delta_{j}$ $j=1,\dots,N_{k}$ is i.i.d. random variable with a certain density $p_{\delta}$. Following a common practice in parameter estimation, we assume that $p_{\delta}$ is a Gaussian density with mean [math] and standard deviation $\sigma_{\delta}\ll 1$. Then, the likelihood takes the form

[TABLE]

With $p_{\bm{\xi}}$ denoting the prior probability density for $\bm{\xi}$, we use Bayes’ rule to obtain a posterior probability density for $\bm{\xi}$ given the observations $\{\hat{u}(x^{(j)})\}_{j=1}^{N_{k}}$:

[TABLE]

The (prior) density $p_{\bm{\xi}}$ of $\bm{\xi}$ is usually chosen to be Gaussian, that is, $\bm{\xi}\sim\mathcal{N}(\bm{\xi}^{o},\theta\mathbf{I})$ with the hyper-parameter $\theta$ controlling the prior variance and $\bm{\xi}^{o}$ being the prior mean. Usually, the prior mean of $\bm{\xi}$ is set to $\bm{0}$. The posterior density of $\bm{\xi}$ with hyper-parameter is rewritten as,

[TABLE]

Then, the maximum a posteriori (MAP) estimator of this posterior distribution is

[TABLE]

Note that (24) is a special case of the minimization problem (21) with $\lambda=\frac{\sigma_{\delta}}{\theta}$.

To generate posterior distributed sampling, we employ the widely used Differential Evolution Adaptive Metropolis (DREAM) MCMC toolbox developed by Vrugt[29]. Given the prior $\bm{\xi}_{0}=\bm{0}$, and $\theta=1$ such that $\bm{\xi}$ is consistent to the form of Gaussian process (13), we generate the samples of $\bm{\xi}$ for the MCMC. For $\sigma_{\delta}$, we assume the errors are small, that is, we choose $\sigma_{\delta}$ to be about $1e-3$. Whenever a sample $\bm{\xi}$ is generated, the surrogate $\tilde{u}_{C}(\bm{x}^{i},\bm{\xi})$ is evaluated and the errors $|\hat{u}(\mathbf{x}^{i})-\tilde{u}_{C}(\mathbf{x}^{i},\bm{\xi})|^{2}$, $i=1,N_{k}$ in the Metropolis Hastings algorithm are computed. Then, the DREAM MCMC toolbox is employed to generate samples of posterior distribution and the MAP probability estimate of this posterior distribution is used as an approximate solution of the optimization problem (21). The numerical examples in Section 6 show the accuracy of this approximation.

5.1 Determining measurement locations of the state $u(\mathbf{x})$

The choice of $u(\mathbf{x})$ measurements $\{\hat{u}(\mathbf{x}^{j})\}_{j=1}^{N_{k}}$ can impact the posterior distributions obtained through MCMC. Because $\sigma^{2}_{u}(\mathbf{x})$, the variance of $u(\mathbf{x},\bm{\xi})$, describes the variability of $u(\mathbf{x},\bm{\xi})$ with respect to the parameters $\bm{\xi}$, we propose the following strategy:

• 1

Approximate $\sigma^{2}_{u}(\mathbf{x})$ by $\tilde{\sigma}^{2}_{u}(\mathbf{x})=\sum_{|\mathbf{i}|=1}^{P}\tilde{c}^{2}_{\mathbf{i}}(\mathbf{x})$.

• 2

Find all local maxima and saddle points on the $\tilde{\sigma}^{2}_{u}({\mathbf{x}})$ denoted by $S:=\{\tilde{\sigma}^{2}_{u}({\mathbf{x_{*}^{j}}})\}_{j=1}^{N_{\ell}}$ in descending order.

• 3

If $N_{\ell}\geq N_{k}$, we choose the first $N_{k}$ locations from ordered set $S$.

• 4

If ${N_{\ell}}<N_{k}$, divide the physical domain $D$ in $N_{k}$ equal blocks. For blocks without saddle points from $S$, find $\sigma^{2}_{u}(\mathbf{x})$ maxima in this block, arrange them in descending order, and choose the first $N_{k}-N_{\ell}$ locations to add to $S$.

The numerical examples show that this strategy outperforms equal or random distribution of $u(\bm{x})$ measurements.

6 Numerical examples

We implement the proposed parameter estimation approach for one- and two-dimensional steady-state diffusion equations. We assume that the “reference” diffusion coefficient is an instance of the known Gaussian process that can be accurately represented with the the $N_{G}$-dimensional KL expansion (7). As a result, the reference $\kappa$ filed lies in the space of the conditional KL expansion (14).

6.1 One-dimensional example

Consider the following one-dimensional steady-state diffusion equation with the unknown coefficient $\hat{\kappa}(x)$ modeled as the random field $\kappa(x,\omega)$:

[TABLE]

where $u_{l}=0$ and $u_{r}=2$. We assume that $\kappa(x,\omega)$ has lognormal distribution (i.e, $g(x,\omega)=\log(\kappa(x,\omega))$ is Gaussian) with mean $\mu_{k}=5.0$, variance $\sigma_{k}=2.5$ and the covariance function

[TABLE]

where the correlation length is $L=0.05$ (i.e., $L$ is 20 times smaller than the domain size).

We find that the finite KL expansion with $N_{G}=25$ terms,

[TABLE]

captures $95\%$ of the spectrum of this field covariance function. Here, the mean $\mu_{g}$ and standard deviation $\sigma_{g}$ of $g(x,\omega)$ are given in terms of $\mu_{g}$ and $\sigma_{g}$ as

[TABLE]

Then, we compute the reference conductivity field ${\hat{\kappa}}(x)$ as a realization of Eq (27) and solve Eq (25) with the finite element method to obtain the reference solution, $\hat{u}(x)$. Next, we select $N_{\kappa}$ values of $\hat{\kappa}(x)$ and $N_{u}$ values of $\hat{u}(x)$ as measurements of these fields and reconstruct $\hat{\kappa}(x)$ as $\tilde{\kappa}(x,\bm{\xi}^{*})$ by solving the minimization problem (24). The conditioned surrogate model $\tilde{u}_{C}$ in (24) is constructed by solving Eq (25) with the finite element method at collocation points. In the finite element solutions, the domain is discretized with $256$ equal size elements. We study the relative error of the inferred conductivity $\tilde{\kappa}(x,\bm{\xi}^{*})$,

[TABLE]

as a function of $N_{\kappa}$, $N_{u}$, and the sampling strategy.

We consider three $\hat{\kappa}(x)$ sampling strategies: Case 1, the observation locations of $\hat{\kappa}(x)$ are chosen randomly; Case 2, the observation locations are equally distributed; and Case 3, the observation locations coincide with locations of local minima and maxima. The referenced field $\kappa(x)$ and the locations of observations for each case are presented in Figure 1a.

We also study the effect of $\hat{u}$ measurement locations on the relative error $\varepsilon(x)$. As suggested in Remark 3, $u(x,\bm{\xi})$ is more sensitive with respect to the random inputs $\bm{\xi}$ where its variance is large. Therefore, it reasonable to assume that the best candidate locations of $\hat{u}$ measurements are those where $\sigma^{2}_{u}(x)$ has local maxima. The $u$ variance $\sigma^{2}_{u}(x)$ can be evaluated with the conditional gPC method using Eq (20). On the other hand, the MCMC optimization algorithm requires $u$ measurements to be uncorrelated. Therefore, we compare the accuracy of the $u$ sampling strategy outlined in Section 5.1 (based on $\sigma^{2}_{u}(x)$ local maxima) with uniform and random sampling strategies.

6.1.1 Case 1: Random-chosen locations of the conductivity measurements

Here we assume that no expert knowledge about the conductivity field is available and the 20 sampling locations of $k$ are chosen randomly, as shown in Figure 1 by red square markers. It can be seen that $\hat{\kappa}$ oscillates significantly between the sampling points and cannot be accurately learned via regression only, as shown in Figure 1b.

To chose the $u$ measurement locations, we compute the variance of $u(x,\bm{\xi})$ using (20) (shown in Figure 2) and then perform the following three tests:

• (a)

We randomly choose six locations of the reference solution $\hat{u}$ for the optimization step as shown in Figure 2 by red “$\square$” markers.

• (b)

We choose six equally distributed locations of $\hat{u}$, as shown in Figure 2 by blue “$\circ$” markers.

• (c)

We choose six locations using algorithm in Section 5.1 as shown in Figure 2 by green “$+$” markers.

Figure 3(a) shows the reference and estimated conductivity field. The associated $\varepsilon(x)$ errors are presented in Figure 3(b). We can see that in all three tests, the errors are largest between the two cross markers where there are no $\hat{\kappa}$ measurements. On the other hand, optimally selected locations of $\hat{u}$ in test (c) improve parameter estimation, i.e., reduce the relative errors $\|\varepsilon\|_{L^{\infty}}$ from more than $60\%$ (for randomly or equally distributed $\hat{u}$ measurements) to less than $14\%$.

6.1.2 Case 2: Observations of the conductivity coefficient are equally distributed

Next, we consider $\kappa$ observation locations equally distributed, as shown in Figure 1a by blue “$\circ$” markers. We test the same three strategies for selecting $u$ measurement locations as in Case 1, including randomly and uniformly distributed locations and locations chosen based on the local maxima of $\sigma^{2}_{u}(x)$. The resulting distributions of the $u$ measurements locations as well as $\sigma^{2}_{u}(x)$ are shown in Figure 4.

• (a)

We randomly choose six locations for the measurements of $\hat{u}$, as shown in Figure 4 by red $\square$ markers.

• (b)

We choose six equally distributed locations, as denoted in Figure 4 by blue $\circ$ markers.

• (c)

We choose six locations using algorithm in Section 5.1, as shown in Figure 4 by green $+$ markers.

Figure 5a shows the reference and estimated conductivity fields and Figure 5b presents the associated relative errors $\varepsilon(x)$. Here, the relative errors in $\hat{\kappa}(x)$ are less than $5\%$ for all considered choices of $u$ measurement locations. The comparison of relative errors in Figure 5(right) indicates that the selection of the $\hat{u}$ measurement locations based on $\sigma^{2}_{u}(x)$ has the smallest $\|\varepsilon(x)\|_{\infty}$ error.

6.1.3 Case 3: Observations of $\kappa$ are colocated with local maxima and minima of $\hat{k}(x)$

Here, we assume that the locations of local maxima and minima of $\hat{\kappa}(x)$ are known (e.g., based on expert knowledge). Then, we take observations at all $14$ local minima and local maxima and one inflection point of $\hat{\kappa}(x)$. Another five observations of $\hat{\kappa}(x)$ are taken at random locations. The selected observations are shown by green “$+$” markers in Figure 1a.

We test the same three strategies for selecting the measurement locations of $\hat{u}(x,\xi)$ as in Cases 1 and 2. Figure 6 shows the conditional variance of $u$ and $u$ measurement locations selected randomly, uniformly, and according to the conditional variance of $u$. Figure 7 shows the reference and estimated $\kappa$ as well as the corresponding errors. For all considered $u$ sampling strategies, $\varepsilon(x)$ are less than $1.2\%$. Among them, the $u$ measurement locations selected according to $\sigma^{2}_{u}(x)$ give the smallest $\|\varepsilon(x)\|_{\infty}$, which is less than $0.6\%$.

6.2 Two-dimensional example with large correlation length of $\kappa$ (smooth $\kappa$ field)

Here, we consider the two-dimensional diffusion equation with the partially observed coefficient $\kappa$ modeled as a random field $\kappa(\mathbf{x},\omega)$,

[TABLE]

subject to the boundary conditions

[TABLE]

where $D=[0,240]\times[0,60]$ and $\mathbf{n}=(0,1)$.

We assume that $\kappa(\mathbf{x},\bm{\xi})$ is a lognormal field with mean $\mu_{\kappa}=5$, standard deviation $\sigma_{\kappa}=2.5$, and $g(\mathbf{x},\bm{\xi})=\log\kappa(x,\bm{\xi})$ is given by the truncated KL expansion:

[TABLE]

where $\mu_{g}$ and $\sigma_{g}$ are computed by (28), and $(\lambda_{i},\epsilon_{i}(\mathbf{x}))$ are the first 25 eigenpairs that in total capture 95% of spectra of the exponential correlation function

[TABLE]

Correlation lengths in the $x_{1}$ and $x_{2}$ directions are $L_{1}=240$ and $L_{2}=100$, which are of the order of the domain size and result in a relatively smoothly varying $\kappa(\mathbf{x},\bm{\xi})$ field shown in Figure 8a.

We choose one realization $\hat{\kappa}(\mathbf{x})=\kappa(\mathbf{x},\bm{\xi}^{*})$ as the reference diffusion coefficient and $\hat{u}(\mathbf{x},\bm{\xi}^{*})$, the corresponding solution of Eqs (30) and (31), as the reference hydraulic head. Finite volume method with uniform mesh ($80\times 20$ elements) is adopted to solve Eqs (30) and (31). We consider two sets of $\hat{\kappa}$ observations, inlcuding a randomly sampled $\hat{\kappa}$ and $\hat{\kappa}$ sampled based on the expert knowledge of local maxima and minima of $\hat{\kappa}$. For each distribution of $\hat{\kappa}$ measurements, we consider two sets of $\hat{u}$ measurements, randomly located and chosen based on $\sigma^{2}_{u}(\mathbf{x})$ as described in Section 5.1. The relative error $\varepsilon$ in the estimated $\hat{\kappa}$ are computed for each case using Eq (29).

6.2.1 Case 1: Randomly distributed locations of $\kappa$ observations

Here, we assume that 20 observations of $\hat{\kappa}$ are randomly distributed in the domain, as shown in Figure 8a. The resulting conditional KL expansion of $\kappa$ has 5 unknown parameters that we estimate using 10 measurements of $\hat{u}$. To study the effect of $\hat{u}$ measurement locations on the $\hat{\kappa}$ estimate errors, we consider measurements of $\hat{u}$ distributed (a) randomly and (b) according to $\sigma^{2}_{u}(\mathbf{x})$. Figure 8b depicts the conditional $\sigma^{2}_{u}$t and the $u$ measurements locations. Figure 9 presents the corresponding relative errors of the estimated ${\hat{\kappa}}$.

In both cases, the relative errors are smaller than 1%. That is no surprise given that the number of measurements of $\kappa$ (20) is close to the number of terms in the unconditional KL representation of $\kappa$. Choosing the location of $\hat{u}$ measurements based on the variance of $u$ reduces $\varepsilon(\mathbf{x})$ by approximately the factor of 3 and the infinity norm, $\|\varepsilon(\mathbf{x})\|_{L^{\infty}}$, by 80% compared to the randomly chosen locations of $\hat{u}$ measurements.

6.2.2 Case 2: Some observations of $\kappa$ are colocated with local maxima and minima of $\kappa$ field

Here, we assume that 9 $\hat{\kappa}$ observations are collocated with local maxima and minima of $\hat{\kappa}$, and the other 11 $\hat{\kappa}$ observations are randomly spread across the domain, as shown in Figure 10a. The resulting conditional KL expansion of $\kappa$ has 5 unknown parameters that we estimate using 10 measurements of $\hat{u}$. As in the previous example, we consider measurements of $\hat{u}$ distributed randomly and according to $\sigma^{2}_{u}(\mathbf{x})$, as shown in Figure 10b. The corresponding relative errors in the estimated $\hat{\kappa}$ field are shown in Figure 11.

We see that for both choices of the $\hat{u}$ measurement locations, $\varepsilon(x)$ is smaller than 1%. Choosing the location of $\hat{u}$ measurements based $\sigma^{2}_{u}(\mathbf{x})$ reduces $\|\varepsilon\|_{L_{\infty}}$ by $25\%$ as compared with the random distribution of $\hat{u}$ measurements. The comparison of Figures 9 and 11 show that the proposed parameter estimation method performs well for both cases.

6.3 Two-dimensional diffusion equation with small correlation length of $\kappa$ (rough $\kappa$ field)

In this example, we assume the exponential covariance function of $g(x,\omega)=\log\kappa(x,\omega)$:

[TABLE]

where $|\mathbf{x}^{(1)}-\mathbf{x}^{(2)}|$ is the Euclidean distance between $\mathbf{x}^{(1)}$ and $\mathbf{x}^{(2)}$ and $L=0.1$. The mean and variance of $\kappa(\mathbf{x},\omega)$ are set to $\mu_{k}=5$ and $\sigma_{k}=2.5$.

We solve the two-dimensional equation (30) on $D=[0,2]\times[0,1]$ subject to the boundary conditions:

[TABLE]

where $\mathbf{n}=(0,1)$.

Here, the ratio of the correlation length to the domain size is on the order of 0.1, which is smaller than in the previous case and results in a rougher $\kappa$. The KL representation of the unconditional random field $g(\mathbf{x},\omega)$ requires $210$ terms to capture 95% of spectrum,

[TABLE]

The reference conductivity field $\kappa(\bm{x},\mathbf{\xi}^{*})$ is constructed as a realization of $\kappa(\bm{x},\mathbf{\xi})$ and is shown in Figure 12a, which is significantly rougher than the field in Figure 8a considered in the previous case.

The corresponding reference field $\hat{u}(\bm{x},\mathbf{\xi}^{*})$ is found as the solution of Eq (30) subject to the boundary conditions (34). This equation is solved using the finite volume method with equal-distanced mesh with $128\times 64$ elements. We assume that 205 measurements of $\hat{\kappa}$ and 10 measurements of $\hat{u}$ are available to estimate the entire field $\hat{\kappa}(x)$. The KL expansion (35), conditioned on 205 measurements, has five random dimensions. As before, we study the effect of $\hat{\kappa}$ and $\hat{u}$ measurements locations distribution on the accuracy of $\kappa$ estimation.

6.3.1 Case 1: Randomly located $\kappa$ observations

Here, we assume that the locations of $\hat{\kappa}$ measurements are randomly distributed, as shown in Figure 12a. Next we consider 10 $\hat{u}$ measurement locations distributed in two ways, including randomly distributed measurements and measurements distributed according to the conditional $\sigma^{2}_{u}(\mathbf{x})$; both distributions are depicted in Figure 12b. Figure 13 presents the corresponding relative error of inferred $\hat{\kappa}(x)$.

The comparison of Figures 13a and b show that choosing $\hat{u}$ measurement locations according to the $u$ variance reduces the $\|\varepsilon\|_{L^{\infty}}$ error by a factor of 10.

6.3.2 Case 2: Some observations of $\kappa$ are colocated with local maxima and minima of $\kappa$ field

Here, we assume that 50 locations of $\hat{\kappa}$ measurements are collocated with local maxima, minima, and saddle points of $\hat{\kappa}$ and, $155$ measurement locations are randomly distributed, as shown in Figure 14a. As before, we consider two choices of $\hat{u}$ measurement locations, including random locations and locations based on the conditional $\sigma^{2}_{u}(\mathbf{x})$ as shown in Figure 14b. Figure 14b also displays the conditional $\sigma^{2}_{u}(\mathbf{x})$ as a function of $\mathbf{x}$. Figure 15 presents the corresponding relative error of the inferred $\hat{\kappa}(\mathbf{x})$.

The comparison of Cases 1 and 2 shows that for randomly located $\hat{\kappa}$ measurements, the selection of $u$ measurement locations based on the conditional $u$ variance is very important – it reduces the $\|\varepsilon(\mathbf{x})\|_{L_{\infty}}$ error by a factor of 10. For specially selected locations of $\hat{\kappa}$, the $\sigma^{2}_{u}(\mathbf{x})$-based selection of $\hat{u}$ measurement locations reduces $\|\varepsilon(\mathbf{x})\|_{L_{\infty}}$ by a factor of 0.4.

Remark 3**.**

In all considered examples, we see a correlation between the error of the estimated $\kappa$ and the variance of $u$. That is, near the area that the gPC solution has large variance, the inferred $\kappa$ usually has larger relative error.

7 Conclusion

In this work, we proposed a conditional KL expansion and gPC surrogate model for estimating space-dependent coefficients in partial differential equation models using measurements of the coefficient and state variable. We demonstrated that the conditional gPC model simplifies the optimization problem and reduces the number of unknown parameters as compared to parameter estimation using traditional (unconditional) gPC surrogates. Specifically, the conditional gPC reduces dimensionality of the parameter space and replaces expensive direct solutions of PDEs with the conditional gPC surrogate. Also, in the absence of measurement error, the estimated parameter field exactly matches the parameter measurements. Furthermore, we proposed using the surrogate model to determine measurement locations based on the variance of state variable conditioned on the coefficient measurements. Specifically, we proposed collocating measurement with the local maxima of the state variable variance predicted by the conditional gPC. We presented one- and two-dimensional examples demonstrating the overall accuracy of the proposed approach (it is more accurate than Gaussian process regression). We also showed that the error in the estimated parameter is smaller when the state variable measurements are chosen based on its variance rather than uniformly or randomly collected. Also, we demonstrated that selecting the coefficient measurement locations based on physical knowledge (i.e., collocating coefficient measurements with local maxima and minima of the coefficient) further reduces error in the estimated coefficient.

Finally, in all considered examples, we see correlation between the error of the estimated coefficient and the variance of the solution; that is, near the area that the gPC solution has large variance, the inferred coefficient has larger relative error. Therefore, the conditional variance of state variables can be used to guide the data acquisition for the coefficient.

In the current work, we only investigated cases where the unknown coefficient lies in a function space of the finite KL modes (i.e., the reference solution is a realization of the known finite-dimensional Gaussian process) and the number of solution measurements is greater than the dimension of the conditional KL representation. Problems with small data sets and coefficients with unknown and/or non-Gaussian distribution will yield additional errors and uncertainty in parameter estimation and will be subject of our future work.

Acknowledgements

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