On the Convergence of the Laplace Approximation and Noise-Level-Robustness of Laplace-based Monte Carlo Methods for Bayesian Inverse Problems

TL;DR

This paper investigates the convergence and noise-level robustness of Laplace approximation-based Monte Carlo methods for Bayesian inverse problems, demonstrating their effectiveness in handling concentrated posteriors with theoretical and numerical support.

Contribution

It introduces the use of Laplace approximation as a base measure for Monte Carlo integration, providing convergence analysis and robustness results for Bayesian inverse problems.

Findings

Laplace-based importance sampling is robust to posterior concentration.

Laplace-based quasi-Monte Carlo methods outperform prior-based methods.

Numerical experiments confirm theoretical robustness and convergence.

Abstract

The Bayesian approach to inverse problems provides a rigorous framework for the incorporation and quantification of uncertainties in measurements, parameters and models. We are interested in designing numerical methods which are robust w.r.t. the size of the observational noise, i.e., methods which behave well in case of concentrated posterior measures. The concentration of the posterior is a highly desirable situation in practice, since it relates to informative or large data. However, it can pose a computational challenge for numerical methods based on the prior or reference measure. We propose to employ the Laplace approximation of the posterior as the base measure for numerical integration in this context. The Laplace approximation is a Gaussian measure centered at the maximum a-posteriori estimate and with covariance matrix depending on the logposterior density. We discuss convergence results of the Laplace approximation in terms of the Hellinger distance and analyze the efficiency of Monte Carlo methods based on it. In particular, we show that Laplace-based importance sampling and Laplace-based quasi-Monte-Carlo methods are robust w.r.t. the concentration of the posterior for large classes of posterior distributions and integrands whereas prior-based importance sampling and plain quasi-Monte Carlo are not. Numerical experiments are presented to illustrate the theoretical findings.

Tables

Table 1. Table 1: Frequently used notation.

Symbol	Meaning	Reference
${\mu }_0$	prior probability measure	(5)
${\pi }_0$	Lebesgue density of the prior: ${\mu }_0\left(dx\right)={\pi }_0\left(x\right)dx$	(5)
$Z_0=1$	hypothetically: Normalization constant of ${\pi }_0$
${\mu }_n$	posterior probability measure	(1), (6)
${\Phi }_n$	(scaled) negative loglikelihood of posterior	(1)
$I_n$	(scaled) negative logdensity of posterior	(7)
${\pi }_n\left(x\right)$	unnormalized Lebesgue density of the posterior	Lemma 1
$Z_n$	normalization constant of ${\pi }_n$	(1)
${ℒ}_{{\mu }_n}$	Laplace approximation of ${\mu }_n$	(9)
${\tilde{I}}_n=T_2{I}_n\left(x_n\right)$	(scaled) negative logdensity of ${ℒ}_{{\mu }_n}$	(10)
${\tilde{\pi }}_n$	unnormalized Lebesgue density of ${ℒ}_{{\mu }_n}$	Lemma 1
${\tilde{Z}}_n$	normalization constant of ${\tilde{\pi }}_n$	(9)

Equations

$${\mu }_n(\mathrm{d}x)=\frac{1}{Z_n}\exp \left(-n{\Phi }_n(x)\right){\mu }_0(\mathrm{d}x),Z_n:={\int }_{{\mathbb{R}}^d}\exp \left(-n{\Phi }_n(x)\right){\mu }_0(\mathrm{d}x),$$

$${\int }_{{\mathbb{R}}^d}f(x)\exp (-n\Phi (x))\mathrm{d}x=f(x_{\star })\frac{(2\pi {)}^{d/2}\exp (-n\Phi (x_{\star }))}{n^{d/2}\sqrt{\det \left({\nabla }^2\Phi (x_{\star })\right)}}\left(1+\mathcal{O}(n^{-1})\right)$$

$$d_{\mathrm{H}}({\mu }_n,{\mathcal{L}}_{{\mu }_n})\in \mathcal{O}(n^{-1/2}).$$

$$n^{\delta }{e}_{n,N}(f)\underset{n\to \infty }{\overset{\mathbb{P}}{\to }}0,\delta <1/2.$$

$$J(n):={\int }_{D}f(x)\exp (-n\Phi (x))\mathrm{d}x,n\in \mathbb{N},$$

$${\delta }_r:=\underset{x\in B_r^c(x_{\star })}{\inf }\Phi (x)-\Phi (x_{\star })>0,$$

$$J(n)=\operatorname{\mathrm{e}}^{-n\Phi(x_{\star})}\,n^{-d/2}\,\left(\sum_{k=0}^{p}c_{k}(f)n^{-{\color[rgb]{0,0,0}k}}+\mathcal{O}\left(n^{-p-1}\right)\right)$$

$$\underset{n\to \infty }{\lim }\frac{J(n)}{c_0(f)\exp (-n\Phi (x_{\star }))n^{-d/2}}=1$$

$$\left|{\int }_{{\mathbb{R}}^d}f(x)\exp (-n\Phi (x))\mathrm{d}x-{\int }_{{\mathbb{R}}^d}f(x)\exp (-\frac{n}{2}\Vert x-x_{\star }{\Vert }_{H_{\star }}^2)\mathrm{d}x\right|\in o(n^{-d/2})$$

$$\underset{n\to \infty }{\lim }\frac{{\int }_{{\mathbb{R}}^d}f(x)\exp (-n\Phi (x))\mathrm{d}x}{{\int }_{{\mathbb{R}}^d}f(x)\exp (-\frac{n}{2}\Vert x-x_{\star }{\Vert }_{H_{\star }}^2)\mathrm{d}x}=1\text{ as }n\to \infty ,$$

$$\underset{n\to \infty }{\lim }{\int }_{{\mathbb{R}}^d}f(x){\mu }_n(\mathrm{d}x)$$

$$=\underset{n\to \infty }{\lim }\frac{c_0(f{\pi }_0)n^{-d/2}}{c_0({\pi }_0)n^{-d/2}}=\frac{c_0(f{\pi }_0)}{c_0({\pi }_0)}=f(x_{\star }).$$

$$\underset{n\to \infty }{\lim }\left|{\int }_{{\mathbb{R}}^d}f(x){\mu }_n(\mathrm{d}x)-{\int }_{{\mathbb{R}}^d}f(x){\mathcal{N}}_{x_{\star },(n{H}_{\star }{)}^{-1}}(\mathrm{d}x)\right|=0,$$

$$\left|{\int }_{{\mathbb{R}}^d}f(x){\mu }_n(\mathrm{d}x)-{\int }_{{\mathbb{R}}^d}f(x){\mathcal{N}}_{x_{\star },(n{H}_{\star }{)}^{-1}}(\mathrm{d}x)\right|\in \mathcal{O}(n^{-1}).$$

$$y=\mathcal{G}(x)+\eta ,$$

$$\Phi (x;y)=\frac{1}{2}|y-\mathcal{G}(x){|}_{\Gamma }^2.$$

$$\mu (\mathrm{d}x)=\frac{1}{Z}\exp (-\Phi (x;y)){\mu }_0(\mathrm{d}x)$$

$${\mu }_0(\mathrm{d}x)={\pi }_0(x)\mathrm{d}x,\text{ and we set }{\mathrm{S}}_0:=\{x\in {\mathbb{R}}^d:{\pi }_0(x)>0\}=\mathrm{supp}{\mu }_0.$$

$${\mu }_n(\mathrm{d}x)=\frac{1}{Z_n}{1}_{{\mathrm{S}}_0}(x)\exp \left(-n{I}_n(x)\right)\mathrm{d}x$$

$$I_n(x):={\Phi }_n(x)-\frac{1}{n}\log {\pi }_0(x).$$

$$I_n(x_n)=0,\nabla I_n(x_n)=0,{\nabla }^2{I}_n(x_n)>0,$$

$${\hat{\Phi }}_n(x):={\Phi }_n(x)-c,{\hat{I}}_n(x):={\hat{\Phi }}_n(x)-\log {\pi }_0(x)$$

$${\mu }_n(\mathrm{d}x)=\frac{1}{{\hat{Z}}_n}\exp \left(-n{\hat{\Phi }}_n(x)\right){\mu }_0(\mathrm{d}x),{\hat{Z}}_n={\int }_{{\mathbb{R}}^d}\exp \left(-n{\hat{\Phi }}_n(x)\right){\mu }_0(\mathrm{d}x),$$

$${\mathcal{L}}_{{\mu }_n}:=\mathcal{N}(x_n,n^{-1}{C}_n),C_n^{-1}:={\nabla }^2{I}_n(x_n).$$

$${\mathcal{L}}_{{\mu }_n}(\mathrm{d}x)=\frac{1}{{\tilde{Z}}_n}\exp \left(-\frac{n}{2}\Vert x-x_n{\Vert }_{C_n^{-1}}^2\right)\mathrm{d}x,{\tilde{Z}}_n:=n^{-d/2}\sqrt{\det (2\pi C_n)},$$

$${\tilde{I}}_n(x):=\frac{1}{2}\Vert x-x_n{\Vert }_{C_n^{-1}}^2=\underset{=0}{\underbrace{I_n(x_n)}}+\underset{=0}{\underbrace{\nabla I_n(x_n{)}^{\top }}}(x-x_n)+\frac{1}{2}\Vert x-x_n{\Vert }_{{\nabla }^2{I}_n(x_n)}^2$$

$${\mu }_n\approx {\mathcal{L}}_{{\mu }_n},\frac{1}{Z_n}\exp (-n{I}_n(x))\approx \frac{1}{{\tilde{Z}}_n}\exp \left(-n{\tilde{I}}_n(x)\right).$$

$${\mu }_0(\mathrm{d}x)$$

$${\mu }_n(\mathrm{d}x)$$

$$=Z_n^{-1}\exp (-n{I}_n(x))1_{S_0}(x)\mathrm{d}x$$

$${\mathcal{L}}_{{\mu }_n}(\mathrm{d}x)$$

$$={\tilde{Z}}_n^{-1}\exp (-\frac{n}{2}\Vert x-x_n{\Vert }_{{\nabla }^2{I}_n(x_n)}^2)\mathrm{d}x$$

$$\left|{\int }_{{\mathbb{R}}^d}f(x){\mu }_n(\mathrm{d}x)-{\int }_{{\mathbb{R}}^d}f(x){\mathcal{L}}_{{\mu }_n}(\mathrm{d}x)\right|\in \mathcal{O}(n^{-1}).$$

Full text

Keywords: Bayesian inverse problems, Laplace approximation, importance sampling, quasi-Monte Carlo, uncertainty quantification

Mathematics Subject Classification: 65M32, 62F15, 60B10, 65C05, 65K10

1 Introduction

The identification of unknown parameters from noisy observations arises in various areas of application, e.g., engineering systems, biological models, environmental systems. In recent years, Bayesian inference has become a popular approach to model inverse problems [39], i.e., noisy observations are used to update the knowledge of unknown parameters from a prior distribution to the posterior distribution. The latter is then the solution of the Bayesian inverse problem and obtained by conditioning the prior distribution on the data. This approach is very appealing in various fields of applications, since uncertainty quantification can be performed, once the prior distribution is updated—barring the fact that Bayesian credible sets are not in a one-to-one correspondence to classical confidence sets, see [7, 40].

To ensure the applicability of the Bayesian approach to computationally demanding models, there has been a lot of research effort towards improved algorithms allowing for effective sampling or integration w.r.t. the resulting posterior measure. For example, the computational burden of expensive forward or likelihood models can be reduced by surrogates or multilevel strategies [27, 14, 20, 35] and for many classical sampling or integration methods such as Quasi-Monte Carlo [12], Markov chain Monte Carlo [6, 33, 42], and numerical quadrature [36, 5] we now know modifications and conditions which ensure a dimension-independent efficiency.

However, a completely different, but very common challenge for many numerical methods has drawn surprisingly less attention so far: the challenge of concentrated posterior measures such as

[TABLE]

Here, $n\gg 1$ and $\mu_{0}$ denotes a reference or prior probability measure on $\mathbb{R}^{d}$ and $\Phi_{n}\colon\mathbb{R}^{d}\to[0,\infty)$ are negative log-likelihood functions resulting, e.g., from $n$ observations.

From a modeling point of view the concentration effect of the posterior is a highly desirable situation due to large data sets and less remaining uncertainty about the parameter to be inferred. From a numerical point of view, on the other hand, this can pose a delicate situation, since standard integration methods may perform worse and worse if the concentration increases due to $n\to\infty$. Hence, understanding how sampling or quadrature methods for $\mu_{n}$ behave as $n\to\infty$ is a crucial task with immediate benefits for purposes of uncertainty quantification. Since small noise yields “small” uncertainty, one might be tempted to consider only optimization-based approaches in order to compute a point estimator (i.e., the maximum a-posteriori estimator) for the unknown parameter which is usually computationally much cheaper than a complete Bayesian inference. However, for quantifying the remaining risk, e.g., computing the posterior failure probability for some quantity of interest, we still require efficient integration methods for concentrated posteriors as $\mu_{n}$. Nonetheless, we will use well-known preconditioning techniques from numerical optimization in order to derive such robust integration methods for the small noise setting.

Numerical methods are often based on the prior $\mu_{0}$, since $\mu_{0}$ is usually a simple measure allowing for direct sampling or explicit quadrature formulas. However, for large $n$ most of the corresponding sample points or quadrature nodes will be placed in regions of low posterior importance missing the needle in the haystack—the minimizers of $\Phi_{n}$. An obvious way to circumvent this is to use a numerical integration w.r.t. another reference measure which can be straightforwardly computed or sampled from and concentrates around those minimizers and shrinks like the posterior measures $\mu_{n}$ as $n\to\infty$. In this paper we consider numerical methods based on a Gaussian approximation of $\mu_{n}$—the Laplace approximation.

When it comes to integration w.r.t. an increasingly concentrated function, the well-known and widely used Laplace’s method provides explicit asymptotics for such integrals, i.e., under certain regularity conditions [44] we have for $n\to\infty$ that

[TABLE]

where $x_{\star}\in\mathbb{R}$ denotes the assumed unique minimizer of $\Phi\colon\mathbb{R}^{d}\to\mathbb{R}$. This formula is derived by approximating $\Phi$ by its second-order Taylor polynomial at $x_{\star}$. We could now use (2) and its application to $Z_{n}$ in order to derive that $\int_{\mathbb{R}^{d}}f(x)\ \mu_{n}(\mathrm{d}x)\to f(x_{\star})$ as $n\to\infty$. However, for finite $n$ this is only of limited use, e.g., consider the computation of posterior probabilities where $f$ is an indicator function. Thus, in practice we still rely on numerical integration methods in order to obtain a reasonable approximation of the posterior integrals $\int_{\mathbb{R}^{d}}f(x)\ \mu_{n}(\mathrm{d}x)$. Nonetheless, the second-order Taylor approximation employed in Laplace’s method provides us with (a guideline to derive) a Gaussian measure approximating $\mu_{n}$.

This measure itself is often called the Laplace approximation of $\mu_{n}$ and will be denoted by $\mathcal{L}_{\mu_{n}}$. Its mean is given by the maximum a-posteriori estimate (MAP) of the posterior $\mu_{n}$ and its covariance is the inverse Hessian of the negative log posterior density. Both quantities can be computed efficiently by numerical optimization and since it is a Gaussian measure it allows for direct samplings and easy quadrature formulas. The Laplace approximation is widely used in optimal (Bayesian) experimental design to approximate the posterior distribution (see, for example, [1]) and has been demonstrated to be particularly useful in the large data setting, see [25, 34] and the references therein for more details. Moreover, in several recent publications the Laplace approximation was already proposed as a suitable reference measure for numerical quadrature [38, 5] or importance sampling [2]. Note that preconditioning strategies based on Laplace approximation are also referred to as Hessian-based strategies due to the equivalence of the inverse covariance and the Hessian of the corresponding optimization problem, cp. [5]. In [38], the authors showed that a Laplace approximation-based adaptive Smolyak quadrature for Bayesian inference with affine parametric operator equations exhibits a convergence rate independent of the size of the noise, i.e., independent of $n$.

This paper extends the analysis in [38] for quadrature to the widely applied Laplace-based importance sampling and Laplace-based quasi-Monte Carlo (QMC) integration.

