Posterior Consistency for Gaussian Process Approximations of Bayesian Posterior Distributions

Introduction

Given a mathematical model of a physical process, we are interested in the inverse problem of determining the inputs to the model given some noisy observations related to the model outputs. Adopting a Bayesian approach , we incorporate our prior knowledge of the inputs into a probability distribution, referred to as the prior distribution, and obtain a more accurate representation of the model inputs in the posterior distribution, which results from conditioning the prior distribution on the observations. Since the posterior distribution is generally intractable, sampling methods such as Markov chain Monte Carlo (MCMC) are typically used to explore it. A major challenge in the application of MCMC methods to problems of practical interest is the large computational cost associated with numerically solving the mathematical model for a given set of the input parameters. Since the generation of each sample by the MCMC method requires a solve of the governing equations, and often millions of samples are required, this process can quickly become very costly.

This drawback of fully Bayesian inference for complex models was recognised several decades ago in the statistics literature, and resulted in key papers which had a profound influence on methodology . These papers advocated the use of a Gaussian process surrogate model to approximate the solution of the governing equations, and in particular the likelihood, at a much lower computational cost. This approximation then results in an approximate posterior distribution, which can be sampled more cheaply using MCMC. However, despite the widespread adoption of the methodology, there has been little analysis of the effect of the approximation on posterior inference. In this work, we study this issue, focussing on the use of Gaussian process emulators as surrogate models. Other choices of surrogate models such as those described in , generalised Polynomial Chaos , sparse grid collocation and adaptive subspace methods might also be studied similarly, but are not considered here. Indeed we note that the paper studied the effect, on the posterior distribution, of stochastic collocation approximation within the forward model and was one of the first papers to address such questions. That paper used the Kullback-Leibler divergence, or relative entropy, to measure the effect on the posterior, and considered finite dimensional input parameter spaces.

The main focus of this work is to analyse the error introduced in the posterior distribution by using a Gaussian process emulator as a surrogate model. The error is measured in the Hellinger distance, which is shown in to be a suitable metric for evaluation of perturbations to the posterior measure in Bayesian inverse problems, including problems with infinite dimensional input parameter spaces. We consider emulating either the parameter-to-observation map or the negative log-likelihood. The convergence results presented in this paper are of two types. In section 3, we present convergence results for simple Gaussian process emulators applied to a general function $f$ satisfying suitable regularity assumptions. In section 4, we prove bounds on the error in the posterior distribution in terms of the error in the Gaussian process emulator. The novel contributions of this work are mainly in section 4. The results in the two sections can be combined to give a final error estimate for the simple Gaussian process emulators presented in section 3. However, the error bounds derived in section 4 are much more general in the sense that they apply to any Gaussian process emulator satisfying the required assumptions. A short discussion on extensions of this work related to Gaussian process emulators used in practice is included in the conclusions in section 6.

We study three different approximations to the posterior distribution. Firstly, we consider using the mean of the Gaussian process emulator as a surrogate model, resulting in a deterministic approximation to the posterior distribution. Our second approximation is obtained by using the full Gaussian process as a surrogate model, leading to a random approximation in which case we study the second moment of the Hellinger distance between the true and the approximate posterior distribution. The uncertainty in the posterior distribution introduced in this way can be thought of representing the uncertainty in the emulator due to the finite number of function evaluations used to construct it. This uncertainty can in applications be large (or comparable) to the uncertainty present in the observations, and a user may want to take this into account to ”inflate” the variance of the posterior distribution. Finally, we construct an alternative deterministic approximation by using the full Gaussian process as surrogate model, and taking the expected value (with respect to the distribution of the surrogate) of the likelihood. It can be shown that this approximation of the likelihood is optimal in the sense that it minimises the $L^{2}$-error . In contrast to the approximation based on only the mean of the emulator, this approximation also takes into account the uncertainty of the emulator, although only in an averaged sense.

The convergence results on Gaussian process regression presented in section 3 are mainly known results from the theory of scattered data interpolation . The error bounds are given in terms of the fill distance of the design points used to construct the Gaussian process emulator, and depend in several ways on the number $K$ of input parameters we want to infer. Firstly, when looking at the error in terms of the number of design points used, rather than the fill distance of these points, the rate of convergence typically deteriorates with the number of parameters $K$. Secondly, the proof of these error estimates requires assumptions on the smoothness of the function being emulated, where the precise smoothness requirements depend on the Gaussian process emulator employed. For emulators based on Matèrn kernels , we require these maps to be in a Sobolev space $H^{s}$, where $s>K/2$. We would like to point out here that it is not necessary for the function being emulated to be in the reproducing kernel Hilbert space (or native space) of the Matèrn kernel used in order to prove convergence (cf Proposition 3.4), but that is suffices to be in a larger Sobolev space in which point evaluations are bounded linear functionals.