Before we investigate the scale invariance or robustness of these methods we examine the behaviour of the Laplace approximation and in particular, the density $\frac{\mathrm{d}\mu_{n}}{\mathrm{d}\mathcal{L}_{\mu_{n}}}$. The reason behind is that, for importance sampling as well as QMC integration, this density naturally appears in the methods, hence, if it deteriorates as $n\to\infty$, this will be reflected in a deteriorating efficiency of the method. For example, for $\Phi_{n}\equiv\Phi$ the density w.r.t. the prior measure $\frac{\mathrm{d}\mu_{n}}{\mathrm{d}\mu_{0}}=\exp(-n\Phi)/Z_{n}$ deteriorates to a Dirac function at the minimizer $x_{\star}$ of $\Phi$ as $n\to\infty$ which causes the shortcomings of Monte Carlo or QMC integration w.r.t. $\mu_{0}$ as $n\to\infty$. However, for the Laplace approximation we show that the density $\frac{\mathrm{d}\mu_{n}}{\mathrm{d}\mathcal{L}_{\mu_{n}}}$ converges $\mathcal{L}_{\mu_{n}}$-almost everywhere to $1$ which in turn results in a robust—and actually improving—performance w.r.t. $n$ of related numerical methods. In summary, the main results of this paper are the following:

Laplace Approximation: Given mild conditions the Laplace approximation $\mathcal{L}_{\mu_{n}}$ converges in Hellinger distance to $\mu_{n}$:

[TABLE]

This result is closely related to the well-known Bernstein–von Mises theorem for the posterior consistency in Bayesian inference [41]. The significant difference here is that the covariance in the Laplace approximation depends on the data and the convergence holds for the particularly observed data whereas in the classical Bernstein–von Mises theorem the covariance is the inverse of the expected Fisher information matrix and the convergence is usually stated in probability. 2. 2.

Importance Sampling: We consider integration w.r.t. measures $\mu_{n}$ as in (1) where $\Phi_{n}(x)=\Phi(x)-\iota_{n}$ for a $\Phi\colon\mathbb{R}^{d}\to[0,\infty)$ and $\iota_{n}\in\mathbb{R}$.

• •

Prior-based Importance Sampling: We consider the case of prior-based importance sampling, i.e., the prior $\mu_{0}$ is used as the importance distribution for computing the expectation of smooth integrands $f\in L^{2}_{\mu_{0}}(\mathbb{R})$. Here, the asymptotic variance w.r.t. such measures $\mu_{n}$ deteriorates like $n^{d/2-1}$.

• •

Laplace-based Importance Sampling. The (random) error $e_{n,N}(f)$ of Laplace-based importance sampling for computing expectations of smooth integrands $f\in L^{2}_{\mu_{0}}(\mathbb{R})$ w.r.t. such measures $\mu_{n}$ using a fixed number of samples $N\in\mathbb{N}$ decays in probability almost like $n^{-1/2}$, i.e.,

[TABLE] 3. 3.

Quasi-Monte Carlo: We focus for the analysis of the quasi-Monte Carlo methods on the bounded case of $\mu_{0}=\mathcal{U}([\frac{1}{2},\frac{1}{2}]^{d})$.

• •

Prior-based Quasi-Monte Carlo: The root mean squared error estimate for computing integrals of the form (2) by QMC using randomly shifted Lattice rules deteriorates like $n^{d/4}$ as $n\to\infty$.

• •

Laplace-based Quasi-Monte Carlo: If the lattice rule is transformed by an affine mapping based on the mean and the covariance of the Laplace approximation, then the resulting root mean squared error decays like $n^{-d/2}$ for integrals of the form (2).

The outline of the paper is as follows: in Section 2 we introduce the Laplace approximation for measures of the form (1) and the notation of the paper. In Section 2.2 we study the convergence of the Laplace approximation. We also consider the case of singular Hessians or perturbed Hessians and provide some illustrative numerical examples. At the end of the section, we shortly discuss the relation to the classical Bernstein–von Mises theorem. The main results about importance sampling and QMC using the prior measure and the Laplace approximation, respectively, are then discussed in Section 3. We also briefly comment on existing results for numerical quadrature and provide several numerical examples illustrating our theoretical findings. The appendix collects the rather lengthy and technical proofs of the main results.

2 Convergence of the Laplace Approximation

We start with recalling the classical Laplace method for the asymptotics of integrals.

Theorem 1** (variant of [44, Section IX.5]).**

Set

[TABLE]

where $D\subseteq\mathbb{R}^{d}$ is a possibly unbounded domain and let the following assumptions hold:

The integral $J(n)$ converges absolutely for each $n\in\mathbb{N}$. 2. 2.

There exists an $x_{\star}$ in the interior of $D$ such that for every $r>0$ there holds

[TABLE]

where $B_{r}(x_{\star})\coloneqq\{x\in\mathbb{R}^{d}\colon\|x-x_{\star}\|\leq r\}$ and $B^{c}_{r}(x_{\star})\coloneqq\mathbb{R}^{d}\setminus B_{r}(x_{\star})$. 3. 3.

In a neighborhood of $x_{\star}$ the function $f:D\to\mathbb{R}$ is $(2p+2)$ times continuously differentiable and $\Phi\colon\mathbb{R}^{d}\to\mathbb{R}$ is $(2p+3)$ times continuously differentiable for a $p\geq 0$, i.e., and the Hessian $H_{\star}\coloneqq\nabla^{2}\Phi(x_{\star})$ is positive definite.

Then, as $n\to\infty$, we have

[TABLE]

where $c_{k}(f)\in\mathbb{R}$ and, particularly, c_{0}(f)={\color[rgb]{0,0,0}\sqrt{\det(2\pi H^{-1}_{\star})}\,f(x_{\star})}.

Remark 1*.*

As stated in [44, Section IX.5] the asymptotic

[TABLE]

with $c_{0}(f)$ is as above, already holds for $f\colon\mathbb{R}^{d}\to\mathbb{R}$ being continuous and $\Phi\colon\mathbb{R}^{d}\to\mathbb{R}$ being twice continuously differentiable in a neighborhood of $x_{\star}$ with positive definite $\nabla^{2}\Phi(x_{\star})$—given that the first two assumptions of Theorem 1 are also satified.

Assume that $\Phi(x_{\star})=0$, then the above theorem and remark imply

[TABLE]

for continuous and integrable $f\colon\mathbb{R}^{d}\to\mathbb{R}$, where $\|\cdot\|_{A}=\|A^{1/2}\cdot\|$ for a symmetric positive definite matrix $A\in\mathbb{R}^{d\times d}$. This is similar to the notion weak convergence (albeit with two different non-static measures). If we additionally claim that $f(x_{\star})>0,$ then also

[TABLE]

which is sort of a relative weak convergence. In other words, the asymptotic behaviour of $\int f\ \operatorname{\mathrm{e}}^{-n\Phi}\ \mathrm{d}x$, in particular, its convergence to zero, is the same as of the integral of $f$ w.r.t. an unnormalized Gaussian density with mean in $x_{\star}$ and covariance $(nH_{\star})^{-1}$.

If we consider now probability measures $\mu_{n}$ as in (1) but with $\Phi_{n}\equiv\Phi$ where $\Phi$ satisfies the assumptions of Theorem 1, and if we suppose that $\mu_{0}$ possesses a continuous Lebesgue density $\pi_{0}\colon\mathbb{R}\to[0,\infty)$ with $\pi_{0}(x_{\star})>0$, then Theorem 1 and Remark 1 will imply for continuous and integrable $f\colon\mathbb{R}^{d}\to\mathbb{R}$ that

[TABLE]

The same reasoning applies to the expectation of $f$ w.r.t. a Gaussian measure $\mathcal{N}(x_{\star},(nH_{\star})^{-1})$ with unnormalized density $\exp(-\frac{n}{2}\|x-x_{\star}\|^{2}_{H_{\star}})$. Thus, we obtain the weak convergence of $\mu_{n}$ to $\mathcal{N}(x_{\star},(nH_{\star})^{-1})$, i.e., for any continuous and bounded $f\colon\mathbb{R}^{d}\to\mathbb{R}$ we have

[TABLE]

where $\mathcal{N}_{x,C}$ is short for $\mathcal{N}(x,C)$. In fact, for twice continously differentiable $f\colon\mathbb{R}^{d}\to\mathbb{R}$ we get by means of Theorem 1 the rate

[TABLE]

Note that due to normalization we do not need to assume $\Phi(x_{\star})=0$ here. Hence, this weak convergence suggests to use $\mathcal{N}_{x_{\star},(nH_{\star})^{-1}}$ as a Gaussian approximation to $\mu_{n}$. In the next subsection we derive similar Gaussian approximation for the general case $\Phi_{n}\not\equiv\Phi$, whereas subsection 2.2 includes convergence results of the Laplace approximation in terms of the Hellinger distance.

Bayesian inference

We present some context for the form of equation (1) in the following. Integrals of the form (1) arise naturally in the Bayesian setting for inverse problems with large amount of observational data or informative data. Note that the mathematical results for the Laplace approximation given in section 2 are derived in a much more general setting and are not restricted to integrals w.r.t. the posterior in the Bayesian inverse framework. We refer to [21, 8] and the references therein for a detailed introduction to Bayesian inverse problems.

Consider a continuous forward response operator $\mathcal{G}\colon\mathbb{R}^{d}\to\mathbb{R}^{K}$ mapping the unknown parameters $x\in\mathbb{R}^{d}$ to the data space $\mathbb{R}^{K}$, where $K\in\mathbb{N}$ denotes the number of observations. We investigate the inverse problem of recovering unknown parameters $x\in\mathbb{R}^{d}$ from noisy observation $y\in\mathbb{R}^{K}$ given by

[TABLE]

where $\eta\sim\mathcal{N}(0,\Gamma)$ is a Gaussian random variable with mean zero and covariance matrix $\Gamma$, which models the noise in the observations and in the model.

The Bayesian approach for this inverse problem of inferring $x$ from $y$ (which is ill-posed without further assumptions) works as follows: For fixed $y\in\mathbb{R}^{K}$ we introduce the least-squares functional (or negative loglikelihood in the language of statistics) $\Phi(\cdot;y):\mathbb{R}^{d}\to\mathbb{R}$ by

[TABLE]

with $|\cdot|_{\Gamma}\coloneqq|\Gamma^{-\frac{1}{2}}\cdot|$ denoting the weighted Euclidean norm in $\mathbb{R}^{K}$. The unknown parameter $x$ is modeled as a $\mathbb{R}^{d}$-valued random variable with prior distribution $\mu_{0}$ (independent of the observational noise $\eta$), which regularizes the problem and makes it well-posed by application of Bayes’ theorem: The pair $(x,y)$ is a jointly varying random variable on $\mathbb{R}^{d}\times\mathbb{R}^{K}$ and hence the solution to the Bayesian inverse problem is the conditional or posterior distribution $\mu$ of $x$ given the data $y$ where the law $\mu$ is given by

[TABLE]

with the normalization constant $Z\coloneqq\int_{\mathbb{R}^{d}}\exp(-\Phi(x;y))\mu_{0}(\mathrm{d}x)$. If we assume a decaying noise-level by introducing a scaled noise covariance $\Gamma_{n}=\frac{1}{n}\Gamma$, the resulting noise model $\eta_{n}\sim N(0,\Gamma_{n})$ yields an $n$-dependent log-likelihood term which results in posterior measures $\mu_{n}$ of the form (1) with $\Phi_{n}(x)=\Phi(x;y)$. Similarly, an increasing number $n\in\mathbb{N}$ of data $y_{1},\ldots,y_{n}\in\mathbb{R}^{k}$ resulting from $n$ observations of $\mathcal{G}(x)$ with independent noises $\eta_{1},\ldots,\eta_{n}\sim N(0,\Gamma)$ yields posterior measures $\mu_{n}$ as in (1) with $\Phi_{n}(x)=\frac{1}{n}\sum_{j=1}^{n}\Phi(x;y_{j})$.

2.1 The Laplace Approximation

Throughout the paper, we assume that the prior measure $\mu_{0}$ is absolutely continuous w.r.t. Lebesgue measure with density $\pi_{0}\colon\mathbb{R}^{d}\to[0,\infty)$, i.e.,

[TABLE]

Hence, also the measures $\mu_{n}$ in (1) are absolutely continuous w.r.t. Lebesgue measure, i.e.,

[TABLE]

where $I_{n}\colon\mathrm{S}_{0}\to\mathbb{R}$ is given by

[TABLE]

In order to define the Laplace approximation of $\mu_{n}$ we need the following basic assumption.

Assumption 1**.**

There holds $\Phi_{n},\pi_{0}\in C^{2}(\mathrm{S}_{0},\mathbb{R})$, i.e., the mappings $\pi_{0},\Phi_{n}\colon\mathrm{S}_{0}\to\mathbb{R}$ are twice continuously differentiable. Furthermore, $I_{n}$ has a unique minimizer $x_{n}\in\mathrm{S}_{0}$ satisfying

[TABLE]

where the latter denotes positive definiteness.

Remark 2*.*

Assuming that $\min_{x\in\mathrm{S}_{0}}I_{n}(x)=0$ is just a particular (but helpful) normalization and in general not restrictive: If $\min_{x\in\mathrm{S}_{0}}I_{n}(x)=c>-\infty$, then we can simply set

[TABLE]

for which we obtain

[TABLE]

and $\min_{x\in\mathrm{S}_{0}}\hat{I}_{n}(x)=\min_{x}\hat{\Phi}_{n}(x)-\frac{1}{n}\log\pi_{0}(x)=0$.

Given Assumption 1 we define the Laplace approximation of $\mu_{n}$ as the following Gaussian measure

[TABLE]

Thus, we have

[TABLE]

where we can view

[TABLE]

as the second-order Taylor approximation $\tilde{I}_{n}=T_{2}I_{n}(x_{n})$ of $I_{n}$ at $x_{n}$. This point of view is crucial for analyzing the approximation

[TABLE]

Notation and recurring equations.

Before we continue, we collect recurring important definitions and where they can be found in Table 1 and provide the following important equations cheat sheet

[TABLE]

2.2 Convergence in Hellinger Distance

By a modification of Theorem 1 for integrals w.r.t. a weight $\operatorname{\mathrm{e}}^{-n\Phi_{n}(x)}$ we may show a corresponding version of (4), i.e., for sufficiently smooth $f\in L^{1}_{\mu_{0}}(\mathbb{R})$

[TABLE]

However, in this section we study a stronger notion of convergence of $\mathcal{L}_{\mu_{n}}$ to $\mu_{n}$, namely, w.r.t. the total variation distance $d_{\text{TV}}$ and the Hellinger distance $d_{\mathrm{H}}$. Given two probability measures $\mu$, $\nu$ on $\mathbb{R}^{d}$ and another probability measure $\rho$ dominating $\mu$ and $\nu$ the total variation distance of $\mu$ and $\nu$ is given by

[TABLE]

and their Hellinger distance by

[TABLE]

It holds true that

[TABLE]

see [17, Equation (8)]. Note that, $d_{\text{TV}}(\mu_{n},\mathcal{L}_{\mu_{n}})\to 0$ implies that $|\int f\mathrm{d}\mu_{n}-\int f\mathrm{d}\mathcal{L}_{\mu_{n}}|\to 0$ for any bounded and continuous $f\colon\mathbb{R}^{d}\to\mathbb{R}$. In order to establish our convergence result, we require almost the same assumptions as in Theorem 1, but now uniformly w.r.t. $n$:

Assumption 2**.**

There holds $\Phi_{n},\pi_{0}\in C^{3}(\mathrm{S}_{0},\mathbb{R})$ for all $n$ and

there exist the limits

[TABLE]

in $\mathbb{R}^{d}$ and $\mathbb{R}^{d\times d}$, respectively, with $H_{\star}$ being positive definite and $x_{\star}$ belonging to the interior of $\mathrm{S}_{0}$. 2. 2.

For each $r>0$ there exists an $n_{r}\in\mathbb{N}$, $\delta_{r}>0$ and $K_{r}<\infty$ such that

[TABLE]

as well as

[TABLE] 3. 3.

There exists a uniformly bounding function $q\colon\mathrm{S}_{0}\to[0,\infty)$ with

[TABLE]

for an $n_{0}\in\mathbb{N}$ such that $q^{1-\epsilon}$ is integrable, i.e., $\int_{\mathrm{S}_{0}}q^{1-\epsilon}(x)\ \mathrm{d}x<\infty$, for an $\epsilon\in(0,1)$.

The only additional assumptions in comparison to the classical convergence theorem of the Laplace method are about the third derivatives of $\pi_{0}$ and $\Phi_{n}$ and the convergence of $x_{n}\to x_{\star}$. We remark that (12) implies

[TABLE]

and, thus, also $\lim_{n\to\infty}C_{n}=H^{-1}_{\star}$.

The uniform lower bound on $I_{n}$ outside a ball around $x_{n}$ as well as the integrable majorant of the unnormalized densities $\operatorname{\mathrm{e}}^{-nI_{n}}\leq 1$ of $\mu_{n}$ can be understood as uniform versions of the first and second assumption of Theorem 1. The third item of Assumption 2 implies the uniform integrability of the $\operatorname{\mathrm{e}}^{-nI_{n}}\leq 1$ and is obviously satisfied for bounded supports $\mathrm{S}_{0}$. However, in the unbounded case it seems to be crucial111We provide a counterexample of nonconcentrating measures $\mu_{n}$ in the case, when the third issue of Assumptions 2 does not hold, in Appendix C. We also state in Appendix C.1 sufficient conditions on $\Phi_{n}$ and $\pi_{0}$ such that the third issue of Assumptions 2 is satisfied. for an increasing concentration of the $\mu_{n}$.

We start our analysis with the following helpful lemma.

Lemma 1**.**

Let Assumption 1 and 2 be satisfied and let $\pi_{n},\tilde{\pi}_{n}\colon\mathbb{R}^{d}\to[0,\infty)$ denote the unnormalized Lebesgue densities of $\mu_{n}$ and $\mathcal{L}_{\mu_{n}}$, respectively, given by

[TABLE]

and

[TABLE]

Then, for any $p\in\mathbb{N}$

[TABLE]

Proof.

We define the remainder term

[TABLE]

i.e., for $x\in\mathrm{S}_{0}$ we have $\frac{\pi_{n}(x)}{\tilde{\pi}_{n}(x)}=\exp(-nR_{n}(x))$. Moreover, note that for $x\in\mathrm{S}_{0}^{c}$ there holds $\pi_{n}(x)=0$. Thus, we obtain

[TABLE]

where we define for a given radius $r>0$

[TABLE]

In Appendix B.1 we prove that

[TABLE]

for $c_{r},c_{r,\epsilon}>0$, which then yields the statement. ∎

Lemma 1 provides the basis for our main convergence theorem.

Theorem 2**.**

Let the assumptions of Lemma 1 be satisfied. Then, there holds

[TABLE]

Proof.

We start with

[TABLE]

For the first term there holds due to Lemma 1

[TABLE]

For the second term on the right-hand side we obtain

[TABLE]

Furthermore, due to Lemma 1 there exists a $c<\infty$ such that

[TABLE]

This yields

[TABLE]

which concludes the proof. ∎

Convergence of other Gaussian approximations.

Let us consider now a sequence of arbitrary Gaussian approximations $\tilde{\mu}_{n}=\mathcal{N}(a_{n},\frac{1}{n}B_{n})$ to the measures $\mu_{n}$ in (1). Under which conditions on $a_{n}\in\mathbb{R}^{d}$ and $B_{n}\in\mathbb{R}^{d\times d}$ do we still obtain the convergence $d_{\text{H}}(\mu_{n},\tilde{\mu}_{n})\to 0$? Of course, $a_{n}\to x_{\star}$ seems to be necessary but how about the covariances $B_{n}$? Due to the particular scaling of $1/n$ appearing in the covariance of $\mathcal{L}_{\mu_{n}}$, one might suppose that for example $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}I_{d})$ or $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}B)$ with an arbitrary symmetric and positive definite (spd) $B\in\mathbb{R}^{d\times d}$ should converge to $\mu_{n}$ as $n\to\infty$. However, since

[TABLE]

and $d_{\text{H}}(\mu_{n},\mathcal{L}_{\mu_{n}})\to 0$, we have

[TABLE]

The following result shows that, in general, $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}I_{d})$ or $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}B)$ do not converge to $\mu_{n}$.

Theorem 3**.**

Let the assumptions of Lemma 1 be satisfied.

For $\tilde{\mu}_{n}\coloneqq\mathcal{N}(x_{n},\frac{1}{n}B_{n})$, $n\in\mathbb{N}$, with spd $B_{n}$, we have that

[TABLE]

If so and if $\|C_{n}-B_{n}\|\in\mathcal{O}(n^{-1})$, then we even have $d_{\mathrm{H}}(\mu_{n},\tilde{\mu}_{n})\in\mathcal{O}(n^{-1/2})$. 2. 2.

For $\tilde{\mu}_{n}\coloneqq\mathcal{N}(a_{n},\frac{1}{n}B_{n})$, $n\in\mathbb{N}$, with $B_{n}$ satisfying (14) and $\|x_{n}-a_{n}\|\in\mathcal{O}(n^{-1})$, we have that $d_{\mathrm{H}}(\mu_{n},{\tilde{\mu}_{n}})\in\mathcal{O}(n^{-1/2})$.

The proof is straightforward given the exact formula for the Hellinger distance of Gaussian measures and can be found in Appendix B.2. Thus, Theorem 3 tells us that, in general, the Gaussian measures $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}I_{d})$ do not converge to $\mu_{n}$ as $n\to\infty$ whereas it is easily seen that $\tilde{\mu}_{n}=\mathcal{N}(x_{n},\frac{1}{n}H_{\star})$, indeed, do converge.

Relation to the Bernstein–von Mises theorem in Bayesian inference.

The Bernstein–von Mises (BvM) theorem is a classical result in Bayesian inference and asymptotic statistics in $\mathbb{R}^{d}$ stating the posterior consistency under mild assumptions [41]. Its extension to infinite-dimensional situations does not hold in general [9, 15], but can be shown under additional assumptions [16, 3, 4, 28]. In order to state the theorem we introduce the following setting: let $Y_{i}\sim\nu_{x_{0}}$, $i\in\mathbb{N}$, be i.i.d. random variables on $\mathbb{R}^{D}$, $d\leq D$, following a distribution $\nu_{x_{0}}(\mathrm{d}y)=\exp(-\ell(y,x_{0}))\bm{1}_{\mathrm{S}_{y}}(y)\mathrm{d}y$ where $\mathrm{S}_{y}\subset\mathbb{R}^{D}$ and where $\ell\colon\mathrm{S}_{y}\times\mathbb{R}^{d}\to[-\ell_{\min},\infty)$ represents the negative log-likelihood function for observing $y\in\mathrm{S}_{y}$ given a parameter value $x\in\mathbb{R}^{d}$. Assuming a prior measure $\mu_{0}(\mathrm{d}x)=\pi_{0}(x)\bm{1}_{\mathrm{S}_{0}}(x)\ \mathrm{d}x$ for the unknown parameter, the resulting posterior after $n$ observations $y_{i}$ of the independent $Y_{i}$, $i=1,\ldots,n$, is of the form (1) with

[TABLE]

We will denote the corresponding posterior measure by $\mu_{n}^{y_{1},\ldots,y_{n}}$ in order to highlight the dependence of the particular data $y_{1},\ldots,y_{n}$. The BvM theorem states now the convergence of this posterior to a sequence of Gaussian measures. This looks very similar to the statement of Theorem 2. However, the difference lies in the Gaussian measures as well as the kind of convergence. In its usual form the BvM theorem states under similar assumptions as for Theorem 2 that there holds in the large data limit

[TABLE]

where $\mu^{Y_{1},\ldots,Y_{n}}_{n}$ is now a random measure depending on the $n$ independent random variables $Y_{1},\ldots,Y_{n}$ and where the convergence in probability is taken w.r.t. randomness of the $Y_{i}$. Moreover, $\hat{x}_{n}=\hat{x}_{n}(Y_{1},\ldots,Y_{n})$ denotes an efficient estimator of the true parameter $x_{0}\in\mathrm{S}_{0}$—e.g., the maximum-likelihood or MAP estimator—and $\mathcal{I}_{x_{0}}$ denotes the Fisher information at the true parameter $x_{0}$, i.e.,

[TABLE]

Now both, the BvM theorem and Theorem 2, state the convergence of the posterior to a concentrating Gaussian measure where the rate of concentration of the latter (or better: of its covariance) is of order $n^{-1}$. Furthermore, also the rate of convergence in the BvM theorem can be shown to be of order $n^{-1/2}$ [19]. However, the main differences are:

• •

The BvM states convergence in probability (w.r.t. the randomness of the $Y_{i}$) and takes as basic covariance the inverse expected Hessian of the negative log likelihood at the data generating parameter value $x_{0}$. Working with this quantity requires the knowledge of the true value $x_{0}$ and the covariance operator is obtained by marginalizing over all possible data outcomes $Y$. This Gaussian measure is not a practical tool to be used but rather a limiting distribution of a powerful theoretical result reconciling Bayesian and classical statistical theory. For this reason, the Gaussian approximation in the statement of the BvM theorem can be thought of as being a “prior” approximation (in the loosest meaning of the word). Usually, a crucial requirement is that the problem is well-specified meaning that $x_{0}$ is an interior point of the prior support $\mathrm{S}_{0}$—although there exist results for misspecified models, see [22]. Here, a BvM theorem is proven without the assumption that $x_{0}$ belongs to the interior of $\mathrm{S}_{0}$. However, in this case the basic covariance is not the Fisher information but the Hessian of the mapping $x\mapsto d_{\text{KL}}(\nu_{0}||\nu_{x})$ evaluated at its unique minimizer where $d_{\text{KL}}(\nu_{0}||\nu_{x})$ denotes the Kullback-Leibler divergence of the data distribution $\nu_{x}$ given parameter $x\in\mathrm{S}_{0}$ w.r.t. the true data distribution $\nu_{0}$.

• •

Theorem 2 states the convergence for given realizations $y_{i}$ and takes the Hessian of the negative log posterior density evaluated at the current MAP estimate $x_{n}$ and the current data $y_{1},\ldots,y_{n}$. This means that we do not need to know the true parameter value $x_{0}$ and we employ the actual data realization at hand rather than averaging over all outcomes. Hence, we argue that the Laplace approximation (as stated in this context) provides a “posterior” approximation converging to the Bayesian posterior as $n\to\infty$.

Also, we require that the limit $x_{\star}=\lim_{n\to\infty}x_{n}$ is an interior point of the prior support $\mathrm{S}_{0}$.

• •

From a numerical point of view, the Laplace approximation requires the computation of the MAP estimate and the corresponding Hessian at the MAP, whereas the BvM theorem employs the Fisher information, i.e. requires an expectation w.r.t. the observable data. Thus, the Laplace approximation is based on fixed and finite data in contrast to the BvM.

The following example illustrates the difference between the two Gaussian measures: Let $x_{0}\in\mathbb{R}$ be an unknown parameter. Consider $n$ measurements $y_{k}\in\mathbb{R}$, $k=1,\ldots,n$, where $y_{k}$ is a realization of

[TABLE]

with $\epsilon_{k}\sim N(0,\sigma^{2})$ i.i.d.. For the Bayesian inference we assume a prior $N(0,\tau^{2})$ on $x$. Then the Bayesian posterior is of the form $\mu_{n}(\mathrm{d}x)\propto\exp(-nI_{n}(x))$ where

[TABLE]

The MAP estimator $x_{n}$ is the Laplace approximation’s mean and can be computed numerically as a minimizer of $I_{n}(x)$. It can be shown that $x_{n}$ converges to $x_{\star}=x_{0}$ for almost surely all realizations $y_{k}$ of $Y_{k}$ due to the strong law of large numbers. Now we take the Hessian (w.r.t. $x$) of $I_{n}$,

[TABLE]

and evaluate it in $x_{n}$ to obtain the covariance of the Laplace approximation, and, thus,

[TABLE]

On the other hand we compute the Gaussian BvM approximation: The Fisher information is given as (recall that $\Phi$ is the loglikelihood term as defined above)

[TABLE]

and hence we get the Gaussian approximation

[TABLE]

Now we clearly see the difference between the two measures and how they will be asymptotically identical, since $x_{n}\to x_{\star}=x_{0}$ due to consistency, $\frac{1}{n}\sum_{k=1}^{n}y_{k}$ converging a.s. to $x_{0}^{3}$ due to the strong law of large numbers, and with the prior-dependent part vanishing for $n\to 0$.

Remark 3*.*

Having raised the issue whether the BvM approximation $\mathcal{N}(\hat{x}_{n},n^{-1}\mathcal{I}^{-1}_{x_{0}})$ or the Laplace one $\mathcal{L}_{\mu_{n}}$ is closer to a given posterior $\mu_{n}$, one can of course ask for the best Gaussian approximation of $\mu_{n}$ w.r.t. a certain distance or divergence. Thus, we mention [26, 31] where such a best approximation w.r.t. the Kullback-Leibler divergence is considered. The authors also treat the case of best Gaussian mixture approximations for multimodal distributions and state a BvM like convergence result for the large data (and small noise) limit. However, the computation of such a best approximation can become costly whereas the Laplace approximation can be obtained rather cheaply.

2.3 The case of singular Hessians

The assumption, that the Hessians $H_{n}=\nabla^{2}\Phi_{n}(x_{n})$ as well as their limit $H_{\star}$ are positive definite, is quite restrictive. For example, for Bayesian inference with more unknown parameters than observational information, this assumption is not satisfied. Hence, we discuss in this subsection the convergence of the Laplace approximation in case of singular Hessians $H_{n}$ and $H_{\star}$. Nonetheless, we assume throughout the section that Assumption 1 is satisfied. This yields that the Laplace approximation $\mathcal{L}_{\mu_{n}}$ is well-defined. This means in particular that we suppose a regularizing effect of the log prior density $\log\pi_{0}$ on the minimization of $I_{n}(x)=\Phi_{n}(x)-\frac{1}{n}\log\pi_{0}(x)$.

We first discuss necessary conditions for the convergence of the Laplace approximation and subsequently state a positive result for Gaussian prior measures $\mu_{0}$.

Necessary conditions.

Let us consider the simple case of $\Phi_{n}\equiv\Phi$, i.e., the probability measures $\mu_{n}$ are given by

[TABLE]

where we assume now that $\Phi\colon\mathrm{S}_{0}\to[c,\infty)$ with $c>-\infty$. Intuitively, $\mu_{n}$ should converge weakly to the Dirac measure $\delta_{\mathcal{M}_{\Phi}}$ on the set

[TABLE]

On the other hand, the associated Laplace approximations $\mathcal{L}_{\mu_{n}}$ will converge weakly to the Dirac measure $\delta_{\mathcal{M}_{\mathcal{L}}}$ in the affine subspace

[TABLE]

Hence, it is necessary for the convergence $\mathcal{L}_{\mu_{n}}\to\mu_{n}$ in total variation or Hellinger distance that $\mathcal{M}_{\Phi}=\mathcal{M}_{\mathcal{L}}$, i.e., that the set of minimizers of $\Phi$ is linear. In order to ensure the latter, we state the following.

Assumption 3**.**

Let $\mathcal{X}\subseteq\mathbb{R}^{d}$ be a linear subspace such that for a projection $\mathrm{P}_{\mathcal{X}}$ onto $\mathcal{X}$ there holds

[TABLE]

and let the restriction $\Phi_{n}\colon\mathcal{X}\to\mathbb{R}$ possess a unique and nondegenerate global minimum for each $n\in\mathbb{N}$.

For the case $\Phi_{n}=\Phi$ this assumption implies, that

[TABLE]

where $\mathcal{X}^{c}$ denotes a complementary subspace to $\mathcal{X}$, i.e., $\mathcal{X}\oplus\mathcal{X}^{c}=\mathbb{R}^{d}$ and $x_{\star}\in\mathcal{X}$ the unique minimizer of $\Phi$ over $\mathcal{X}$. Besides that, Assumption 3 also yields that $x^{\top}H_{n}x=0$ iff $x\in\mathcal{X}^{c}$. Hence, this also holds for the limit $H_{\star}=\lim_{n\to\infty}H_{n}$ and we obtain

[TABLE]

Moreover, since Assumption 3 yields

[TABLE]

where $x_{\mathcal{X}}\coloneqq\mathrm{P}_{\mathcal{X}}x$ and $x_{c}\coloneqq\mathrm{P}_{\mathcal{X}^{c}}x=x-x_{\mathcal{X}}$, the marginal of $\mu_{n}$ coincides with the marginal of $\mu_{0}$ on $\mathcal{X}^{c}$. Hence, the Laplace approximation can only converge to $\mu_{n}$ in total variation or Hellinger distance if this marginal is Gaussian. We, therefore, consider the special case of Gaussian prior measures $\mu_{0}$.

Remark 4*.*

Please note that, despite this to some extent negative result for the Laplace approximation for singular Hessians, the preconditioning of sampling and quadrature methods via the Laplace approximation may still lead to efficient algorithms in the small noise setting. The analysis of Laplace approximation-based sampling methods, as introduced in the next section, in the underdetermined case will be subject to future work.

Convergence for Gaussian prior $\mu_{0}$.