The remainder of this paper is organised as follows. In section 2, we set up the Bayesian inverse problem of interest. We then recall some results on Gaussian process regression in section 3. The heart of the paper is section 4, where we introduce the different approximations to the posterior and perform an error analysis. Our theoretical results are confirmed on a simple model problem in section 5, and some conclusions are finally given in section 6.

Bayesian Inverse Problems

We make the following assumption on the regularity of the parameter-to-observation map $\mathcal{G}$.

Note that in Assumption 2.2, the smoothness requirement on $\mathcal{G}$ becomes stronger as $K$ increases. The reason for this is that in order to apply the results in section 3, we require $\mathcal{G}$ to be in a Sobolev space in which point evaluations are bounded linear functionals. The second part of Assumption 2.2 is mainly included to define the constant $C_{\mathcal{G}}$, since the fact that $\sup_{u\in X}\|\mathcal{G}(u)\|$ is finite follows from the continuity of $\mathcal{G}$ and the compactness of $X$.

Gaussian Process Regression

Typical choices of the mean function $m$ include the zero function and polynomials . A family of covariance functions $k$ frequently used in applications are the Matèrn covariance functions , given by

where $\Gamma$ denotes the Gamma function, $B_{\nu}$ denotes the modified Bessel function of the second kind and $\nu,\lambda$ and $\sigma_{k}^{2}$ are positive parameters. The parameter $\lambda$ is referred to as the correlation length, and governs the length scale at which ${f_{0}}(u)$ and ${f_{0}}(u^{\prime})$ are correlated. The parameter $\sigma_{k}^{2}$ is referred to as the variance, and governs the magnitude of ${f_{0}}(u)$. Finally, the parameter $\nu$ is referred to as the smoothness parameter, and governs the regularity of ${f_{0}}$ as a function of $u$. As the limit when $\nu\rightarrow\infty$, we obtain the Gaussian covariance

Now suppose we are given data in the form of a set of distinct design points $U:=\{u^{n}\}_{n=1}^{N}\subseteq X$, together with corresponding function values

Conditioning the Gaussian process (3.1) on the known values $f(U)$, we hence obtain another Gaussian process $f_{N}$, known as the predictive process. We have

where the vector of coefficients is given by $\alpha=K(U,U)^{-1}f(U)$. Concerning the predictive covariance $k_{N}$, we note that $k_{N}(u,u)<k(u,u)$ for all $u\in X$, since $K(U,U)^{-1}$ is positive definite. Furthermore, we also note that $k_{N}(u^{n},u^{n})=0$, for $n=1,\dots,N$, since $k(u^{n},U)^{T}\;K(U,U)^{-1}\;k(u^{n},U)=k(u^{n},u^{n})$.

For stationary covariance functions $k(u,u^{\prime})=k(\|u-u^{\prime}\|)$, the predictive mean is a radial basis functions interpolant of $f$, and we can make use of results from the radial basis function literature to investigate the behaviour of $m_{N}^{f}$ and $k_{N}$ as $N\rightarrow\infty$. Before we do this, in subsection 3.2, we recall some results on native spaces (also know as reproducing kernel Hilbert spaces) in subsection 3.1.

We recall the notion of the reproducing kernel Hilbert space corresponding to the kernel $k$, usually referred to as the native space of $k$ in the radial basis function literature.

for all $u\in X$, $k(u,u^{\prime})$, as a function of its second argument, belongs to $H_{k}$,

for all $u\in X$ and $f\in H_{k}$, $\langle f,k(u,\cdot)\rangle_{H_{k}}=f(u)$.

By the Moore-Aronszajn Theorem , a unique RKHS exists for each symmetric, positive definite kernel $k$. Furthermore, this space can be constructed using Mercer’s Theorem , and it is equal to the Cameron-Martin space of the covariance operator $C$ with kernel $k$. For covariance kernels of Matèrn type, the native space is isomorphic to a Sobolev space .

Let $k_{\nu,\lambda,\sigma_{k}^{2}}$ be a Matèrn covariance kernel as defined in (3.2). Then the native space $H_{k_{\nu,\lambda,\sigma_{k}^{2}}}$ is equal to the Sobolev space $H^{\nu+K/2}(X)$ as a vector space, and the native space norm and the Sobolev norm are equivalent.

Native spaces for more general kernels, including non-stationary kernels, are analysed in . For stationary kernels, the native space can generally be characterised by the rate of decay of the Fourier transform of the kernel. The native space of the Gaussian kernel (3.3), for example, consists of functions whose Fourier transform decays exponentially, and is hence strictly contained in the space of analytic functions. Proposition 3.2 shows that as a vector space, the native space of the Matèrn kernel $k_{\nu,\lambda,\sigma_{k}^{2}}$ is fully determined by the smoothness parameter $\nu$. The parameters $\lambda$ and $\sigma_{k}^{2}$ do, however, influence the constants in the norm equivalence of the native space norm and the standard Sobolev norm.

2 Radial basis function interpolation

For stationary covariance functions $k(u,u^{\prime})=k(\|u-u^{\prime}\|)$, the predictive mean is a radial basis functions interpolant of $f$. In fact, it is the minimum norm interpolant ,

Given the set of design points $U=\{u^{n}\}_{n=1}^{N}\subseteq X$, we define the fill distance $h_{U}$, separation radius $q_{U}$ and mesh ratio $\rho_{U}$ by

The fill distance is the maximum distance any point in $X$ can be from $U$, and the separation radius is half the smallest distance between any two distinct points in $U$. The mesh ratio provides a measure of how uniformly the design points $U$ are distributed in X. We have the following theorem on the convergence of $m_{N}^{f}$ to $f$ .

for all sets $U$ with $h_{U}$ sufficiently small.

Proposition 3.3 assumes that the function $f$ is in the RKHS of the kernel $k$. Convergence estimates for a wider class of functions can be obtained using interpolation in Sobolev spaces .

for all sets $U$ with $h_{U}$ and $\rho_{U}$ sufficiently small.

We would like to point out here that in practice, it is much more informative to obtain convergence rates in terms of the number of design points $N$ rather than their associated fill distance $h_{U}$. This is of course possible in general, but the precise relation between $N$ and $h_{U}$ will depend on the specific choice of design points $U$. For uniform tensor grids $U$, the fill distance $h_{U}$ is of the order $N^{-1/K}$ (cf section 5). This suggests a strong dependence on the input dimension $K$ of the convergence rate in terms of the number of design points $N$.

Convergence of the predictive variance $k_{N}(u,u)$ follows under the assumptions of Proposition 3.3 or Proposition 3.4 using the relation in Proposition 3.5 below. This was already noted, without proof, in ; we give a proof here for completeness.

Suppose $m_{N}^{f}$ and $k_{N}$ are given by (3.6). Then

The final equality follows from the Cauchy-Schwarz inequality, which becomes an equality when the two functions considered are linearly dependent. By Definition 3.1, we then have

The second string of equalities, appearing in the middle part of the proof Proposition 3.5, might appear counter-intuitive at first glance in that the left-most quantity is a norm squared of quantities which scale like $k$, whilst the right-most quantity scales like $k$ itself. However, the space $H_{k}$ itself depends on the kernel $k$, and scales inversely proportional to $k$, explaining that the identity is indeed dimensionally correct.

(Exponential convergence for the Gaussian kernel) The RKHS corresponding to the Gaussian kernel (3.3) is no longer isomorphic to a Sobolev space; it is contained in $H^{\tau}(X)$, for any $\tau<\infty$. For functions $f$ in this RKHS, Gaussian process regression with the Gaussian kernel converges exponentially in the fill distance $h_{U}$. For more details, see .

(Regression with non-zero mean) If in (3.1) we use a non-zero mean $m(\cdot)$, the formula for the predictive mean $m_{N}^{f}$ changes to

Approximation of the Bayesian posterior distribution

In this section, we analyse the error introduced in the posterior distribution $\mu^{y}$ when we use a Gaussian process emulator to approximate the parameter-to-observation map $\mathcal{G}$ or the negative log-likelihood $\Phi$. The aim is to show convergence, in a suitable sense, of the approximate posterior distributions to the true posterior distribution as the number of observations $N$ tends to infinity. For a given approximation $\mu^{y,N}$ of the posterior distribution $\mu^{y}$, we will focus on bounding the Hellinger distance between the two distributions, which is defined as