A useful feature of Gaussian prior measures $\mu_{0}$ is that the Laplace approximation possesses a convenient representation via its density w.r.t. $\mu_{0}$.

Proposition 1** (cf. [43, Proposition 1]).**

Let Assumption 1 be satisfied and $\mu_{0}$ be Gaussian. Then there holds

[TABLE]

where $T_{2}\Phi_{n}(\cdot;x_{n})$ denotes the Taylor polynomial of order 2 of $\Phi_{n}$ at the point $x_{n}\in\mathbb{R}^{d}$.

In fact, the representation (17) does only hold for prior measures $\mu_{0}$ with Lebesgue density $\pi_{0}\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ satisfying $\nabla^{3}\log\pi_{0}\equiv 0$.

Corollary 1**.**

Let Assumption 1 be satisfied and $\mu_{0}$ be Gaussian. Further, let Assumption 3 hold true and assume that the restriction $\Phi_{n}\colon\mathcal{X}\to\mathbb{R}$ and the marginal density $\pi_{0}$ on $\mathcal{X}$ satisfy Assumption 2 on $\mathcal{X}$. Then the approximation result of Theorem 2 holds.

Proof.

By using Proposition 1, we can express the Hellinger distance $d_{\text{H}}(\mu_{n},\mathcal{L}_{\mu_{n}})$ as follows

[TABLE]

We use now the decomposition $\mathbb{R}^{d}=\mathcal{X}\oplus\mathcal{X}^{c}$ with $x\coloneqq x_{\mathcal{X}}+x_{c}$ for $x\in\mathbb{R}^{d}$ with $x_{\mathcal{X}}\in\mathcal{X}$ and $x_{c}\in\mathcal{X}^{c}$. We note, that due to Assumption 3, we have that

[TABLE]

We then obtain by disintegration and denoting $\tilde{\Phi}_{n}(x)\coloneqq T_{2}\Phi_{n}(x;x_{n})=\tilde{\Phi}_{n}(x_{\mathcal{X}})$

[TABLE]

where $\mu_{0}(\mathrm{d}x_{\mathcal{X}})$ denotes the marginal of $\mu_{0}$ on $\mathcal{X}$. Since $\Phi_{n}$ and $I_{n}(x_{\mathcal{X}})=\Phi_{n}(x_{\mathcal{X}})-\frac{1}{n}\log\pi_{0}(x_{\mathcal{X}})$, where $\pi_{0}(x_{\mathcal{X}})$ denotes the Lebesgue density of the marginal $\mu_{0}(\mathrm{d}x_{\mathcal{X}})$, satisfy the assumptions of Theorem 2 on $\mathrm{S}_{0}\cap\mathcal{X}=\mathcal{X}$, the statement follows. ∎

We provide some illustrative examples for the theoretical results stated in this subsection.

Example 1** (Divergence of the Laplace approximation in the singular case).**

We assume a Gaussian prior $\mu_{0}=N(0,I_{2})$ on $\mathbb{R}^{2}$ and $\Phi(x)=\|y-F(x)\|^{2}$ where

[TABLE]

We plot the Lebesgue densities of the resulting $\mu_{n}$ and $\mathcal{L}_{\mu_{n}}$ for $n=128$ in the left and middle panel of Figure 1. The red line in both plots indicate the different sets

[TABLE]

around which $\mu_{n}$ and $\mathcal{L}_{\mu_{n}}$, respectively, concentrate as $n\to\infty$. As $\mathcal{M}_{\Phi}\neq\mathcal{M}_{\mathcal{L}}$, we observe no convergence of the Laplace approximation as $n\to\infty$, see the right panel of Figure 1. Here, the Hellinger distance is computed numerically by applying a tensorized trapezoidal rule on a suffieciently large subdomain of $\mathbb{R}^{2}$.

Example 2** (Convergence of the Laplace approximation in the singular case in the setting of Corollary 1).**

Again, we suppose a Gaussain prior $\mu_{0}=N(0,I_{2})$ and $\Phi$ in the form of $\Phi(x)=\|y-F(x)\|^{2}$ with

[TABLE]

Thus, the invariant subspace is $\mathcal{X}^{c}=\{x\in\mathbb{R}^{2}\colon x_{1}=x_{2}\}$. In the left and middle panel of Figure 2 we present the Lebesgue densities of $\mu_{n}$ and its Laplace approximation $\mathcal{L}_{\mu_{n}}$ for $n=25$ and by the red line the sets $\mathcal{M}_{\Phi}=\mathcal{M}_{\mathcal{L}}=x_{\star}+\mathcal{X}^{c}$. We observe the convergence guaranteed by Corollary 1 in the right panel of Figure 2 where we can also notice a preasymptotic phase with a shortly increasing Hellinger distance. Such a preasmyptotic phase is to be expected due to $d_{\text{H}}(\mu_{n},\mathcal{L}_{\mu_{n}})\in\mathcal{O}(n^{-1/2})+\mathcal{O}(\operatorname{\mathrm{e}}^{-n\delta_{r}}n^{d/2})$ as shown in the proof of Theorem 2.

3 Robustness of Laplace-based Monte Carlo Methods

In practice, we are often interested in expectations or integrals of quantities of interest $f\colon\mathbb{R}^{d}\to\mathbb{R}$ w.r.t. $\mu_{n}$ such as

[TABLE]

For example, in Bayesian statistics the posterior mean ($f(x)=x$) or posterior probabilities ($f(x)=\mathbf{1}_{A}(x)$, $A\in\mathcal{B}(\mathbb{R}^{d})$) are desirable quantities. Since $\mu_{n}$ is seldom given in explicit form, numerical integration must be applied for approximating such integrals. To this end, since the prior measure $\mu_{0}$ is typically a well-known measure for which efficient numerical quadrature methods are available, the integral w.r.t. $\mu_{n}$ is rewritten as two integrals w.r.t. $\mu_{0}$

[TABLE]

If then a quadrature rule such as $\int_{\mathbb{R}^{d}}g(x)\ \mu_{0}(\mathrm{d}x)\approx\frac{1}{N}\sum_{i=1}^{N}w_{i}\ g(x_{i})$ is used, we end up with an approximation

[TABLE]

This might be a good approximation for small $n\in\mathbb{N}$. However, as soon as $n\to\infty$ the likelihood term $\exp(-n\Phi_{n}(x_{i}))$ will deteriorate and this will be reflected by a deteriorating efficiency of the quadrature scheme—not in terms of the convergence rate w.r.t. $N$, but w.r.t. the constant in the error estimate, as we will display later in examples.

If the Gaussian Laplace approximation $\mathcal{L}_{\mu_{n}}$ of $\mu_{n}$ is used as the prior measure for numerical integration instead of $\mu_{0}$, we get the following approximation

[TABLE]

where $\pi_{n}$ and $\tilde{\pi}_{n}$ denote the unnormalized Lebesgue density of $\mu_{n}$ and $\mathcal{L}_{\mu_{n}}$, respectively. This time, we can not only apply well-known quadrature and sampling rules for Gaussian measures, but moreover, we also know due to Lemma 1, that the ratio $\frac{\pi_{n}(x)}{\tilde{\pi}_{n}(x)}$ converges in mean w.r.t. $\mathcal{L}_{\mu_{n}}$ to $1$. Hence, we do not expect a deteriorating efficiency of the numerical integration as $n\to\infty$. On the contrary, as we subsequently discuss for several numerical integration methods, their efficiency for a finite number of samples $N\in\mathbb{N}$ will even improve as $n\to\infty$ if they are based on the Laplace approximation $\mathcal{L}_{\mu_{n}}$.

For the sake of simplicity, we consider the simple case of $\Phi_{n}\equiv\Phi+\text{const}$ in the following presentation—nonetheless, the presented results can be extended to the general case given appropriate modifications of the assumptions. Thus, we consider probability measures $\mu_{n}$ of the form

[TABLE]

where we assume that $\Phi$ satisfies the assumptions of Theorem 1. However, when dealing with the Laplace approximation of $\mu_{n}$ and, particularly, with the ratios of the corresponding normalizing constants, it is helpful to use the following representation

[TABLE]

where $\iota_{n}\coloneqq\min_{x\in\mathrm{S}_{0}}\Phi(x)-\frac{1}{n}\log\pi_{0}(x)$ and $Z_{n}=\operatorname{\mathrm{e}}^{n\iota_{n}}\int_{\mathbb{R}^{d}}\operatorname{\mathrm{e}}^{-n\Phi(x)}\pi_{0}(x)\ \mathrm{d}x$. By this construction the resulting $I_{n}(x)\coloneqq\Phi_{n}(x)-\frac{1}{n}\log\pi_{0}(x)$ satisfies $I_{n}(x_{n})=0$ as required in Assumption 1 for the construction of the Laplace approximation $\mathcal{L}_{\mu_{n}}$. Note, that for $\Phi_{n}=\Phi-\iota_{n}$ the Assumption 1 and 2 imply the assumptions of Theorem 1 for $f=\pi_{0}$ and $p=0$.

Preliminaries.

Before we start analyzing numerical methods based on the Laplace approximation as their reference measure, we take a closer look at the details of the asymptotic expansion for integrals provided in Theorem 1 and their implications for expectations w.r.t. $\mu_{n}$ given in (22).

The coefficients: The proof of Theorem 1 in [44, Section IX.5] provides explicit expressions222There is a typo in [44, Section IX.5] stating that the sum in (23) is taken over all $\bm{\alpha}\in\mathbb{N}_{0}^{d}$ with $|\bm{\alpha}|=k$. for the coefficients $c_{k}\in\mathbb{R}$ in the asymptotic expansion

[TABLE]

namely—given that $f\in C^{2p+2}(D,\mathbb{R})$ and $\Phi\in C^{2p+3}(D,\mathbb{R})$—that

[TABLE]

where for $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})$ we have $|\bm{\alpha}|=\alpha_{1}+\cdots+\alpha_{d}$, $\bm{\alpha}!=\alpha_{1}!\cdots\alpha_{d}!$, $D^{\bm{\alpha}}=D^{\alpha_{1}}_{x_{1}}\cdots D^{\alpha_{d}}_{x_{d}}$ and

[TABLE]

with $h\colon\Omega\to U(x_{\star})$ being a diffeomorphism between $0\in\Omega\subset\mathbb{R}^{d}$ and a particular neighborhood $U(x_{\star})$ of $x_{\star}$ mapping $h(0)=x_{\star}$ and such that $\det(\nabla h(0))=1$. The diffeomorphism $h$ is specified by the well-known Morse’s Lemma and depends only on $\Phi$. In particular, if $\Phi\in C^{2p+3}(D,\mathbb{R})$, then $h\in C^{2p+1}(\Omega,U(x_{\star}))$. For the constants $\kappa_{\bm{\alpha}}=\kappa_{\alpha_{1}}\cdots\kappa_{\alpha_{d}}\in\mathbb{R}$ we have $\kappa_{\alpha_{i}}=0$ if $\alpha_{i}$ is odd and $\kappa_{\alpha_{i}}=(2/\lambda_{i})^{(\alpha_{i}+1)/2}\Gamma((\alpha_{i}+1)/2)$ otherwise with $\lambda_{i}>0$ denoting the $i$th eigenvalue of $H_{\star}=\nabla^{2}\Phi(x_{\star})$. Hence, we get

[TABLE] 2. 2.

The normalization constant of $\mu_{n}$: Theorem 1 implies that if $\pi_{0}\in C^{2}(\mathbb{R}^{d};\mathbb{R})$ and $\Phi\in C^{3}(\mathbb{R}^{d},\mathbb{R})$, then

[TABLE]

Hence, we obtain for the normalizing constant $Z_{n}$ in (22) that

[TABLE]

If we compare this to the normalizing constant $\tilde{Z}_{n}=n^{-d/2}\sqrt{\det(2\pi C_{n})}$ of its Laplace approximation we get

[TABLE]

We now show that

[TABLE]

First, we get due to $C_{n}\to H_{\star}^{-1}$ that {\color[rgb]{0,0,0}\sqrt{\det(C_{n})}\to\frac{1}{\sqrt{\det(H_{\star})}}} as $n\to\infty$. Moreover,

[TABLE]

Since $x_{n}\to x_{\star}$ continuity implies $\pi_{0}(x_{n})\to\pi_{0}(x_{\star})$ as $n\to\infty$. Besides that, the strong convexity of $\Phi$ in a neighborhood of $x_{\star}$—due to $\nabla^{2}\Phi(X_{\star})>0$ and $\Phi\in C^{3}(\mathbb{R}^{d},\mathbb{R})$—implies that for a $c>0$

[TABLE]

also known as Polyak–Łojasiewicz condition. Because of

[TABLE]

since $\nabla I_{n}(x_{n})=0$, we have that $|\Phi(x_{n})-\Phi(x_{\star})|\in\mathcal{O}(n^{-2})$, and hence,

[TABLE]

This yields (26). 3. 3.

The expectation w.r.t. $\mu_{n}$: The expectation of a $f\in L^{1}_{\mu_{0}}(\mathbb{R})$ w.r.t. $\mu_{n}$ is given by

[TABLE]

If $f,\pi_{0}\in C^{2}(\mathbb{R}^{d},\mathbb{R})$ and and $\Phi\in C^{3}(\mathbb{R}^{d},\mathbb{R})$, then we can apply the asymptotic expansion above to both integrals and obtain

[TABLE]

If $f,\pi_{0}\in C^{4}(\mathbb{R}^{d};\mathbb{R})$ and $\Phi\in C^{5}(\mathbb{R}^{d},\mathbb{R})$, then we can make this more precise by using the next explicit terms in the asymptotic expansions of both integrals, apply the rule for the division of (asymptotic) expansions (cf. [29, Section 1.8]) and obtain ${\bm{\mathbb{E}}}_{\mu_{n}}\left[f\right]=f(x_{\star})+\tilde{c}_{1}(f,\pi_{0})n^{-1}+\mathcal{O}(n^{-2})$ where $\tilde{c}_{1}(f,\pi_{0})=\frac{1}{c_{0}(\pi_{0})}c_{1}(f\pi_{0})-\frac{c_{1}(\pi_{0})}{c^{2}_{0}(\pi_{0})}c_{0}(f\pi_{0})$. 4. 4.

The variance w.r.t. $\mu_{n}$: The variance of a $f\in L^{2}_{\mu_{0}}(\mathbb{R})$ w.r.t. $\mu_{n}$ is given by

[TABLE]

If $f,\pi_{0}\in C^{2}(\mathbb{R}^{d};\mathbb{R})$ and $\Phi\in C^{3}(\mathbb{R}^{d},\mathbb{R})$, then we can exploit the result for the expectation w.r.t. $\mu_{n}$ from above and obtain

[TABLE]

If $f,\pi_{0}\in C^{4}(\mathbb{R}^{d},\mathbb{R})$ and $\Phi\in C^{5}(\mathbb{R}^{d},\mathbb{R})$, then a straightforward calculation—see Appendix B.5— using the explicit formulas for $c_{1}(f^{2}\pi_{0})$ and $c_{1}(f\pi_{0})$ as well as $\nabla h(0)=I$ yields

[TABLE]

Hence, the variance $\operatorname{Var}_{\mu_{n}}(f)$ decays like $n^{-1}$ provided that $\nabla f(x_{\star})\neq 0$—otherwise it decays (at least) like $n^{-2}$.

Remark 5*.*

As already exploited above, the assumptions of Theorem 1 imply that $\Phi$ is strongly convex in a neighborhood of $x_{\star}=\lim_{n\to\infty}x_{n}$, where $x_{n}=\operatornamewithlimits{argmin}_{x\in\mathrm{S}_{0}}\Phi(x)-\frac{1}{n}\log\pi_{0}(x)$. This yields $|\Phi(x_{n})-\Phi(x_{\star})|\in\mathcal{O}(n^{-2})$, and thus

[TABLE]

3.1 Importance Sampling

Importance sampling is a variant of Monte Carlo integration where an integral w.r.t. $\mu$ is rewritten as an integral w.r.t. a dominating importance distribution $\mu\ll\nu$, i.e.,

[TABLE]

The integral appearing on the righthand side is then approximated by Monte Carlo integration w.r.t. $\nu$: given $N$ independent draws $x_{i}$, $i=1,\ldots,N$, according to $\nu$ we estimate

[TABLE]

Often the density or importance weight function $w=\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\colon\mathbb{R}^{d}\to[0,\infty)$ is only known up to a normalizing constant $\tilde{w}\propto\frac{\mathrm{d}\mu}{\mathrm{d}\nu}$. In this case, we can use self-normalized importance sampling

[TABLE]