As proven in [15, Lemma 6.12 and 6.14], the Hellinger distance provides a bound for the Total Variation distance

and for $f\in L^{2}_{\mu^{y}}(X)\cap L^{2}_{\mu^{y,N}}(X)$, the Hellinger distance also provides a bound on the error in expected values

Suppose $\sup_{u\in X}\|\mathcal{G}(u)-m^{\mathcal{G}}_{N}(u)\|$ and $\sup_{u\in X}|\Phi(u)-m^{\Phi}_{N}(u)|$ converge to 0 as $N$ tends to $\infty$, and assume $\sup_{u\in X}\|\mathcal{G}(u)\|\leq C_{\mathcal{G}}$. Then there exist positive constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

where $C_{1}$ is independent of $U$ and $N$.

Since $\sup_{u\in X}|\Phi(u)|$ is bounded when $\sup_{u\in X}\|\mathcal{G}(u)\|$ is bounded, the fact that every convergent sequence is bounded again gives

We would like to point out here that the assumptions in Lemma 4.1 can be relaxed to assuming that the sequences $\sup_{u\in X}\|\mathcal{G}(u)-m^{\mathcal{G}}_{N}(u)\|$ and $\sup_{u\in X}|\Phi(u)-m^{\Phi}_{N}(u)|$ are bounded, since this is sufficient to prove the result.

Under the Assumptions of Lemma 4.1, there exist constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

For the first term, we use the local Lipschitz continuity of the exponential function, together with the equality $a^{2}-b^{2}=(a-b)(a+b)$ and the reverse triangle inequality to bound

As in equation (4.1), the first supremum can be bounded independently of $U$ and $N$, from which it follows that

for a constant $C$ independent of $U$ and $N$. For the second term, a very similar argument, together with Lemma 4.1 and Jensen’s inequality, shows

for a constant $C$ independent of $U$ and $N$.

Using Lemma 4.1 and Jensen’s inequality, we furthermore have

for a constant $C$ independent of $U$ and $N$. ∎

We remark here that Theorem 4.2 does not make any assumptions on the predictive means $m_{N}^{\mathcal{G}}$ and $m_{N}^{\Phi}$ other than the requirement that $\sup_{u\in X}\|\mathcal{G}(u)-m^{\mathcal{G}}_{N}(u)\|$ and $\sup_{u\in X}|\Phi(u)-m^{\Phi}_{N}(u)|$ converge to 0 as $N$ tends to $\infty$. Whether the predictive means are defined as in (3.6), or are derived by alternative approaches to Gaussian process regression , does not affect the conclusions of Theorem 4.2. Under Assumption 2.2, we can combine Theorem 4.2 with Proposition 3.3 (or Proposition 3.4) with $\beta=0$ to obtain error bounds in terms of the fill distance of the design points.

Suppose $m_{N}^{\Phi}$ and $m_{N}^{\mathcal{G}^{j}}$, $j=1,\dots,J$, are defined as in (3.6), with Matèrn kernel $k=k_{\nu,\lambda,\sigma_{k}^{2}}$. Suppose Assumption 2.2 holds with $s=\nu+K/2$, and the assumptions of Proposition 3.3 and Theorem 4.2 are satisfied. Then there exist constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

If Assumption 2.2 holds only for some $s<\nu+K/2$, an analogue of Corollary 4.3 can be proved using Proposition 3.4 with $\beta=0$. As already discussed in section 3.2, translating convergence rates in terms of the fill distance $h_{U}$ into rates in terms of the number of points $N$ typically leads to a strong dependence on the input dimension $K$. For uniform tensor grids $U$, the rates of convergence in $N$ predicted by Corollary 4.3 are given in Table 1.

2 Approximations based on the predictive process

Alternative to the mean-based approximations considered in the previous section, we now consider approximations to the posterior distribution $\mu^{y}$ obtained using the full predictive processes $\mathcal{G}_{N}$ and $\Phi_{N}$. In contrast to the mean, the full Gaussian processes also carry information about the uncertainty in the emulator due to only using a finite number of function evaluations to construct it.

Deterministic approximations of the posterior distribution $\mu^{y}$ can now be obtained by taking the expected value with respect to the predictive processes $\mathcal{G}_{N}$ and $\Phi_{N}$. This results in the marginal approximations

Firstly, we recall the following classical results from the theory of Gaussian measures on Banach spaces .

(Fernique’s Theorem) Let $E$ be a separable Banach space and $\nu$ a centred Gaussian measure on $(E,\mathcal{B}(E))$. If $\lambda,r>0$ are such that

Recall that, as in (3.1), $\Phi_{0}$ and $\mathcal{G}^{j}_{0}$ denote the initial Gaussian process models for $\Phi$ and $\mathcal{G}^{j}$, respectively, and, as in (3.5), $\Phi_{N}$ and $\mathcal{G}^{j}_{N}$ denote the conditioned Gaussian process models for $\Phi$ and $\mathcal{G}^{j}$, respectively.

From Jensen’s inequality, it then follows that

To determine $C_{1}$, we use the triangle inequality to bound, for any $1\leq p<\infty$,

The first factor can be bounded independently of $U$ and $N$ using the triangle inequality, together with $\sup_{u\in X}\|\mathcal{G}(u)\|\leq C_{\mathcal{G}}$ and $\sup_{u\in X}\left\|\mathcal{G}(u)-m^{\mathcal{G}}_{N}(u)\right\|\rightarrow 0$ as $N\rightarrow\infty$. For the second factor, we use Fernique’s Theorem (Proposition 4.4). First, we note that (using independence)

We would like to point out here that the assumption that $\sup_{u\in X}k_{N}(u,u)$ converges to 0 as $N$ tends to infinity in Lemma 4.7 is crucial in order to enable the choice of any $1\leq p<\infty$. This is related to the fact that the parameter $\lambda$ needs to be sufficiently small compared to $\sup_{u\in X}k_{N}(u,u)$ in order to satisfy the assumptions of Fernique’s Theorem.

In Lemma 4.7, we supposed that the assumptions of the Sudakov-Fernique inequality hold, for $g=\Phi_{0}$ and $f=\Phi_{N}-m_{N}^{\Phi}$, and for $g=\mathcal{G}^{j}_{0}$ and $f=\mathcal{G}^{j}_{N}-m_{N}^{\mathcal{G}^{j}}$, for $j\in\{1,\dots,J\}$. This is an assumption on the predictive variance $k_{N}$. In the following Lemma, we prove this assumption for the predictive variance given in (3.6).

Suppose the predictive variance $k_{N}$ is given by (3.6). Then the assumptions of the Sudakov-Fernique inequality hold, for $g=\Phi_{0}$ and $f=\Phi_{N}-m_{N}^{\Phi}$, and for $g=\mathcal{G}^{j}_{0}$ and $f=\mathcal{G}^{j}_{N}-m_{N}^{\mathcal{G}^{j}}$, for $j\in\{1,\dots,J\}$.

since the matrix $K(U,U)^{-1}$ is positive definite. ∎

We are now ready to prove bounds on the approximation error in the posterior distributions.

Under the assumptions of Lemma 4.7, there exist constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

For the first term, we use the (in)equalities $a-b=(a^{2}-b^{2})/(a+b)$ and $(\sqrt{a}+\sqrt{b})^{2}\geq a+b$, for $a,b>0$, to derive

For the first factor, using the convexity of $1/x$ on $(0,\infty)$, together with Jensen’s inequality, we have for all $u\in X$ the bound

As in the proof of Lemma 4.7, it then follows by Fernique’s Theorem that the right hand side can be bounded by a constant independent of $U$ and $N$.

For the second factor in the bound on $\frac{Z}{2}I$, the linearity of expectation, the local Lipschitz continuity of the exponential function, the equality $a^{2}-b^{2}=(a-b)(a+b)$, the reverse triangle inequality and Hölder’s inequality with conjugate exponents $p=(1+\delta)/\delta$ and $q=1+\delta$ give

for any $\delta>0$. The supremum in the above expression can be bounded by a constant independent of $U$ and $N$ by Fernique’s Theorem as in the proof of Lemma 4.7, since $\sup_{u\in X}\|\mathcal{G}(u)\|\leq C_{\mathcal{G}}<\infty$. It follows that there exists a constant $C$ independent of $U$ and $N$ such that

For the second term in the bound on the Hellinger distance, we have