As for Monte Carlo, there holds a strong law of large numbers (SLLN) for self-normalized importance sampling, i.e.,

[TABLE]

where $X_{i}\sim\nu$ are i.i.d. , which follows from the ususal SLLN and the continuous mapping theorem. Moreover, by the classical central limit theorem (CLT) and Slutsky’s theorem also a similar statement holds for self-normalized importance sampling: given that

[TABLE]

we have

[TABLE]

Thus, the asymptotic variance $\sigma^{2}_{\mu,\nu}(f)$ serves as a measure of efficiency for self-normalized importance sampling. To ensure a finite $\sigma^{2}_{\mu,\nu}(f)$ for many functions of interest $f$, e.g., bounded $f$, the importance distribution $\nu$ has to have heavier tails than $\mu$ such that the ratio $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}$ belongs to $L^{2}_{\nu}(\mathbb{R})$, see also [32, Section 3.3]. Moreover, if we even have $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\in L^{\infty}_{\nu}(\mathbb{R})$ we can bound

[TABLE]

i.e., the ratio between the asymptotic variance of importance sampling w.r.t. $\nu$ and plain Monte Carlo w.r.t. $\mu$ can be bounded by the $L^{\infty}_{\nu}$- or supremum norm of the importance weight $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}$.

For the measures $\mu_{n}$ a natural importance distribution (called $\nu$ above) which allows for direct sampling are the prior measure $\mu_{0}$ and the Gaussian Laplace approximation $\mathcal{L}_{\mu_{n}}$. We study the behaviour of the resulting asymptotic variances $\sigma^{2}_{\mu_{n},\mu_{0}}(f)$ and $\sigma^{2}_{\mu_{n},\mathcal{L}_{\mu_{n}}}(f)$ in the following.

Prior importance sampling.

First, we consider $\mu_{0}$ as importance distribution. For this choice the importance weight function $w_{n}\coloneqq\frac{\mathrm{d}\mu_{n}}{\mathrm{d}\mu_{0}}$ is given by

[TABLE]

with $\Phi_{n}(x)=\Phi(x)-\iota_{n}$, see (22). Concerning the bound in (31) we immediately obtain for sufficiently smooth $\pi_{0}$ and $\Phi$ by (25), assuming w.l.o.g. $\min_{x}\Phi(x)=\Phi(x_{\star})=0$, that

[TABLE]

explodes as $n\to\infty$.

Of course, that is just the deterioration of an upper bound, but in fact we can prove the following rather negative result where we use the notation $g(n)\sim h(n)$ for the asymptotic equivalence of functions of $n$, i.e., $g(n)\sim h(n)$ iff $\lim_{n\to\infty}\frac{g(n)}{h(n)}=1$.

Lemma 2**.**

Given $\mu_{n}$ as in (22) with $\Phi$ satisfying the assumptions of Theorem 1 for $p=1$ and $\pi_{0}\in C^{4}(\mathbb{R}^{d},\mathbb{R})$ with $\pi_{0}(x_{\star})\neq 0$, we have for any $f\in C^{4}(\mathbb{R}^{d},\mathbb{R})\cap L^{1}_{\mu_{0}}(\mathbb{R})$ with $\nabla f(x_{\star})\neq 0$ that

[TABLE]

which yields $\frac{\sigma^{2}_{\mu_{n},\mu_{0}}(f)}{\operatorname{Var}_{\mu_{n}}(f)}\sim\tilde{c}_{f}n^{d/2}$ for another $\tilde{c}_{f}>0$.

Proof.

W.l.o.g. we may assume that $f(x_{\star})=0$, since $\sigma^{2}_{\mu_{n},\mu_{0}}(f)=\sigma^{2}_{\mu_{n},\mu_{0}}(f-c)$ for any $c\in\mathbb{R}$. Moreover, for simplicity we assume w.l.o.g. that $\Phi(x_{\star})=0$. We study

[TABLE]

by analyzing the growth of the numerator and denominator w.r.t. $n$. Due to the preliminaries presented above we know that $\operatorname{\mathrm{e}}^{-2n\iota_{n}}Z^{2}_{n}=c_{0}^{2}n^{-d}+\mathcal{O}(n^{-d-1})$ with c_{0}={\color[rgb]{0,0,0}(2\pi)^{d/2}\,\pi_{0}(x_{\star})/\sqrt{\det(H_{\star})}}>0. Concerning the numerator we start with decomposing

[TABLE]

where this time

[TABLE]

We derive now asymptotic expansions of these terms based on Theorem 1. It is easy to see that the assumptions of Theorem 1 are also fulfilled when considering integrals w.r.t. $\operatorname{\mathrm{e}}^{-2n\Phi}$. We start with $J_{1}$ and obtain due to $f(x_{\star})=0$ that

[TABLE]

where $c^{\prime}_{1}(f^{2}\pi_{0})\in\mathbb{R}$ is the same as $c_{1}(f^{2}\pi_{0})$ in (23) but for $2\Phi$ instead of $\Phi$.

Next, we consider $J_{2}$ and recall that due to $f(x_{\star})=0$ we have ${\bm{\mathbb{E}}}_{\mu_{n}}\left[f\right]\in\mathcal{O}(n^{-1})$, see (27). Furthermore, $f(x_{\star})=0$ also implies $\int_{\mathrm{S}_{0}}f(x)\pi_{0}(x)\,\operatorname{\mathrm{e}}^{-2n\Phi(x)}\ \mathrm{d}x\in\mathcal{O}(n^{-1-d/2})$, see Theorem 1. Thus, we have

[TABLE]

Finally, we take a look at $J_{3}$. By Theorem 1 we have $\int_{\mathrm{S}_{0}}\exp(-2n\Phi(x))\ \mu_{0}(\mathrm{d}x)\in\mathcal{O}(n^{-d/2})$ and, hence, obtain

[TABLE]

Hence, $J_{1}$ has the dominating power w.r.t. $n$ and we have that

[TABLE]

At this point, we remark that due to the assumption $\nabla f(x_{\star})\neq 0$ we have $c^{\prime}_{1}(f^{2}\pi_{0})\neq 0$: we know by (24) that $c^{\prime}_{1}(f^{2}\pi_{0})=\frac{1}{2}\sum_{j=1}^{d}\kappa_{2\bm{e}_{j}}D^{2\bm{e}_{j}}F(0)$ where $F(x)=\pi_{0}(h^{\prime}(x))f^{2}(h^{\prime}(x))\det(\nabla h^{\prime}(x))$ and $h^{\prime}$ denotes the diffeomorphism for $2\Phi$ appearing in Morse’s lemma and mapping [math] to $x_{\star}$; applying the product formula and using $f(x_{\star})=0$ as well as $\det(\nabla h^{\prime}(0))=1$ we get that $D^{2\bm{e}_{j}}F(x_{\star})=\pi_{0}(x_{\star})D^{2\bm{e}_{j}}(f^{2}(h^{\prime}(x_{\star})))$; similarly, we get using $f(x_{\star})=0$ that $D^{2\bm{e}_{j}}(f^{2}(h^{\prime}(x_{\star})))=2|\bm{e}_{j}^{\top}\nabla h^{\prime}(0)\nabla f(x_{\star})|^{2}$; since $h^{\prime}$ is a diffeomorphishm $\nabla h^{\prime}(0)$ is regular and, thus, $c^{\prime}_{1}(f^{2}\pi_{0})\neq 0$. The statement follows now by

[TABLE]

and by recalling that $\operatorname{Var}_{\mu_{n}}(f)\sim cn^{-1}$ because of $\nabla f(x_{\star})\neq 0$, see (29). ∎

Thus, Lemma 2 tells us that the asymptotic variance of importance sampling for $\mu_{n}$ with the prior $\mu_{0}$ as importance distribution grows like $n^{d/2-1}$ as $n\to\infty$ for a large class of integrands. Hence, its efficiency deteriorates like $n^{d/2-1}$ for $d\geq 3$ as the target measures $\mu_{n}$ become more concentrated.

Laplace-based importance sampling.

We now consider the Laplace approximation $\mathcal{L}_{\mu_{n}}$ as importance distribution which yields the following importance weight function

[TABLE]

with $R_{n}(x)=I_{n}(x)-\tilde{I}_{n}(x)=I_{n}(x_{n})-\frac{1}{2}\|x-x_{n}\|^{2}_{C_{n}^{-1}}$ for $x\in\mathrm{S}_{0}$. In order to ensure $w_{n}\in L^{2}_{\mathcal{L}_{\mu_{n}}}(\mathbb{R})$ we need that

[TABLE]

Despite pathological counterexamples a sound requirement for $w_{n}\in L^{2}_{\mathcal{L}_{\mu_{n}}}(\mathbb{R})$ is that

[TABLE]

for example by assuming that there exist $\delta,c_{1}>0$, $c_{0}>0$, and $n_{0}\in\mathbb{N}$ such that

[TABLE]

If the Lebesgue density $\pi_{0}$ of $\mu_{0}$ is bounded, then (33) is equivalent to the existence of $n_{0}$ and a $\tilde{c}_{0}$ such that

[TABLE]

Unfortunately, condition (33) is not enough to ensure a well-behaved asymptotic variance $\sigma^{2}_{\mu_{n},\mathcal{L}_{\mu_{n}}}(f)$ as $n\to\infty$, since

[TABLE]

Although, we know due to (26) that $\frac{\tilde{Z}_{n}}{Z_{n}}\to 1$ as $n\to\infty$, the supremum norm of the importance weight $w_{n}$ of Laplace-based importance sampling will explode exponentially with $n$ if $\min_{x}R_{n}(x)<0$. This can be sharpened to proving that even the asymptotic variance of Laplace-based importance sampling w.r.t. $\mu_{n}$ as in (22) deteriorates exponentially as $n\to\infty$ for many functions $f\colon\mathbb{R}^{d}\to\mathbb{R}$ if

[TABLE]

by means of Theorem 1 applied to

[TABLE]

This means, except when $\Phi$ is basically strongly convex, the asymptotic variance of Laplace-based important sampling can explode exponentially or not even exist as $n$ increases. However, in the good case, so to speak, we obtain the following.

Proposition 2**.**

Consider the measures $\mu_{n}$ as in (22) with $\Phi_{n}=\Phi-\iota_{n}$ and $\pi_{0}$ satisfying Assumption 1 and 2. If there exist an $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$ we have

[TABLE]

then for any $f\in L^{2}_{\mu_{0}}(f)$

[TABLE]

Proof.

The assumption (34) ensures that $R_{n}(x)=I_{n}(x)-\tilde{I}_{n}(x)\geq 0$ for each $x\in\mathrm{S}_{0}$. Thus,

[TABLE]

and the assertion follows by (31) and the fact that $\lim_{n\to\infty}\frac{\tilde{Z}_{n}}{Z_{n}}=1$ due to (26). ∎

Condition (34) is for instance satisfied, if $I_{n}$ is strongly convex with a constant $\gamma\geq\lambda_{\min}(\nabla^{2}I_{n}(x_{n}))$ where the latter denotes the smallest eigenvalue of the positive definite Hessian $\nabla^{2}I_{n}(x_{n})$. However, this assumption or even (34) is quite restrictive and, probably, hardly fulfilled for many interesting applications. Moreover, the success in practice of Laplace-based importance sampling is well-documented. How come that despite a possible infinite asymptotic variance Laplace-based importance sampling performs that well? In the following we refine our analysis and exploit the fact that the Laplace approximation concentrates around the minimizer of $I_{n}$. Hence, with an increasing probability samples drawn from the Laplace approximation are in a small neighborhood of the minimizer. Thus, if $I_{n}$ is, e.g., only locally strongly convex—which the assumptions of Theorem 2 actually imply—then with a high probability the mean squared error might be small. We clarify these arguments in the following and present a positive result for Laplace-based importance sampling under mild assumptions but for a weaker error criterion than the decay of the mean squared error.

First we state a concentration result for $N$ samples drawn from $\mathcal{L}_{\mu_{n}}$ which is an immediate consequence of Proposition 4.

Proposition 3**.**

Let $N\in\mathbb{N}$ be arbitrary and let $X^{(n)}_{i}\sim\mathcal{L}_{\mu_{n}}$ be i.i.d. where $i=1,\ldots,N$. Then, for a sequence of radii $r_{n}\geq r_{0}n^{-q}>0$, $n\in\mathbb{N}$, with $q\in(0,1/2)$ we have

[TABLE]

Remark 6*.*

In the following we require expectations w.r.t. restrictions of the measures $\mu_{n}$ in (22) to shrinking balls $B_{r_{n}}(x_{n})$. To this end, we note that the statements of Theorem 1 also hold true for shrinking domains $D_{n}=B_{r_{n}}(x_{\star})$ with $r_{n}=r_{0}n^{-q}$ as long as $q<1/2$. This can be seen from the proof of Theorem 1 in [44, Section IX.5]. In particular, all coefficients in the asymptotic expansion for $\int_{D_{n}}f(x)\exp(-n\Phi(x))\mathrm{d}x$ with sufficiently smooth $f$ are the same as for $\int_{D}f(x)\exp(-n\Phi(x))\mathrm{d}x$ and the difference between both integrals decays for increasing $n$ like $\exp(-cn^{\epsilon})$ for an $\epsilon>0$ and $c>0$. Concerning the balls $B_{r_{n}}(x_{n})$ with decaying radii $r_{n}=r_{0}n^{-q}$, $q\in[0,1/2)$, we have due to $\|x_{n}-x_{\star}\|\in\mathcal{O}(n^{-1})$—see Remark 5—that $B_{r_{n}/2}(x_{\star})\subset B_{r_{n}}(x_{n})\subset B_{2r_{n}}(x_{\star})$ for sufficiently large $n$. Thus, the facts for $\mu_{n}$ as in (22) stated in the preliminaries before Section 3.1 do also apply to the restrictions of $\mu_{n}$ to $B_{r_{n}}(x_{n})$ with $r_{n}=r_{0}n^{-q}$, $q\in[0,1/2)$. In particular, the difference between ${\bm{\mathbb{E}}}_{\mu_{n}}\left[f\right]$ and ${\bm{\mathbb{E}}}_{\mu_{n}}\left[f\ |\ B_{r_{n}}(x_{n})\right]$ decays faster than any negative power of $n$ as $n\to\infty$.

The next result shows that the mean absolute error of the Laplace-based importance sampling behaves like $n^{-(3q-1)}$ conditional on all $N$ samples belonging to shrinking balls $B_{r_{n}}(x_{n})$ with $r_{n}=r_{0}n^{-q}$, $q\in(1/3,1/2)$.

Lemma 3**.**

Consider the measures $\mu_{n}$ in (22) and suppose they satisfy the assumptions of Theorem 2. Then, for any $f\in C^{2}(\mathbb{R}^{d},\mathbb{R})\cap L^{2}_{\mu_{0}}(\mathbb{R})$ there holds for the error

[TABLE]

of the Laplace-based importance sampling with $N\in\mathbb{N}$ samples that

[TABLE]

where $r_{n}=r_{0}n^{-q}$ with $q\in(1/3,1/2)$.

Proof.

We start with

[TABLE]

The second term decays subexponentially w.r.t. $n$, see Remark 6. Hence, it remains to prove that

[TABLE]

To this end, we write the self-normalizing Laplace-based importance sampling estimator as

[TABLE]

where we define

[TABLE]

and recall that $w_{n}$ is as in (32) and $\tilde{w}_{n}(x)=\exp(-nR_{n}(x))$. Notice that

[TABLE]

Let us denote the event that $X^{(n)}_{1},\ldots,X^{(n)}_{N}\in B_{r_{n}}(x_{n})$ by $A_{n,N}$ for brevity. Then,

[TABLE]

The first term in the last line can be bounded by the conditional variance of $S_{n,N}$ given $X^{(n)}_{1},\ldots,X^{(n)}_{N}\in B_{r_{n}}(x_{n})$, i.e., by Jensen’s inequality we obtain

[TABLE]

see Remark 6 and the preliminaries before Subsection 3.1. Thus,

[TABLE]

and it remains to study if ${\bm{\mathbb{E}}}\left[\left.\left|\left(H_{n,N}-1\right)S_{n,N}\right|\ \right|\ A_{n,N}\right]\in\mathcal{O}(n^{-(3q-1)})$. Given that $X^{(n)}_{1},\ldots,X^{(n)}_{N}\in B_{r_{n}}(x_{n})$ we can bound the values of the random variable $H_{n,N}$ for sufficiently large $n$: first, we have $Z_{n}/\tilde{Z}_{n}=1+\mathcal{O}(n^{-1})$, see (26), and second

[TABLE]