Using the linearity of expectation, Tonelli’s Theorem and Jensen’s inequality, we have

which can now be bounded as before. The first claim of the theorem now follows by Lemma 4.7.

The first factor can again be bounded using Jensen’s inequality,

which as in the proof of Lemma 4.7, can be bounded by a constant independent of $U$ and $N$ by Fernique’s Theorem. For the second factor in the bound on $\frac{Z}{2}I$, the linearity of expectation and the local Lipschitz continuity of the exponential function give

For the second term in the bound on the Hellinger distance, the linearity of expectation, Tonelli’s Theorem and Jensen’s inequality give

which can now be bounded as before. The second claim of the theorem then follows by Lemma 4.7. ∎

Similar to Theorem 4.2, Theorem 4.9 provides error bounds for general Gaussian process emulators of $\mathcal{G}$ and $\Phi$. An example of a Gaussian process emulator that satisfies the assumptions of Theorem 4.9 is the emulator defined by (3.6), however, other choices are possible. As in Corollary 4.3, we can now combine Assumption 2.2, Theorem 4.9 and Proposition 3.3 with $\beta=0$ to derive error bounds in terms of the fill distance.

Suppose $\mathcal{G}_{N}$ and $\Phi_{N}$ are defined as in (3.6), with Matèrn kernel $k=k_{\nu,\lambda,\sigma_{k}^{2}}$. Suppose Assumption 2.2 holds with $s=\nu+K/2$, and the assumptions of Proposition 3.3 and Theorem 4.9 are satisfied. Then there exist constants $C_{1},C_{2},C_{3}$ and $C_{4}$, independent of $U$ and $N$, such that

The first term can be bounded by using Assumption 2.2, Proposition 3.2 and Proposition 3.3,

for a constant $C$ independent of $U$ and $N$. The second term can be bounded by using Assumption 2.2, Proposition 3.2, Proposition 3.3, Proposition 3.5, the linearity of expectation and the Sobolev Embedding Theorem

for a constant $C$ independent of $U$ and $N$. The claim of the corollary then follows. ∎

If Assumption 2.2 holds only for some $s<\nu+K/2$, an analogue of Corollary 4.10 can be proved using Proposition 3.4 with $\beta=0$.

Under the Assumptions of Lemma 4.7, there exist constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

For the first term, Tonelli’s Theorem, the local Lipschitz continuity of the exponential function, the equality $a^{2}-b^{2}=(a-b)(a+b)$, the reverse triangle inequality and Hölder’s inequality with conjugate exponents $p=(1+\delta)/\delta$ and $q=1+\delta$ give

for any $\delta>0$. The supremum in the above bound can be bounded independently of $U$ and $N$ by Fernique’s Theorem as in the proof of Lemma 4.7. It follows that there exists a constant $C$ independent of $U$ and $N$ such that

For the second term in the bound on the Hellinger distance, we have

By Jensen’s inequality and the same argument as above, we have

Together with Tonelli’s Theorem and Hölder’s inequality with conjugate exponents $p=(1+\delta)/\delta$ and $q=1+\delta$, we then have

for any $\delta>0$. The supremum in the bound above can be bounded independently of $U$ and $N$ by Lemma 4.7 and Fernique’s Theorem. The first claim of the Theorem then follows.

Together with Tonelli’s Theorem and Hölder’s inequality with conjugate exponents $p=(1+\delta)/\delta$ and $q=1+\delta$, we then have

for any $\delta>0$. The first expected value in the bound above can be bounded independently of $U$ and $N$ by Lemma 4.7. The second claim of the Theorem then follows. ∎

Similar to Theorem 4.2 and Theorem 4.9, Theorem 4.11 provides error bounds for general Gaussian process emulators of $\mathcal{G}$ and $\Phi$. As a particular example, we can take the emulators defined by (3.6). We can now combine Assumption 2.2, Theorem 4.11 and Proposition 3.3 with $\beta=0$ to derive error bounds in terms of the fill distance.

Suppose $\mathcal{G}_{N}$ and $\Phi_{N}$ are defined as in (3.6), with Matèrn kernel $k=k_{\nu,\lambda,\sigma_{k}^{2}}$. Suppose Assumption 2.2 holds with $s=\nu+K/2$, and the assumptions of Proposition 3.3 and Theorem 4.11 are satisfied. Then there exist constants $C_{1},C_{2},C_{3}$ and $C_{4}$, independent of $U$ and $N$, such that