Since $|R_{n}(x)|\leq c_{3}\|x-x_{n}\|^{3}$ for $|x-x_{n}|\leq r_{n}$ due to the local boundedness of the third derivative of $I_{n}$ and $r_{n}=r_{0}n^{-q}$, we have that

[TABLE]

where $c>0$. Thus, there exist $\alpha_{n}\leq 1\leq\beta_{n}$ with $\alpha_{n}=\operatorname{\mathrm{e}}^{-cn^{1-3q}}(1+\mathcal{O}(n^{-1}))$ and $\beta_{n}\sim\operatorname{\mathrm{e}}^{cn^{1-3q}}(1+\mathcal{O}(n^{-1}))$ such that

[TABLE]

Since $\operatorname{\mathrm{e}}^{\pm cn^{1-3q}}(1+\mathcal{O}(n^{-1}))=1\pm cn^{1-3q}+\mathcal{O}(n^{-1})$ we get that for sufficiently large $n$ there exists a $\tilde{c}>0$ such that

[TABLE]

Hence,

[TABLE]

since ${\bm{\mathbb{E}}}\left[|S_{n,N}|\ |\ A_{n,N}\right]\leq{\bm{\mathbb{E}}}_{\mu_{n}}\left[|f|\ |\ B_{r_{n}(x_{n})}\right]$ is uniformly bounded w.r.t. $n$. This concludes the proof. ∎

We now present our main result for the Laplace-based importance sampling which states that the corresponding error decays in probability to zero as $n\to\infty$ and the order of decay is arbitrary close to $n^{-1/2}$.

Theorem 4**.**

Let the assumptions of Lemma 3 be satisfied. Then, for any $f\in C^{2}(\mathbb{R}^{d},\mathbb{R})\cap L^{2}_{\mu_{0}}(\mathbb{R})$ and each sample size $N\in\mathbb{N}$ the error $e_{n,N}(f)$ of Laplace-based importance sampling satisfies

[TABLE]

Proof.

Let $0\leq\delta<1/2$ and $\epsilon>0$ be arbitrary. We need to show that

[TABLE]

Again, let us denote the event that $X^{(n)}_{1},\ldots,X^{(n)}_{N}\in B_{r_{n}}(x_{n})$ by $A_{n,N}$ for brevity. By Proposition 3 we obtain for radii $r_{n}=r_{0}n^{-q}$ with $q\in(1/3,1/2)$ that

[TABLE]

The second term on the righthand side in the last line obviously tends to $1$ exponentially as $n\to\infty$. Thus, it remains to prove that

[TABLE]

To this end, we apply a conditional Markov inequality for the positive random variable $e_{n,N}(f)$, i.e.,

[TABLE]

where we used Lemma 3. Choosing $q\in(1/3,1/2)$ such that $q>\frac{1+\delta}{3}\in[1/3,1/2)$ yields the statement. ∎

3.2 Quasi-Monte Carlo Integration

We now want to approximate integrals as in (20) w.r.t. measures $\mu_{n}(\mathrm{d}x)\propto\exp(-n\Phi(x))\mu_{0}(\mathrm{d}x)$ as in (22) by Quasi-Monte Carlo methods. These will be used to estimate the ratio $Z^{\prime}_{n}/Z_{n}$ by separately approximating the two integrals $Z^{\prime}_{n}$ and $Z_{n}$ in (20). The preconditioning strategy using the Laplace approximation will be explained exemplarily for Gaussian and uniform priors, two popular choices for Bayesian inverse problems.

We start the discussion by first focusing on a uniform prior distribution $\mu_{0}=\mathcal{U}([-\frac{1}{2},\frac{1}{2}]^{d})$. The integrals $Z^{\prime}_{n}$ and $Z_{n}$ are then

[TABLE]

where we set $\Theta_{n}(x)\coloneqq\exp(-n\Phi(x))$ for brevity.

We consider Quasi-Monte Carlo integration based on shifted Lattice rules: an $N$-point Lattice rule in the cube $[-\frac{1}{2},\frac{1}{2}]^{d}$ is based on points

[TABLE]

where $z\in\{1,\ldots,N-1\}^{d}$ denotes the so-called generating vector, $\Delta$ is a uniformly distributed random shift on $[-\frac{1}{2},\frac{1}{2}]^{d}$ and $\operatorname{frac}$ denotes the fractional part (component-wise). These randomly shifted points provide unbiased estimators

[TABLE]

of the two integrals $Z^{\prime}_{n}$ and $Z_{n}$ in (35). Under the assumption that the quantity of interest $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is linear and bounded, we can focus in the following on the estimation of the normalization constant $Z_{n}$, the results can be then straightforwardly generalized to the estimation of $Z^{\prime}_{n}$. For the estimator $Z_{n,QMC}$ we have the following well-known error bound.

Theorem 5**.**

[12, Thm. 5.10]** Let $\gamma=\{\gamma_{\bm{\nu}}\}_{{\bm{\nu}}\subset\{1,\ldots,d\}}$ denote POD (product and order dependent) weights of the form $\gamma_{\bm{\nu}}=\alpha_{|{\bm{\nu}}|}\prod_{j\in{\bm{\nu}}}\beta_{j}$ specified by two sequences $\alpha_{0}=\alpha_{1}=1,\alpha_{2},\ldots\geq 0$ and $\beta_{1}\geq\beta_{2}\geq\ldots>0$ for ${\bm{\nu}}\subset\{1,\ldots,d\}$ and $|{\bm{\nu}}|=\#{\bm{\nu}}$. Then, a randomly shifted Lattice rule with $N=2^{m},m\in\mathbb{N}$, can be constructed via a component-by-component algorithm with POD weights at costs of $\mathcal{O}(dN\log N+d^{2}N)$ operations, such that for sufficiently smooth $\Theta\colon[-\frac{1}{2},\frac{1}{2}]^{d}\to[0,\infty)$

[TABLE]

for $\kappa\in(1/2,1]$ with

[TABLE]

and $\zeta(a)\coloneqq\sum_{k=1}^{\infty}k^{-a}$.

The norm $\|\Theta_{n}\|_{\gamma}$ in the convergence analysis depends on $n$, in particular, it can grow polynomially w.r.t. the concentration level $n$ of the measures $\mu_{n}$ as we state in the next result.

Lemma 4**.**

Let $\Phi\colon\mathbb{R}^{d}\to[0,\infty)$ satisfy the assumptions of Theorem 1 for $p=2d$. Then, for the norm $\|\Theta_{n}\|_{\gamma}$ in the error bound in Theorem 5 there holds

[TABLE]

The proof of Lemma 4 is rather technical and can be found in Appendix B.3. We remark that Lemma 4 just tells us that the root mean squared error estimate for QMC integration based on the prior measure explodes like $n^{d/4}$. This does in general not indicate that the error itself explodes; in fact the QMC integration error for the normalization constant is bounded by $1$ in our setting. Nonetheless, Lemma 4 indicates that a naive Quasi-Monte Carlo integration based on the uniform prior $\mu_{0}$ is not suitable for highly concentrated target or posterior measures $\mu_{n}$. We subsequently propose and study a Quasi-Monte Carlo integration based on the Laplace approximation $\mathcal{L}_{\mu_{n}}$.

Laplace-based Quasi-Monte Carlo.

To stabilize the numerical integration for concentrated $\mu_{n}$, we propose a preconditioning based on the Laplace approximation, i.e., an affine rescaling according to the mean and covariance of $\mathcal{L}_{\mu_{n}}$. In the uniform case, the functionals $I_{n}$ are independent of $n$. The computation of the Laplace approximation requires therefore only one optimization to solve for $x_{n}=x_{\star}=\operatornamewithlimits{argmin}_{x\in[-\frac{1}{2},\frac{1}{2}]^{d}}\Phi(x)$. In particular, the Laplace approximation of $\mu_{n}$ is given by $\mathcal{L}_{\mu_{n}}=\mathcal{N}(x_{\star},\frac{1}{n}H_{\star}^{-1})$ where $H_{\star}$ denotes the positive definite Hessian $\nabla^{2}\Phi(x_{\star})$. Hence, $H_{\star}$ allows for an orthogonal diagonalization $H_{\star}=QDQ^{\top}$ with orthogonal matrix $Q\in\mathbb{R}^{d\times d}$ and diagonal matrix $D=\operatorname{diag}(\lambda_{1},\ldots\lambda_{d})\in\mathbb{R}^{d\times d}$, $\lambda_{1}\geq\ldots\geq\lambda_{d}>0$.

We now use this decomposition in order to construct an affine transformation which reverses the increasing concentration of $\mu_{n}$ and yields a QMC approach robust w.r.t. $n$. This transformation is given by

[TABLE]

where $\tau\in(0,1)$ is a truncation parameter. The idea of the transformation $g_{n}$ is to zoom into the parameter domain and thus, to counter the concentration effect. The domain will then be truncated to $G_{n}\coloneqq g_{n}([-\frac{1}{2},\frac{1}{2}]^{d})\subset[-\frac{1}{2},\frac{1}{2}]^{d}$ and we consider

[TABLE]

The determinant of the Jacobian of the transformation $g_{n}$ is given by $\det(\nabla g_{n}(x))\equiv C_{\text{trans},n}=\left(\frac{2|\ln\tau|}{n}\right)^{\frac{d}{2}}\sqrt{\det(H_{\star})}\sim c_{\tau}n^{-d/2}$. We will now explain how the parameter $\tau$ effects the truncation error. For given $\tau\in(0,1)$, the Laplace approximation is used to determine the truncation effect:

[TABLE]

Thus, since due to the concentration effect of the Laplace approximation we have $\int_{\mathrm{S}_{0}}\ \mathcal{L}_{\mu_{n}}(\mathrm{d}x)\to 1$ exponentially with $n$, we get

[TABLE]

thus, the truncation error $\int_{\mathrm{S}_{0}\setminus G_{n}}\ \mathcal{L}_{\mu_{n}}(\mathrm{d}x)$ becomes arbitrarly small for sufficiently small $\tau\ll 1$, since $\mathrm{erf}(t)\to 1$ as $t\to 1$. If we apply now QMC integration using shifted Lattice rule in order to compute the integral over $[-\frac{1}{2},\frac{1}{2}]^{d}$ on the righthand side of (38), we obtain the following estimator for $\hat{Z}_{n}$ in (38):

[TABLE]

with $x_{i}$ as in (36). Concerning the norm $\|\Theta_{n}\circ g_{n}\|_{\gamma}$ appearing in the error bound for $|\hat{Z}_{n}-\hat{Z}_{n,QMC}|$ we have now the following result.

Lemma 5**.**

Let $\Phi\colon\mathbb{R}^{d}\to[0,\infty)$ satisfy the assumptions of Theorem 1 for $p=2d$. Then, for the norm $\|\Theta_{n}\circ g_{n}\|_{\gamma}$ with $g_{n}$ as above there holds

[TABLE]

Again, the proof is rather technical and can be found in Appendix B.4. This proposition yields now our main result.

Corollary 2**.**

Given the assumptions of Lemma 5, a randomly shifted lattice rule with $N=2^{m},m\in\mathbb{N}$, can be constructed via a component-by-component algorithm with product and order dependent weights at costs of $\mathcal{O}(dN\log N+d^{2}N)$ operations, such that for $\kappa\in(1/2,1]$

[TABLE]

with constants $c_{1},c_{2},c_{3}>0$ independent of $n$ and $h(\tau)=1-\mathrm{erf}(0.5\sqrt{|\ln\tau}|)^{d}$.

Proof.

The triangle inequality leads to a separate estimation of the domain truncation error of the integral w.r.t. the posterior and the QMC approximation error, i.e.

[TABLE]

The second term on the right hand side corresponds to the QMC approximation error. Thus, Theorem 5 and Lemma 5 imply

[TABLE]

where the term $n^{-\frac{d}{2}}$ is due to $C_{\text{trans},n}\sim c_{\lambda}n^{-\frac{d}{2}}$. The domain truncation error can be estimated similar to the proof of Lemma 1:

[TABLE]

where $\tilde{Z}_{n}=n^{-\frac{d}{2}}\sqrt{\det(2\pi H_{\star}^{-1})}$. The result follows by the proof of Lemma 1. ∎

Remark 7*.*

In the case of Gaussian priors, the transformation simplifies to $w=g_{n}(x)=x_{\star}+n^{\frac{1}{2}}QD^{-\frac{1}{2}}x$ due to the unboundedness of the prior support. However, to show an analogous result to Corollary 2, uniform bounds w.r. to $n$ on the norm of the mixed first order derivatives of the preconditioned posterior density $\Theta_{n}(g_{n}(T^{-1}x))$ in a weighted Sobolev space, where $T^{-1}$ denotes the inverse cumulative distribution function of the normal distribution, need to be proven. See [23] for more details on the weighted space setting in the Gaussian case. Then, a similar result to Corollary 2 follows straightforwardly from [23, Thm 5.2]. The numerical experiments shown in subsection 3.3 suggest that we can expect a noise robust behavior of Laplace-based QMC methods also in the Gaussian case. This will be subject to future work.

Remark 8*.*

Note that the QMC analysis in Theorem 5 can be extended to an infinite dimensional setting, cp. [23] and the references therein for more details. This opens up the interesting possibility to generalize the above results to the infinite dimensional setting and to develop methods with convergence independent of the number of parameters and independent of the measurement noise. Furthermore, higher order QMC methods can be used for cases with smooth integrands, cp. [13, 11, 10], leading to higher convergence rates than the first order methods discussed here. In the uniform setting, it has been shown in [38] that the assumptions on the first order derivatives (and also higher order derivatives) of the transformed integrand are satisfied for Bayesian inverse problems related to a class of parametric operator equations, i.e., the proposed approach leads to a robust performance w.r.t. the size of the measurement noise for integrating w.r.t. posterior measure resulting from this class of forward problems. The theoretical analysis of this setting will be subject to future work.

Remark 9* (Numerical Quadrature).*

Higher regularity of the integrand allows to use higher order methods such as sparse quadrature and higher order QMC methods, leading to faster convergence rates. In the infinite dimensional Bayesian setting with uniform priors, we refer to [36, 37] for more details on sparse quadrature for smooth integrands. In the case of uniform priors, the methodology introduced above can be used to bound the quadrature error for the preconditioned integrand by the truncation error and the sparse grid error.

3.3 Examples

In this subsection we present two examples illustrating our previous theoretical results for importance sampling and quasi-Monte Carlo integration based on the prior measure $\mu_{0}$ and the Laplace approximation $\mathcal{L}_{\mu_{n}}$ of the target measure $\mu_{n}$. Both examples are Bayesian inference or inverse problems where the first one uses a toy forwad map and the second one is related to inference for a PDE model.

3.3.1 Algebraic Bayesian inverse problem

We consider inferring $x\in[-\frac{1}{2},\frac{1}{2}]^{d}$ for $d=1,2,3,4$ based on a uniform prior $\mu_{0}=\mathcal{U}([-\frac{1}{2},\frac{1}{2}]^{d})$ and a realisation of the noisy observable of $Y=F(X)+\varepsilon_{n}$ where $X\sim\mu_{0}$ and the noise $\varepsilon_{n}\sim N(0,n^{-1}\Gamma_{d})$, $\Gamma_{d}=0.1I_{d}$, are independent, and $F(x)=(F_{1}(x),\ldots,F_{d}(x))$ with

[TABLE]

for $x=(x_{1},\ldots,x_{d})$. The resulting posterior measure $\mu_{n}$ on $[-\frac{1}{2},\frac{1}{2}]^{d}$ are of the form (22) with

[TABLE]

We used $y=G(0.25\cdot{\bm{1}}_{d})$ throughout where ${\bm{1}}_{d}=(1,\ldots,1)\in\mathbb{R}^{d}$. We then compute the posterior expectation of the quantity of interest $f(x)=x_{1}+\cdots+x_{d}$. To this end, we employ importance sampling and quasi-Monte Carlo integration based on $\mu_{0}$ and the Laplace approximation $\mathcal{L}_{\mu_{n}}$ as outlined in the precious subsections. We compare the output of these methods to a reference solution obtained by a brute-force tensor grid trapezoidal rule for integration. In particular, we estimate the root mean squared error (RMSE) of the methods and how it evolves as $n$ increases.