If Assumption 2.2 holds only for some $s<\nu+K/2$, an analogue of Corollary 4.12 can be proved using Proposition 3.4 with $\beta=0$.

We furthermore have the following result on a generalised total variation distance , defined by

Under the Assumptions of Lemma 4.7, there exist constants $C_{1}$ and $C_{2}$, independent of $U$ and $N$, such that

Numerical Examples

We consider the model inverse problem of determining the diffusion coefficient of an elliptic partial differential equation (PDE) in divergence form from observation of a finite set of noisy continuous functionals of the solution. This type of equation arises, for example, in the modelling of groundwater flow in a porous medium. We consider the one-dimensional model problem

where the coefficient $\kappa$ depends on parameters $u=\{u_{j}\}_{j=1}^{K}\in^{K}$ through the linear expansion

In this setting the forward map $G:^{K}\rightarrow H^{1}_{0}(D)$, defined by $G(u)=p$, is an analytic function . Since the observation operator $\mathcal{O}$ is linear and bounded, Assumption 2.2 is satisfied for any $s>K/2$.

Unless stated otherwise, we will throughout this section approximate the solution $p$ by standard, piecewise linear, continuous finite elements on a uniform grid with mesh size $h=1/32$. The corresponding approximate forward map, denoted by $G_{h}$, is also an analytic function of $u$ , and Assumption 2.2 is satisfied for any $s>K/2$ also for $G_{h}$. By slight abuse of notation, we will denote the posterior measure corresponding to the forward map $G_{h}$ by $\mu^{y}$, and use this as our reference measure. The error induced by the finite element approximation will be ignored.

The emulators $\mathcal{G}_{N}$ and $\Phi_{N}$ are computed as described in section 3.2, with mean and covariance kernel given by (3.6). In the Gaussian process prior (3.1), we choose $m\equiv 0$ and $k=k_{\nu,1,1}$, a Matèrn kernel with variance $\sigma_{k}^{2}=1$, correlation length $\lambda=1$ and smoothness parameter $\nu$.

For a given approximation $\mu^{y,N}$ to $\mu^{y}$, we will compute twice the Hellinger distance squared,

The integral over $^{K}$ is approximated by a randomly shifted lattice rule with product weight parameters $\gamma_{j}=1/j^{2}$ . The generating vector for the rule used is available from Frances Kuo’s website (http://web.maths.unsw.edu.au/$\sim$fkuo/) as “lattice-39102-1024-1048576.3600”. For the marginal and random approximations, the expected value over the Gaussian process is approximated by Monte Carlo sampling, using the MATLAB command mvnrnd.

2 Marginal approximations

3 Random approximations

Conclusions and further work

Gaussian process emulators are frequently used as surrogate models. In this work, we analysed the error that is introduced in the Bayesian posterior distribution when a Gaussian process emulator is used to approximate the forward model, either in terms of the parameter-to-observation map or the negative log-likelihood. We showed that the error in the posterior distribution, measured in the Hellinger distance, can be bounded in terms of the error in the emulator, measured in a norm dependent on the approximation considered.

An issue that requires further consideration is the efficient emulation of vector-valued functions. A simple solution, employed in this work, is to emulate each entry independently. In many applications, however, it is natural to assume that the entries are correlated, and a better emulator could be constructed by including this correlation in the emulator. Furthermore, there are still a lot of open questions about how to do this optimally . Also the question of scaling the Gaussian process methodology to high dimensional input spaces remains open. The current error bounds from scattered data approximation employed in this paper feature a strong dependence on the input dimension $K$, yielding poor convergence estimates in high dimensions.

Another important issue is the selection of the design points used to construct the Gaussian process emulator, also known as experimental design. In applications where the posterior distribution concentrates with respect to the prior, it might be more efficient to choose design points that are somehow adapted to the posterior measure instead of space-filling designs that have a small fill distance. For example, we could use the sequential designs in . It would be interesting to prove suitable error bounds in this case, maybe using ideas from .

In practical applications of Gaussian process emulators, such as in , the derivation of the emulator is often more involved than the simple approach presented in section 3. The hyper-parameters in the covariance kernel of the emulator are often unknown, and there is often a discrepancy between the mathematical model of the forward map and the true physical process, known as model error. These are both important issues for which the assumptions in our error bounds have not yet been verified.

References