Results for importance sampling: In order to be sufficiently close to the asymptotic limit, we use $N=10^{5}$ samples for self-noramlized importance sampling. We run $1,000$ independent simulations of the importance sampling integration and compute the resulting empirical RMSE. In Figure 3 we present the results for increasing $n$ and various $d$ for prior-based and Laplace-based importance sampling. We obtain a good match to the theoretical results, i.e., the RMSE for choosing the prior measure as importance distribution behaves like $n^{d/4-1/2}$ in accordance to Lemma 2. Besides that the RMSE for choosing the Laplace approximation as importance distribution decays like $n^{1/2}$ after a preasymptotic phase. This is relates to the statement of Theorem 4 where we have shown that the absolute error333We have also computed the empirical $L^{1}$-error which showed a similar behaviour as the RMSE. decays in probability like $n^{1/2}$. Note that the assumptions of Proposition 2 are not satisfied for this example for all $d=1,2,3,4$.

Results for quasi-Monte Carlo: We use $N=2^{10}$ quasi-Monte Carlo points for prior- and Laplace-based QMC. For the Laplace-case we employ a truncation parameter of $\tau=10^{-6}$ and discard all transformed points outside of the domain $[-\frac{1}{2},\frac{1}{2}]^{d}$. Again, we run $1,000$ random shift simulations for both QMC methods and estimate the empirical RMSE. However, for QMC we report the relative RMSE, since, for example, the decay of the normalization constant $Z_{n}\in\mathcal{O}(n^{-d/2})$ dominates the growth of the absolute error of prior QMC integration for the normalization constant. In Figure 4 and 5 we display the resulting relative RMSE for the quantity related integral $Z_{n}^{\prime}$, the normalization constant $Z_{n}$, i.e.,

[TABLE]

and the resulting ratio $\frac{Z^{\prime}_{n}}{Z_{n}}$ for increasing $n$ and various $d$ for prior-based and Laplace-based QMC. For prior-based QMC we observe for dimensions $d\geq 2$ a algebraic growth of the relative error. In the previous subsection we have proven that the corresponding classical error bound will explode which does not necessary imply that the error itself explodes—as we can see for $d=1$. However, this simple example shows that also the error will often grow algebraically with increasing $n$. For the Laplace-based QMC we observe on the other hand in Figure 5 a decay of the relative empirical RMSE. By Corollary 2 we can expect that the relative errors stay bounded as $n\to\infty$. This provides motivation for a further investigation. In particular, we will analysize the QMC ratio estimator for $\frac{Z^{\prime}_{n}}{Z_{n}}$ in a future work.

3.3.2 Bayesian inference for an elliptic PDE

In the following we illustrate the preconditioning ideas from the previous section by Bayesian inference with differential equations. To this end we consider the following model parametric elliptic problem

[TABLE]

with $f(t)=100\cdot t$, $t\in[0,1]$, and diffusion coefficient

[TABLE]

where $\psi_{k}(t)=\frac{0.1}{k}\sin(k\pi t)$ and the $x_{k}\in\mathbb{R}$, $k=1,\ldots,d$, are to be inferred by noisy observations of the solution $q$ at certain points $t_{j}\in[0,1]$. For $d=1,2$ these observations are taken at $t_{1}=0.25$ and $t_{2}=0.75$ and for $d=3$ they are taken at $t_{j}\in\{0.125,0.25,0.375,0.6125,0.75,0.875\}$. We suppose an additive Gaussian observational noise $\epsilon\sim\mathcal{N}(0,\Gamma_{n})$ with noise covariance $\Gamma_{n}=n^{-1}\Gamma_{obs}$ and $\Gamma_{obs}\in\mathbb{R}^{2\times 2}$ or $\Gamma_{obs}\in\mathbb{R}^{7\times 7}$, respectively, specified later on. In the following we place a uniform and a Gaussian prior $\mu_{0}$ on $\mathbb{R}^{d}$ and would like to integrate w.r.t. the resulting posterior $\mu_{n}$ on $\mathbb{R}^{d}$ which is of the form (22) with

[TABLE]

where $F\colon\mathbb{R}^{d}\to\mathbb{R}^{2}$ for $d=1,2$, and $F\colon\mathbb{R}^{d}\to\mathbb{R}^{7}$ for $d=3$, respectively, denotes the mapping from the coefficients $x\coloneqq(x_{k})_{k=1}^{d}$ to the observations of the solution $q$ of the elliptic problem above and the vector $y\in\mathbb{R}^{2}$ or $y\in\mathbb{R}^{7}$, respectively, denotes the observational data resulting from $y=F(x)+\epsilon$ with $\epsilon$ as above. Our goal is then to compute the posterior expectation (i.e., w.r.t. $\mu_{n}$) of the following quantity of interest $f\colon\mathbb{R}^{d}\to\mathbb{R}$: $f(u)$ is the value of the solution $q$ of the elliptic problem at $t=0.5$.

Uniform prior.

We place a uniform prior $\mu_{0}=\mathcal{U}([-\frac{1}{2},\frac{1}{2}]^{d})$ for $d=1,2$ or $3$ and choose $\Gamma_{obs}=0.01I_{2}$ for $d=1,2$ and $\Gamma_{obs}=0.01I_{7}$ for $d=3$. We display the resulting posteriors $\mu_{n}$ for $d=2$ in Figure 6 illustrating the concentration effect of the posterior for various values of the noise scaling $n$ and the resulting transformed posterior with $\Phi\circ g_{n}$ based on Laplace approximation. The truncation parameter is set to $\tau=10^{-6}$. We observe the almost quadratic behavior of the preconditioned posterior, as expected from the theoretical results.

We are now interested in the numerical performance of the Importance Sampling and QMC method based on the prior distribution compared to the performance of the preconditioned versions based on Laplace approximation. The QMC estimators are constructed by an off-the-shelf generating vector (order-2 randomly shifted weighted Sobolev space), which can be downloaded from https://people.cs.kuleuven.be/~dirk.nuyens/qmc-generators/ (exod2_base2_m20_CKN.txt). The reference solution used to estimate the error is based on a (tensor grid) trapezoidal rule with $10^{6}$ points in 1D, $4\cdot 10^{6}$ points in 2D in the original domain, i.e., the truncation error is quantified and in the transformed domain in 3D with $10^{6}$ points. Figure 7 illustrates the robust behavior of the preconditioning strategy w.r.t. the scaling $1/n$ of the observational noise. Though we know from the theory that in the low dimensional case (1D, 2D), the importance sampling method based on the prior is expected to perform robust, we encounter numerical instabilities due to the finite number of samples used for the experiments. The importance sampling results are based on $10^{6}$ sampling points, the QMC results on $8192$ shifted lattice points with $2^{6}$ random shifts.

Figure 8 shows the RMSE of the normalization constant $Z_{n}$ using the preconditioned QMC approach with respect to the noise scaling $1/n$. We observe a numerical confirmation of the predicted dependence of the error w.r.t. the dimension (cp. Corollary 2).

We remark that the numerical experiments for the ratio suggest a behavior $n^{-1/2}$, i.e., a rate independent of the dimension $d$, of the RMSE for the preconditioned QMC approach, cp. Figure 7. This will be subject to future work.

Gaussian prior.

We choose as prior $\mu_{0}=\mathcal{N}(0,I_{2})$ for the coefficients $x=(x_{1},x_{2})\in\mathbb{R}^{2}$ for $\hat{u}_{2}$ in the elliptic problem above. For the noise covariance we set this time $\Gamma_{obs}=I_{2}$. The performance of the prior based and preconditioned version of Importance Sampling is depicted in Figure 9. Clearly, the Laplace approximation as a preconditioner improves the convergence behavior; we observe a robust behavior w.r.t. the noise level.

The convergence of the QMC approach is depicted in Fig 10, showing a consistent behavior with the considerations in the previous section.

4 Conclusions and Outlook to Future Work

This paper makes a number of contributions in the development of numerical methods for Bayesian inverse problems, which are robust w.r.t. the size of the measurement noise or the concentration of the posterior measure, respectively. We analyzed the convergence of the Laplace approximation to the posterior distribution in Hellinger distance. This forms the basis for the design of variance robust methods. In particular, we proved that Laplace-based importance sampling behaves well in the small noise or large data size limit, respectively. For uniform priors, Laplace-based QMC methods have been developed with theoretically and numerically proven errors decaying with the noise level or concentration of the measure, repestively. Some future directions of this work include the development of noise robust Markov chain Monte Carlo methods and the combination of dimension independent and noise robust strategies. This will require the study of the Laplace approximation in infinite dimensions in a suitable setting. Finally, we could study in more details the error in the ratio estimator using Laplace-based QMC methods. The use of higher order QMC methods has been proven to be a promising direction for a broad class of Bayesian inverse problems and the design of noise robust versions is an interesting and potentially fruitful research direction.

Acknowledgements.

The authors are grateful to the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Uncertainty quantification for complex systems: theory and methodologies when work on this paper was undertaken. This work was supported by: EPSRC grant numbers EP/K032208/1 and EP/R014604/1. Moreover, the authors would like to thank the anonymous referees and Remo Kretschmann for helpful comments which improved the paper significantly. Furthermore, BS was supported by the DFG project 389483880 and by the DFG RTG1953 ”Statistical Modeling of Complex Systems and Processes”.

Appendix A Concentration of the Laplace Approximation

Due to the well-known Borell-TIS inequality for Gaussian measures on Banach spaces [24, Chapter 3] we obtain the following useful concentration result for the Laplace approximation.

Proposition 4**.**

Let Assumption 1 and Assumption 2 be satisfied. Then, for any $r>0$ there exists an $c_{r}>0$ such that

[TABLE]

where $B^{c}_{r}(x_{n})=\{x\in\mathbb{R}^{d}\colon\|x-x_{n}\|>r\}$.

Proof.

Let $X_{n}\sim\mathcal{N}\left(0,n^{-1}C_{n}\right)$. Then

[TABLE]

We now use the well-known concentration of Gaussian measures, namely,

[TABLE]

where $\sigma_{n}^{2}\coloneqq\sup_{\|x\|\leq 1}{\bm{\mathbb{E}}}\left[|x^{\top}X_{n}|^{2}\right]$, see [24, Chapter 3]. There holds $\sigma_{n}^{2}=\lambda_{\text{max}}(n^{-1}C_{n})=\|n^{-1}C_{n}\|$ and we get

[TABLE]

Due to Assumption 2, i.e., $C_{n}\to H^{-1}_{\star}>0$, there exists a finite constant $0<c$ such that $\|C_{n}\|\geq c$ for all $n\in\mathbb{N}$. Analogously, there exist a constant $K<\infty$ such that $\mathrm{tr}\,(C_{n})\leq K$ for all $n$. The latter implies

[TABLE]

Hence, for an arbitrary $r$ let $n_{0}$ be such that ${\bm{\mathbb{E}}}\left[\|X_{n}\|\right]\leq r/2$ for all $n\geq n_{0}$, which yields

[TABLE]

∎

Appendix B Proofs

We collect the rather technical proofs in this appendix.

B.1 Proof of Lemma 1

Recall that we want to bound

[TABLE]

where $R_{n}(x)\coloneqq I_{n}(x)-\tilde{I}_{n}(x)=I_{n}(x)-I_{n}(x_{n})-\frac{1}{2}\|x-x_{n}\|^{2}_{C^{-1}_{n}}$ and $r>0$ is an at the moment arbitrary radius which will be specified in the following first paragraph.

Bounding $J_{1}$.

Due to $\Phi_{n},\log\pi_{0}\in C^{3}(\mathrm{S}_{0};\mathbb{R})$, we have that for any $x\in B_{r}(x_{n})$ there exists a $\xi_{x,n}\in B_{r}(x_{n})$ such that

[TABLE]

Moreover, since $x_{n}\to x_{\star}$ there exists an $0<r_{0}<\infty$ such that

[TABLE]

Hence, the local uniform boundedness of $\|\nabla^{3}I_{n}(\cdot)\|$, see Assumption 2, implies the existence of a finite $c_{3}>0$ such that for sufficiently large $n$, i.e., $n\geq n_{r}$, we have

[TABLE]

Thus, we obtain, due to $|\operatorname{\mathrm{e}}^{-t}-1|=1-\operatorname{\mathrm{e}}^{-t}\leq\operatorname{\mathrm{e}}^{t}-1$ for $t\geq 0$,

[TABLE]

which yields

[TABLE]

Now, since $C^{-1}_{n}\to H_{\star}>0$, there exists for sufficiently large $n$ a $\gamma>0$ such that

[TABLE]

Hence, for $x\in B_{r}(x_{n})$, i.e., $\|x-x_{n}\|\leq r$, we get

[TABLE]

By choosing $r\coloneqq\frac{\gamma}{2c_{3}}$ we obtain further

[TABLE]

Let us now introduce the auxiliary Gaussian measure $\nu_{n}\coloneqq\mathcal{N}\left(0,\frac{1}{n\gamma}I\right)$ with which we get

[TABLE]

There holds now

[TABLE]

due to the continuity of the determinant and $H_{\star}>0$. Moreover, since $1-\operatorname{\mathrm{e}}^{-t}\leq t$ for $t\geq 0$ we obtain with $\xi\sim\mathcal{N}(0,I)$ that

[TABLE]

This yields $J_{1}(n)\in\mathcal{O}(n^{-p/2})$ for the particular choice $r=\frac{\gamma}{2c_{3}}$. In the following two paragraphs we will use exactly this particular radius.

Bounding $J_{0}$.

Due to Assumption 2, we have $x_{n}\to x_{\star}$ as $n\to\infty$. Hence, there exists an $n_{0}<\infty$ such that $\|x_{n}-x_{\star}\|\leq r/2$ for all $n\geq n_{0}$. This implies by Assumption 2

[TABLE]

and, hence, $\mathrm{S}_{0}^{c}\subseteq B^{c}_{r/2}(x_{n})$. By Proposition 4, we obtain

[TABLE]

with a $0<c_{r/2}<\infty$.

Bounding $J_{2}$.

We divide the set $B_{r}^{c}(x_{n})\cap\mathrm{S}_{0}$ into several subsets in order to bound $J_{2}(n)$. First, we define

[TABLE]

i.e., $B^{c}_{r}(x_{n})\cap\mathrm{S}_{0}=\mathrm{P}_{n}\dot{\cup}\mathrm{P}^{c}_{n}$, and notice that

[TABLE]

Hence, due to Propostion 4 there exists a $c_{r}\in(0,\infty)$ such that

[TABLE]

[TABLE]

We now prove that

[TABLE]

with $\epsilon\in(0,1)$ as in Assumption 2. To this end, we observe that due to $I_{n}(x)\geq\delta_{r}$ for all $x\in\mathrm{P}^{c}_{n}\subset B^{c}_{r}(x_{n})$ the functions

[TABLE]

converge pointwise to zero as $n\to\infty$. Moreover, $I_{n}(x)\geq\delta_{r}$ yields

[TABLE]

with the bounding function $q$ introduced in Assumption 2. Thus, since $q^{1-\epsilon}$ is integrable we obtain by Lebesgue’s dominated convergence theorem

[TABLE]

as $n\to\infty$. Hence, (41) holds. Since $\widetilde{Z}_{n}\in\mathcal{O}(n^{-d/2})$ we get in summary,

[TABLE]

with $c_{r,\epsilon}:=\min\{c_{r},\epsilon\delta_{r}\}>0$.

B.2 Proof of Theorem 3

A straightforward calculation, see also [30, Exercise 1.14], yields that for Gaussian measures $\mathcal{N}_{a,Q}\coloneqq\mathcal{N}(a,Q)$, $\mathcal{N}_{b,Q}\coloneqq\mathcal{N}(b,Q)$ and $\mathcal{N}_{a,R}\coloneqq\mathcal{N}(a,R)$ we have

[TABLE]

By (13) we obtain for $\tilde{\mu}_{n}\coloneqq\mathcal{N}(x_{n},\frac{1}{n}B_{n})$ that

[TABLE]

Due to the local Lipschitz continuity of the determinant and $\|C_{n}^{-1}-H_{\star}\|\to 0$ we obtain the first statement for $\tilde{\mu}_{n}\coloneqq\mathcal{N}(x_{n},\frac{1}{n}B_{n})$. Furthermore, by the triangle inequality

[TABLE]

and exploiting the local Lipschitz continuity of $f(x)=\frac{1}{x}$ and of the determinant, we obtain

[TABLE]

with a generic constant $c>0$. Moreover, due to the local Lipschitz continuity of the square root of a matrix, we get

[TABLE]

where we used that $\|C_{n}^{-1/2}\|\to\|H_{\star}^{1/2}\|$. Furthermore, we get analogously that $\|I-C_{n}^{1/2}B_{n}^{-1/2}\|\leq c\|B_{n}-C_{n}\|$ using that $\|B_{n}^{-1/2}\|\to\|H_{\star}^{1/2}\|$. Thus, the second statement of the first item follows by

[TABLE]

The second item follows by applying the triangle inequality, expressing the Hellinger distance between $\mathcal{N}(x_{n},\frac{1}{n}B_{n})$ and $\mathcal{N}(a_{n},\frac{1}{n}B_{n})$ by

[TABLE]

and estimating

[TABLE]

where we used the fact that the spd matrices $B_{n}$ converge to the spd matrix $H_{\star}$, hence, the sequence of the smallest eigenvalue of $B_{n}$ is bounded away from zero.

B.3 Proof of Lemma 4

For the following proof we use the famous Faa di Bruno-formula for higher order derivatives of compositions given in [18]. To this end, let $v\colon\mathbb{R}^{d}\to\mathbb{R}$ and $u\colon\mathbb{R}\to\mathbb{R}$ be sufficiently smooth functions and define $w\coloneqq u\circ v$, i.e., $w\colon\mathbb{R}^{d}\to\mathbb{R}$. For a subset ${\bm{\nu}}\subset\{1,\ldots,d\}$ we consider the partial derivative $\frac{\partial^{|{\bm{\nu}}|}w}{\partial x_{{\bm{\nu}}}}$ where we set

[TABLE]

We then obtain (see [18])

[TABLE]

where $\Pi({\bm{\nu}})$ denotes the set of all partitions $P$ of the set ${\bm{\nu}}\subset\{1,\ldots,d\}$, $B\in P$ refers to running through the elements or blocks of the partition $P$ with $|B|$ denoting the cardinality of $B\subset{\bm{\nu}}$ and $|P|$ the number of blocks in $P$, and the same notational convention for $\partial x_{B}$ as above, i.e.,

[TABLE]

By the application of the multivariate Faa di Bruno-formula to

[TABLE]

we obtain for ${\bm{\nu}}\subset\{1,\ldots,d\}$ that

[TABLE]

By setting $x_{-{\bm{\nu}}}\coloneqq x_{\{1,\ldots,d\}\setminus{\bm{\nu}}}$ we get

[TABLE]

where we shortened $\Pi_{d}\coloneqq\Pi(\{1,\ldots,d\})$. We will now investigate, how—i.e., to which power w.r.t. $n$—

[TABLE]

decays as $n\to\infty$. To this end, we write

[TABLE]

and apply in the following Laplace’s method in order to derive the asymptotics of

[TABLE]

Since in the considered setting of a uniform prior on $[-\frac{1}{2},\frac{1}{2}]^{d}$ we have $I_{n}(x)=\Phi(x)$ for $x\in[-\frac{1}{2},\frac{1}{2}]^{d}$, there holds that $x_{n}=x_{\star}$, $\Phi(x_{n})=0$, and, by construction, also $\nabla\Phi(x_{n})=0$. The latter may cause a faster decay of $\int_{[-\frac{1}{2},\frac{1}{2}]^{d}}h_{P,\tilde{P}}\operatorname{\mathrm{e}}^{-2n\Phi}\mathrm{d}x$ than the usual $n^{-d/2}$ depending on the partitions $P,\tilde{P}$. For example, let $P=\tilde{P}=\{\{1\},\ldots,\{d\}\}$, i.e., $P$ and $\tilde{P}$ consist only of single blocks $B=\{i\}$, $i=1,\ldots,d$, then

[TABLE]

Exploiting (24) for the coefficients in the asymptotic expansion of $\int_{[-\frac{1}{2},\frac{1}{2}]^{d}}h_{P,\tilde{P}}\operatorname{\mathrm{e}}^{-2n\Phi}\mathrm{d}x$ one can calculate that for these particular partitions $P=\tilde{P}=\{\{1\},\ldots,\{d\}\}$ we have

[TABLE]

but

[TABLE]

since $\frac{\partial^{2}}{\partial x^{2}_{j}}\Phi(x_{\star})\neq 0$ due to $\nabla^{2}\Phi(x_{\star})$ being positive definite. Hence, for these partitions $|P|=|\tilde{P}|=d$ we get

[TABLE]

We can extend this reasoning to arbitrary partitions $P,\tilde{P}\in\Pi_{d}$. To this end, let $|P|_{1}\coloneqq|\{B\in P\colon|B|=1\}|$ denote the number of single blocks $|B|=1$ in $P\in\Pi_{d}$. Then, we know by the definition of $h_{P,\tilde{P}}$ that $h_{P,\tilde{P}}$ posseses a zero of order $|P|_{1}+|\tilde{P}|_{1}$ in $x_{\star}$. This in turn, implies that the first $\left\lfloor\frac{|P|_{1}+|\tilde{P}|_{1}}{2}\right\rfloor$ coefficients in the asymptotic expansion of $\int_{[-\frac{1}{2},\frac{1}{2}]^{d}}h_{P,\tilde{P}}\operatorname{\mathrm{e}}^{-2n\Phi}\mathrm{d}x$ are zero, hence,

[TABLE]

Thus, for arbitrary $P,\tilde{P}\in\Pi_{d}$ we have

[TABLE]

If we maximize the exponent on the righthand side we get that

[TABLE]

where the maximum is obtained, e.g., for the above choice of $P=\tilde{P}=\{\{1\},\ldots,\{d\}\}$. This means that certain summands in $F_{d}(n)$ grow like $n^{d/2}$ whereas the other ones grow slower w.r.t. $n$. Thus, we get that $F_{d}(n)\sim cn^{d/2}$ which yields the statement.

B.4 Proof of Lemma 5

Since the transformation

[TABLE]

is linear, we have for $j\in\{1,\ldots,d\}$

[TABLE]

where $A_{\cdot j}$ denotes the $j$th column of $A$. Thus, for a ${\bm{\nu}}\subset\{1,\ldots,d\}$ we get

[TABLE]

where the last term on the righthand side denotes the application of the multilinear form $(\nabla^{|{\bm{\nu}}|}\Phi)(g_{n}(x))\colon\mathbb{R}^{d\times|{\bm{\nu}}|}\to\mathbb{R}$ to the $|{\bm{\nu}}|$ arguements $A_{\cdot\nu_{j}}\in\mathbb{R}^{d}$. To keep the notation short, we denote by $\nabla^{|{\bm{\nu}}|}\Phi(g_{n}(x))[A_{\bm{\nu}}]$ the term $(\nabla^{|{\bm{\nu}}|}\Phi)(g_{n}(x))[A_{\cdot\nu_{1}},\ldots,A_{\cdot\nu_{|{\bm{\nu}}|}}]$. By Faa di Bruno we obtain now for any ${\bm{\nu}}\subset\{1,\ldots,d\}$ that

[TABLE]

which yields

[TABLE]

where

[TABLE]

Note, that we can bound

[TABLE]

where we set the “norm” of the multilinear form $\nabla^{|B|}\Phi(g_{n}(x))\colon\mathbb{R}^{d\times|B|}\to\mathbb{R}$ as

[TABLE]

Setting $c_{A}\coloneqq\prod_{j=1}^{d}\|A_{j}\|$ we get $\left|\prod_{B\in P}\nabla^{|B|}\Phi(g_{n}(x))[A_{B}]\right|\leq c_{A}\prod_{B\in P}\left\|\nabla^{|B|}\Phi(g_{n}(x))\right\|$ and obtain by multiplication—omitting the integral domains for a moment—that

[TABLE]

where we used the Cauchy–Schwarz and Jensen’s inequality in the second line. We apply Laplace’s method in order to examine the $L^{2}$-norm above:

[TABLE]

where $h_{P}=\prod_{B\in P}\left\|\left(\nabla^{|B|}\Phi\right)\circ g_{n}\right\|$. By the substitution $y\coloneqq g_{n}(x)=x_{\star}+n^{-1/2}Ax$ we get

[TABLE]

where $C_{J_{Trans}}(n)=\det(n^{-1/2}A)=n^{-d/2}\det(A)$. Since also $2\Phi$ satisfies the assumptions of Theorem 1 and $\Phi(x_{\star})=0$ we obtain

[TABLE]

However, since $\nabla\Phi(x_{\star})=0$ there holds $h_{P}^{2}(x_{\star})=0$ if there exists a single block $|B|=1$ in $P$ which then yields to a decay of the $L^{2}$-norm as least as fast as $n^{-1}$. In particular, denoting by $|P|_{1}$ the number of single blocks in $P$ we obtain by the same reasoning as in Section B.3 that

[TABLE]

where $c_{P}=c_{|P|_{1}}(h^{2}_{P})$ as in (23). Hence,

[TABLE]

Similarly to Section B.3 we can derive $\max_{P\in\Pi({\bm{\nu}})}|P|-|P|_{1}/2=|{\bm{\nu}}|/2$ which yields that $F_{\bm{\nu}}(n)\in\mathcal{O}(1)$ as $n\to\infty$ and concludes the proof.

B.5 Proof of Equation (29)

We show in the following that for $f,\pi_{0}\in C^{4}(\mathbb{R}^{d},\mathbb{R})$ and $\Phi\in C^{5}(\mathbb{R}^{d},\mathbb{R})$ we have

[TABLE]

To this end, we use $\operatorname{Var}_{\mu_{n}}(f)={\bm{\mathbb{E}}}_{\mu_{n}}\left[f^{2}\right]-{\bm{\mathbb{E}}}_{\mu_{n}}\left[f\right]^{2}$ and the asymptotics

[TABLE]

where $\tilde{c}_{1}(f,\pi_{0})=\frac{1}{c_{0}(\pi_{0})}c_{1}(f\pi_{0})-\frac{c_{1}(\pi_{0})}{c^{2}_{0}(\pi_{0})}c_{0}(f\pi_{0})$, analogously for $\tilde{c}_{1}(f^{2},\pi_{0})$. Thus, we have

[TABLE]

Since $\operatorname{Var}_{\mu_{n}}(f)=\operatorname{Var}_{\mu_{n}}(f-\mathrm{const})$ we may assume w.l.o.g. that $f(x_{\star})=0$. Hence,

[TABLE]

where the last equality follows due to $c_{0}(f^{2}\pi_{0})\propto f^{2}(x_{\star})=0$. We recall that

[TABLE]

where

[TABLE]

with $h\colon\Omega\to U(x_{\star})$ denoting the diffeomorphism between $0\in\Omega\subset\mathbb{R}^{d}$ and a particular neighborhood $U(x_{\star})$ of $x_{\star}$ mapping $h(0)=x_{\star}$ and such that $\det(\nabla h(0))=1$. Moreover, as outlined in [44, Section IX.2], we have that $\nabla h(0)^{\top}H_{\star}\nabla h(0)=\Lambda$ where $\Lambda:=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{d})$ denotes the diagonal matrix containing the eigenvalues $\lambda_{i}$ of $H_{\star}$. This is equivalent to

[TABLE]

which we will use later on. Furthermore, since $\kappa_{\bm{\alpha}}=\kappa_{\alpha_{1}}\cdots\kappa_{\alpha_{d}}\in\mathbb{R}$ and $\kappa_{\alpha_{i}}=(2/\lambda_{i})^{(\alpha_{i}+1)/2}\Gamma((\alpha_{i}+1)/2)$ for even $\alpha_{i}$, we get

[TABLE]

since $\kappa_{\bm{0}}=\det(2\pi H^{-1}_{\star})^{1/2}$. Thus,

[TABLE]

We now compute $\frac{\partial^{2}}{\partial x_{i}^{2}}\widetilde{F}_{0}(0)$. To this end, let us introduce $G_{0}(x):=\pi_{0}(h(x))\,\det(\nabla h(x))$, i.e., $\widetilde{F}_{0}(x)=f^{2}(h(x))G_{0}(x)$ as well as $G_{0}(0)=\pi_{0}(x_{\star})$. We first derive by the product and chain rule

[TABLE]

For the second derivative $\frac{\partial^{2}}{\partial x_{i}^{2}}\widetilde{F}_{0}(x)$ evaluated at $x=0$ we can exploit once more that w.l.o.g. $f(x_{\star})=0$ which then yields

[TABLE]

i.e., all other terms which we obtain by applying the product (and chain) rule are zero, since they contain $f(h(0))=f(x_{\star})=0$ as a factor. In summary, we end up with

[TABLE]

This concludes the proof.

Appendix C Visualization of the necessity of the third term in Assumption 2

Consider the following setting: Let $\pi_{0}(x)\propto\exp(-x^{2})$, $x\in\mathbb{R}$, (with suitable normalization) and let $\Phi_{n}\colon\mathbb{R}\to\mathbb{R}$ be a smooth function satisfying the following conditions:

• •

For some $0<r_{0}<\frac{1}{2}$ we have $\Phi_{n}(x)=x^{2}$ for $|x|\leq r_{0}$.

• •

For some $0<\delta<1$ and $r_{n}:=\sqrt{n}>r_{0}$ as well as $r_{n}^{+}:=r_{n}+\exp(n\delta)$ we have for all $r_{n}<|x|\leq r_{n}^{+}$ that

[TABLE]

• •

For a $r^{\prime}_{n}>r_{n}^{+}$ we have $\Phi_{n}(x)\equiv 0$ for all $|x|\geq r^{\prime}_{n}$.

• •

On $r_{0}<|x|<r_{n}$ and $r_{n}^{+}<|x|<r^{\prime}_{n}$ the function $\Phi$ is smooth and monotone.

This means that $I_{n}(x)=\Phi_{n}(x)-\frac{1}{n}\log(\pi_{0}(x))$ has the following properties:

$x_{n}:=\operatornamewithlimits{argmin}I_{n}=0$ with $I_{n}(x_{n})=0$ and $I_{n}(x)=\frac{n+1}{n}x^{2}$ für $|x|\leq r_{0}$. 2. 2.

For $r_{n}<|x|\leq r_{n}^{+}$ we have $I_{n}(x)\equiv\delta$. 3. 3.

The quadratic approximation used in the Laplace approximation is then $\widetilde{I}_{n}(x)=\frac{n+1}{n}x^{2}$ for all $x$ which yields for $r_{n}<|x|\leq r_{n}^{+}$

[TABLE]

Note that Assumption 1 and items 1 and 2 of Assumption 2 are fulfilled. However, we will now show that the resulting sequence of measures $\mu_{n}$ with normalized Lebesgue density $\frac{1}{Z_{n}}\exp(-nI_{n})$ is not tight which by Prokhorov’s theorem precludes the existence of a weakly converging subsequence. In particular, the $\mu_{n}$ are not a Cauchy sequence in total variation or Hellinger distance. To this end, we first obtain for the normalizing constants $Z_{n}$ that

[TABLE]

Moreover, we have

[TABLE]

Hence, we get for the ball $B_{r_{0}}(x_{n})=\{x\colon|x|<r_{0}\}$

[TABLE]

This means, that for $n\to\infty$ the fixed ball $B_{r_{0}}(x_{n})$ around global peak of the density of $\mu_{n}$ has $\mu_{n}$-probability mass less than $\frac{1}{2}$. In other words, the measures $\mu_{n}$ do not concentrate around $x_{n}=0$. Additionally, we can assume that $\Phi_{n}$ is such that for a constant $c<\infty$ we have

[TABLE]

and, thus, $Z_{n}\leq 3+c$. But then, for all $n\in\mathbb{N}$

[TABLE]

which means that the sequence $\{\mu_{n}\}_{n}$ is not tight. This should demonstrate sufficiently why the third term in Assumption 2 is needed for a suitable concentration of the measures $\mu_{n}$.

C.1 A sufficient condition for the third term in Assumption 2

We provide now a useful sufficient condition on the functions $\Phi_{n}$ and $\pi_{0}$ such that the third term in Assumption 2 is fulfilled. In order to satisfy $I_{n}(x_{n})=0$ as required in Assumption 1 we set

[TABLE]

consider $\widehat{\Phi}_{n}(x):=\Phi_{n}(x)-\iota_{n}$ and represent the measures $\mu_{n}$ by

[TABLE]

with $I_{n}(x)=\widehat{\Phi}_{n}(x)-\frac{1}{n}\log\pi_{0}(x)$ in accordance to Remark 2. We then have

[TABLE]

By the first term in Assumption 2, we can conlude that $\pi_{0}(x_{n})\to\pi_{0}(x_{\star})>0$, i.e., the third term of Assumption 2 follows if there exists an $\epsilon>0$ such that

[TABLE]

and if

[TABLE]

Note that for $\Phi_{n}\equiv\Phi$ we have $\Phi(x_{n})-\inf_{\mathrm{S}_{0}}\Phi\in\mathcal{O}(n^{-2})$, since then $x_{\star}=\lim_{n\to\infty}x_{n}$ is a minimizer of $\Phi$ and $\Phi$ is strongly convex in a neighborhood of $x_{\star}$, cf. the paragraph “Preliminaries” in Section 3 and Remark 5.