Dimension free non-asymptotic bounds on the accuracy of high dimensional Laplace approximation

This note attempts to revisit the classical results on Laplace approximation in a modern non-asymptotic and dimension free form. Such an extension is motivated by applications to high dimensional statistical and optimization problems. The established results provide explicit non-asymptotic bounds on the quality of a Gaussian approximation of the posterior distribution in total variation distance in terms of the so called \emph{effective dimension} \( p_G \). This value is defined as interplay between information contained in the data and in the prior distribution. In the contrary to prominent Bernstein - von Mises results, the impact of the prior is not negligible and it allows to keep the effective dimension small or moderate even if the true parameter dimension is huge or infinite. We also address the issue of using a Gaussian approximation with inexact parameters with the focus on replacing the Maximum a Posteriori (MAP) value by the posterior mean and design the algorithm of Bayesian optimization based on Laplace iterations. The results are specified to the case of nonlinear inverse problem.


Introduction
High dimensional Laplace approximation has recently gained an increasing attention in connection with Bayesian inference for complicated nonlinear parametric models such as nonlinear inverse problems and Deep Neuronal Networks. Laplace approximation is obtained by replacing a log density with its second order Taylor approximation around the point of maximum. This leads to a Gaussian measure centered at the maximum with a covariance corresponding to the Hessian of the negative log-density (see, e.g., (Bishop, 2006, Section 4.4)). The asymptotic behavior of the parametric Laplace approximation in the small noise or large data limit has been studied extensively in the past (see, e.g., Wong (2001)). The asymptotic approximation of general integrals of the form e λf (x) g(x) dx by Laplace's method is presented in Olver (1974); Wong (2001). Non-asymptotic error bounds for the Laplace approximation can be found in Olver (1968) for the univariate case and in McClure and Wong (1983); Inglot and Majerski (2014) for the multivariate case. Lapiński (2019) studied the Laplace approximation error and its convergence in the limit λ → ∞ in the multivariate case when the function f depends on λ . Coefficients in the asymptotic expansion of the approximated integral are given in Nemes (2013).
Laplace approximation is an important step in establishing the prominent Bernsteinvon Mises (BvM) Theorem that quantifies the convergence of the scaled posterior distribution toward a Gaussian distribution in the large data or small noise limit. Parametric BvM theory is well-understood Van der Vaart (2000); Le Cam (2012). Modern applications with a high dimensional parameter space and limited sample size pose new questions and identify new issues in study of applicability and accuracy of Laplace approximation.
We refer to Lu (2017) for a study a parametric BvM theorem for nonlinear Bayesian inverse problems with an increasing number of parameters. A number of papers discuss the BvM phenomenon for nonlinear inverse problems; see e.g. Nickl (2020); Monard et al. (2019); Giordano and Kekkonen (2020), where the convergence is quantified in a distance that metrizes the weak convergence. Schillings et al. (2020) showed that the Laplace approximation error in Hellinger distance converges to zero in the order of the noise level.
The recent paper Helin and Kretschmann (2022) provides a finite sample error of Laplace approximation for the total variation distance with an explicit dependence on the dimen-sion and of the nonlinearity of the forward mapping. The Laplace approximation is also widely utilized for different purposes in computational Bayesian statistics; see e.g. Rue et al. (2009). 1.1 Motivation: Gaussian approximation of the posterior As one of the main motivation for this study, consider the problem of Bayesian inference for the log-likelihood function L(θ) = L(Y , θ) with data Y , a parameter θ ∈ Ê p and a Gaussian prior π ∼ N (θ 0 , G −2 ) . Here G −2 is a symmetric positive definite matrix in Ê p . Then the posterior density π G (·) of θ given Y is proportional to the product e L(θ) e − G(θ−θ 0 ) 2 /2 : where the sign ∝ means equality up to a normalizing multiplicative constant. Assume that the penalized maximum likelihood estimator (pMLE) θ G is well defined: Clearly θ G is also maximizer of π G (θ) . That is why it is often referred to as maximum a posteriori (MAP) estimator. Let also the log-likelihood function L(θ) be twice differentiable and weakly concave. Define Assuming the latter expression to be positive definite for all considered θ , define also its square root D G (θ) = G (θ) . We use the shortcut D G = D G ( θ G ) . Laplace's approximation means that the posterior distribution π G is close to the Gaussian distribution N ( θ G , D −2 G ) . A closely related Bernstein -von Mises phenomenon claims an approximation of the posterior by N ( θ, D −2 ) , where θ is the MLE and D 2 = = −∇ 2 L(θ * ) is the Fisher information matrix for the true parameter value θ * ; see e.g. Schillings et al. (2020) for a detailed discussion in context on nonlinear inverse problems. The mentioned results provide an efficient tool for Bayesian uncertainty quantification and constructing the elliptic credible sets as level sets of the approximating Gaussian distribution; see Helin and Kretschmann (2022) or Reich and Rozdeba (2020) for applications to drift and diffusion estimation.

Motivation: Gradient free optimization
Suppose that the point of maximum x * of the function ℓ(x) is not known and has to be evaluated numerically. We also assume that the function ℓ(·) is sufficiently smooth and can be efficiently computed at any point x , however, its gradient is not available. Nesterov and Spokoiny (2017) offered a powerful gradient free method based on averaging the exponent e ℓ(x) w.r.t. to some Gaussian distribution for x . Here we aim at reconsidering this idea within the Bayesian optimization framework. Namely, we intend to design a procedure based on Laplace iterations with an updated Gaussian prior at each step. Let π 0 be a starting Gaussian prior π 0 ∼ N (x 0 , G −2 0 ) . The corresponding posterior is defined by normalizing the product e ℓ(x) π 0 (x) . Due to the Laplace approximation result, this posterior is nearly Gaussian with the mean x 1 = argmax ℓ(x) − G 0 x 2 /2 ; see e.g. Theorem 2.1 below. This leads to the general idea of Bayesian optimization: use the posterior mean/variance to update the prior and repeat the Laplace approximation step. One may use standard Monte-Carlo or quasi Monte-Carlo methods for numerical approximation of the posterior mean; see e.g. Schillings et al. (2020) in context of Bayesian inverse problems.

Challenges
In spite of numerous existing results on Laplace approximation, some questions are still open. A small list of relevant issues is given below.
Prior impact in high dimension Modern applications in statistics and optimization force to extend the scope of applicability of the classical results on Laplace approximation by including the cases of a very large dimension even exceeding the sample size. It appears that in the case of a high-dimensional parameter, the prior is not washing out from the posterior. In the contrary to the classical theory, it becomes extremely important. It is of crucial importance to understand well the impact of the prior and its proper choice.
Smoothness conditions Standard smoothness conditions in terms of a uniform bound on the third derivative of f could be too restrictive and hard to check.
The use of posterior mean Another issue is possibility of using a Gaussian approximation of the posterior with inexact parameters. Indeed, the pMLE/MAP θ G is hard to compute, the classical MCMC based Bayesian computations deliver an estimate of the posterior mean and of posterior covariance. A big question is to justify the use of these quantities for uncertainty quantification or Bayesian optimization.
Non-concavity The classical results of Laplace approximation are stated under the condition that the function f (·) is strongly globally concave. In many applications including nonlinear inverse problems Nickl (2020); Schillings et al. (2020); Helin and Kretschmann (2022) or Deep Neuronal Networks such an assumption appears to be unrealistic.

This paper's contributions
This paper aims at reconsidering the classical results on Laplace approximation and to address the above issues. Below the list of the most important achievements in the paper.
Effective dimension and dimension free guarantees We introduce the notion of effective dimension p G of the problem which can be small of moderate even for huge parameter dimension p . The value p G is defined by an interplay between the information delivered by the data and information contained in the prior; see Section 2.1.2 for more details. Further we establish explicit non-asymptotic and dimension free guarantees for the accuracy of Gaussian approximation of the posterior in total variation (TV) distance in terms of effective dimension; see Theorem 2.1. In the case when the nonpenalized log-likelihood function grows linearly with the sample size n , the quality of Laplace approximation is of order p 3 G /n . It can be improved to p 3 G /n if instead of TV-distance, we limit ourselves to the class of centrally symmetric sets. The proofs combine classical variational arguments with sharp bounds for Gaussian quadratic forms.
Conditions require that f is strongly concave and locally smooth with a uniform bound on the third Gateaux derivative of f in a local vicinity of the point of maximum.
Critical dimension The result of Theorem 2.1 helps to address the issue of critical dimension for applicability of Laplace approximation: the relation p 3 G ≪ n between the sample size n and the effective dimension p G is sufficient for our main results. The result on concentration of the posterior only requires p G ≪ n .
Posterior mean in place of MAP The use of the posterior mean x in place of the MAP x * is justified by Theorem 2.10 under the same condition p 3 G ≪ n . The obtained results justify a Bayesian optimization procedure based on iterated Laplace approximation with a sequentially adjusted prior.
Non-concavity Section 3 indicates for the special case of nonlinear inverse problem how the assumption of strong concavity of f (·) can be relaxed using a "warm start" condition which means that a good starting guess is available.
The paper is organized as follows. Section 2 presents general bounds on the error of Laplace approximation and also discusses the use of posterior mean in place of posterior mode as well as the Bayesian optimization procedure based on Laplace iterations. Section 3 specifies the results to the case on nonlinear inverse problem. Proofs and technical tools are collected in the Appendix.

Dimension free bounds for Laplace approximation
Here we present some general results on accuracy of Laplace approximation.

Setup and conditions
Let f (x) be a function in a high-dimensional Euclidean space Ê p such that e f (x) dx = C < ∞ , where the integral sign without limits means the integral over the whole space Ê p . Then f determines a distribution È f with the density C −1 e f (x) . Let x * be a point of maximum: We also assume that f (·) is smooth, more precisely, three or even four time differentiable.
We aim at approximating the measure È f by a Gaussian measure N (x * , D −2 ) . Given a function g(·) , define its expectation w.r.t. È f after centering at x * : (2.1) A Gaussian approximation I D (g) for I(g) reads: The choice of the distance between È f and N (x * , D −2 ) depends on the considered class of functions g . The most strong total variation distance can be obtained as the supremum of |I(g) − I D (g)| over all measurable functions g(·) with |g(u)| ≤ 1 : The results can be substantially improved if only centrally symmetric functions g(·) with g(x) = g(−x) are considered. Obviously Moreover, as x * = argmax x f (x) , it holds ∇f (x * ) = 0 and Implicitly we assume that the negative Hessian D 2 = −∇ 2 f (x * ) is sufficiently large in the sense that the Gaussian measure N (0, D −2 ) concentrates on a local set U . This allows to use a local Taylor expansion for f ( is also strongly concave, then the mass of È f of the complement of U is exponentially small yielding the desirable Laplace approximation.
Motivated by applications to statistical inference, we consider f in a special form for a symmetric p -matrix G 2 ≥ 0 . Here ℓ(·) stands for the log-likelihood function while the quadratic penalty G(x − x 0 ) 2 /2 corresponds to a Gaussian prior N (x 0 , G −2 ) . We also assume that ℓ(·) is concave with D 2 def = −∇ 2 ℓ(x) ≥ 0 . Then clearly In typical asymptotic setups, the log-likelihood function ℓ(x) grows with the sample size or inverse noise variance, while the prior precision matrix G 2 is kept fixed. Decomposition (2.5) is of great importance for obtaining the dimension free results. The main reason is that the quadratic penalty does not affect smoothness properties of the function ℓ(·) but greatly improves the quadratic approximation term. Now we state precise conditions.
More specifically, we assume the following condition.
(C 0 ) There exists an operator G 2 ≤ −∇ 2 f (x * ) in Ê p such that the function If ℓ(·) in decomposition (2.5) is concave then this condition is obviously fulfilled.
More generally, if ℓ(·) in (2.5) is weakly concave, so that ℓ( plays an important role in our conditions and results. Remark 2.1. The condition of strong concavity of f on the whole space Ê p can be too restrictive; see the example of nonlinear inverse problem in Section 3. This condition can be replaced by its local version: there exists a set X 0 such that the Gaussian prior N (x 0 , G −2 ) concentrates on X 0 with a high probability and the maximizer x * G belongs to X 0 . In all the results, the integral over Ê p has to replaced by the integral over X 0 .

Effective dimension
With decomposition (2.6) in mind, we use another notation for D 2 = −∇ 2 f (x * ) : Also we write I G (g) instead of I D (g) in (2.2). The effective dimension p G is given by Of course, p G ≤ p but a proper choice of the penalty G 2 in (2.5) allows to avoid the "curse of dimensionality" issue and ensure a small or moderate effective dimension p G even for p large or infinite. The value p G helps to describe a local vicinity U around x * such that the most of mass of È f concentrates on U ; see Section A.4. Namely, let us fix some ν < 1 , e.g. ν = 2/3 , and some x > 0 ensuring that e −x is our significance level. Define (2.8)

Local smoothness conditions
Let p ≤ ∞ and let f (·) be a three times continuously differentiable function on Ê p . We fix a reference point x and local region around x given by the local set U ⊂ Ê p from (2.8). Also consider the second order Taylor approximation , u ⊗2 and similarly the third order expansion and introduce the remainders The use of the Taylor formula allows to bound It is worth noting that the quadratic penalty and the quadratic term in definition of the values δ k (x, u) cancels, k ≥ 3 . Local smoothness of f (·) or, equivalently, of ℓ(·) , at x will be measured by the value ω(x) : We also denote Our results apply under the condition ω ≪ 1 . Local concentration of the measure È f requires ω ≤ 1/3 ; see Proposition A.11. The main results about Gaussian approximation of È f are valid under a stronger condition ω p G ≤ 2/3 with the effective dimension p G from (2.7).

Error bounds for Laplace approximation
Our first result describes the quality of approximation of the measure È f by the Gaussian measure N (x * , D −2 G ) with mean x and the covariance D −2 G in total variation distance. In all our result, the value x is fixed to ensure that e −x is negligible. First we present the general results which will be specified later under the self-concordance condition.

General bounds
Theorem 2.1. Suppose (C 0 ) . Let also p G be defined by (2.7) and r G and U by (2.8).
Introduce the following conditions.

Corollary 2.3. Under the conditions of Theorem 2.2, it holds for
.

Bounds under self-concordance
Often the smoothness properties of f are described in terms of the third derivative.
Our approach does not require a bounded third or fourth full dimensional derivative, instead we consider directional Gateaux derivatives. Let ℓ(·) be three times continuously differentiable and ∇ 3 ℓ(·) stand for the third derivative. Then the Taylor expansion of the third order implies ω ≤ α/3 with Moreover, in many applications, the function ℓ(·) is of the form −ℓ(x) = nh(x) for a fixed smooth function h(·) satisfying the following condition: for a strongly convex function h(·) and (S 4 ) −ℓ(·) = nh(·) for a strongly convex function h(·) and This is a local version of the so called self-concordance condition; see Nesterov (1988).
Under this condition, we can easily bound the values δ k (x * , u) for k = 3, 4 , ω and α : Finally we state the bounds under (S 3 ) and (S 4 ) .

Critical dimension
Here we briefly discuss the important issue of critical dimension. Theorem 2.1 states concentration of È f under the condition ω ≤ 1/3 and Gaussian approximation under the stronger condition ω p G ≤ 2/3 . Moreover, self-concordance (S 3 ) and Lemma A.12 enables us to bound ω p G /n and ω p G p 3 G /n . Hence, concentration of È f requires the critical dimension condition p G ≪ n , while Gaussian approximation applies under p 3 G ≪ n . We see that there is a gap between these conditions. We guess that in the region n 1/3 p G n , a non-Gaussian approximation of the posterior is possible.
Bochkina and Green (2014) provides examples of non-Gaussian posterior limits for nonregular models.

Kullback-Leibler divergence
Theorem 2.1 through 2.4 quantify the approximation È f ≈ N (x * , D −2 G ) in the total variation distance. Another useful characteristic could be the Kullback-Leibler (KL) divergence between È f and N (x * , D −2 G ) . The KL divergence K (È 1 , È 2 ) = 1 log(dÈ 1 /dÈ 2 ) is asymmetric, K (È 1 , È 2 ) = K (È 2 , È 1 ) with few exceptions like the case of Gaussian measures È 1 and È 2 . Moreover, K (È 1 , È 2 ) can explode if È 2 is not absolutely continuous w.r.t. È 1 . We present two bounds for each ordering. For ease of presentation, we limit ourselves to the case when either (L 3 ) or (S 3 ) meets. (2.18) Now we briefly discuss the value K (È G , È f ) . We already know that È f concentrates on U and can be well approximated by È G on U . However, this does not guarantee a small value of K (È G , È f ) . It can even explode if e.g. È f has a compact support. In fact, the log-density of

Mean and MAP
Here we present the bound on |I(g) − I G (g)| for the case of a linear vector function (2.20) Now we specify the result for the special choice Q = D .
2. An interesting question is whether the result of Theorem 2.7 or Corollary 2.8 applies with Q = D G . This issue is important in connection to inexact Laplace approximation; see the next section. The answer is negative. The problem is related to the last term C e −x in the right hand-side of (2.20). The constant C involves the moments of Q(X − x * ) 2 which explode for Q = D G and p = ∞ .

Inexact approximation and the use of posterior mean
Now we change the setup. Namely, we suppose that the true maximizer x * of the function f is not available, but x is somehow close to the point of maximum x * .
Similarly, the negative Hessian D 2 G (x * ) = −∇ 2 f (x * ) is hard to obtain and we use a proxy H 2 . We already know that È f can be well approximated by N (x * , D −2 ) with . This section addresses the question whether N (x, H −2 ) can be used instead. Here we may greatly benefit from the fact that Theorem 2.1 provides a bound on the total variation distance between È f and N (x * , D −2 Therefore, it suffices to bound the TV-distance between two Gaussian distributions: N (x * , D −2 G ) naturally appears in the Laplace approximation, the second one is a proxy used instead. Pinsker's inequality provides an upper bound: for any two measures P, Q where K (P, Q) is the Kullback-Leibler divergence between P and Q . The KL-divergence between two Gaussians has a closed form: However, Pinsker's inequality is only a general upper bound which applied to any two distributions P and Q . If P and Q are Gaussian, it might be too rough. Particularly the use of tr B 2 is disappointing, this quantity is full dimensional even if each of D −2 G and H −2 has a bounded trace. Also dependence on to the class B el (Ê p ) of elliptic sets A in Ê p of the form A = u ∈ Ê p : Q(u − x) ≤ r for some linear mapping Q : Ê p → Ê q , x ∈ Ê p , and r > 0 . Given two symmetric q -matrices Σ 1 , Σ 2 and a vector a ∈ Ê q , define where λ 1 is the vector of eigenvalues of Σ 1 arranged in the non-increasing order and similarly for λ 2 . Götze et al. (2019) stated the following bound for γ 1 ∼ N (0, Σ 1 ) and Fr /3 , k = 1, 2 . By the Weilandt-Hoffman inequality, λ 1 − λ 2 1 ≤ Σ 1 − Σ 2 1 , see e.g. Markus (1964). Here M 1 = tr |M | = j |λ j (M )| for a symmetric matrix M with eigenvalues λ j (M ) .
Furthermore, for X ∼ È f and γ ∼ N (0, I I p ) , any linear mapping Q : As a special case, we consider the use of the posterior mean Theorem 2.10. Assume the conditions of Theorem 2.7 and Theorem 2.9. Then it holds for any linear mapping Q : The same bound holds with Q = n 1/2 I I p in place of D provided that D 2 ≥ C 0 n I I p for some fixed C 0 > 0 . We may conclude that the use of posterior mean x in place of the posterior mode x * is justified under the same condition on critical dimension p 3 G ≪ n as required for the main result about Gaussian approximation.

Bayesian optimization and iterated Laplace approximation
Suppose that the point of maximum x * of the function ℓ(x) is not known and has to be evaluated numerically. We also assume that the function ℓ(·) is sufficiently smooth and can be efficiently computed at any point x , however, its gradient is not available.
Nesterov and Spokoiny (2017) offered a powerful gradient free method based on averaging the exponent e ℓ(x) w.r.t. to some Gaussian distribution for x . Inspired by the results of previous sections, we propose here a modified version of Nesterov and Spokoiny (2017) based on iterated Laplace approximations. Starting from some π 0 = N (x 0 , G −2 0 ) , we iteratively update the Gaussian prior π k ∼ N (x k , G −2 k ) and try to numerically assess the corresponding posterior obtained by normalization of e f k (x) for f k (x) = ℓ(x) − G k (x − x k ) 2 /2 . This posterior is known to be nearly Gaussian by Theorem 2.1, its mean x k+1 can be efficiently estimated by Monte-Carlo or quasi Monte-Carlo sampling; see e.g. Schillings et al. (2020). We use x k+1 for building the Gaussian prior for the step k + 1 . Moreover, due to our results, the value x k+1 is nearly the maximizer of the corresponding quadratic approximation of f k (x) in the vicinity of x k . One can say that computing the posterior mean x k+1 mimics well a step of the second order Newton-Raphson optimization method but does not require computing the gradient and the inverse Hessian matrix. The prior precision matrix G 2 k can be taken in the form λ −1 k G 2 0 , where λ k replaces the step size, λ k → 0 as k increases. The procedure can be written as follows.
The starting guess is important as well the choice of the step multiplier a . Moreover, one can use the variable multiplier a k and incorporate the posterior covariance for variance reduction schemes. One can make a rigorous analysis of the convergence of this algorithm similarly to Nesterov and Spokoiny (2017).
The proposed procedure is closely related to ensemble Kalman filtering technique; see e.g. Schillings and Stuart (2017), Reich (2022) and references therein.
Algorithm 1 Laplace iterations 1. Start with k = 0 ; Fix x 0 and G 2 0 ; 2. Draw a sample (x (m) ) m≤M from N (x k , G −2 k ) . For each m ≤ M , compute w (m) = e ℓ(x (m) ) and update the sums W k = m w (m) and S k = m w (m) x (m) .
3. Compute the value x k+1 = S k /W k as the next prior mean; update the posterior precision matrix G 2 k+1 = aG 2 k . 4. Increase k by one and repeat from Step 2. 5. Iterate until convergence.

Laplace approximation for non-linear inverse problem
Let m(x) = m i (x), i ≤ n ∈ Ê n be a nonlinear mapping (operator) of the source signal x ∈ Ê p to the target space Ê n . We consider the problem of inverting the relation z = m(x) : given an image vector z ∈ Ê n , recover the corresponding source x ∈ Ê p . This leads to the nonlinear least square problem of maximizing the function ℓ(·) of the form Given a Gaussian prior N (x 0 , G −2 ) , we consider Laplace's approximation for the penalized function Laplace's approximation requires weak concavity of ℓ(x) ; see However, such a condition seems can be hard to ensure for all x ∈ Ê p unless m(·) is linear. Instead we restrict ourselves to some subset X 0 ⊂ Ê p which can be viewed as a concentration set of the Gaussian measure N (x 0 , G −2 ) . As the prior mass of the complement of X 0 is exponentially small, when sampling from the prior N (x 0 , G −2 ) , one would need exponentially many samples to hit any set A outside of X 0 . Therefore, we restrict the prior to the set X 0 by skipping all draws of x with x ∈ X 0 . In what follows we also assume that x 0 is a reasonable guess ensuring that the target x * G belongs to X 0 . Then the result about Laplace approximation applies even after restricting the parameter space to X 0 . This local set X 0 is defined as elliptic concentration set for the Gaussian prior N (x 0 , G −2 ) in the form where Q is a linear operator in Ê p and r 0 is fixed to ensure that the most of prior mass is within X 0 . Theorem B.4 suggests X 0 in form of a ball around x 0 with the radius of order √ tr G −2 . A proper choice of G 2 ≥ I I p has to ensure a small value of tr(G −2 ) . Note that the operator Q can be scaled by any factor. Then the radius r 0 in (3.5) will be scaled correspondingly.
Introduce for any x ∈ X 0 a p -symmetric matrix This is the first term in expansion (3.3) corresponding to a linear approximation of m(·) at x . Injectivity of m(·) means thatD 2 (x) is positive definite and well conditioned. If n −1D2 (x 0 ) ≤ G 2 then the choice Q =D(x 0 ) in (3.4) is our alternative to Q = I I p .
Below we assume the following regularity condition.
(∇m) Let the local set X 0 be defined by (3.4) and (3.5) with Q =D 0 def =D(x 0 ) . For some fixed C 2 , C n , and any x ∈ X 0 , u ∈ Ê p , it holds

7)
and with some C 0 ≥ 1D Condition (∇m) only requires some local regularity of the m i (x) 's and injectivity of the gradient mapping ∇m(x) within X 0 . The shape of X 0 is defined by the prior covariance G −2 which has to be selected to ensure that the set X 0 is a local set. More precisely, we connect the size of X 0 and local regularity of m(·) by the following condition.
(r 0 ) With X 0 defined by (3.4) and (3.5) and C n , C 0 , and C 2 from (∇m) , it holds To get a feeling of this condition, consider the case withD 2 = nI I p . Then (3.5) yields r 2 0 ≍ n tr G −2 and (3.9) requires tr G −2 ≪ 1 , that is, the prior N (x 0 , G 2 ) is supported to a small vicinity of x 0 . Now we introduce the key condition of "warm start" which informally means a good choice of the starting point x 0 ensuring that the local vicinity X 0 of x 0 contains a point Later we show that (∇m) , (r 0 ) , and (x 0 ) ensure for some δ ≤ 3/4 , and hence, ℓ(x) is strongly concave for x ∈ X 0 . Now we consider the point x * G from (3.2) and its local vicinity. It is important to secure that x * G is within X 0 . This appears to be automatically fulfilled if the true solution x * = m −1 (z) is within X 0 ; see Lemma A.14. Now we are back to the setup with a strongly concave function ℓ(x) for x ∈ X 0 . Similarly to the general case, define The effective dimension is given by p G = tr(D 2 D −2 G ) ; see (2.7). The local vicinity U G of x * G can be defined by the rule (2.8). Alternatively one can useD 2 =D 2 (x * G ) instead of D 2 : The radius r has to be adjusted to ensure that the such defined vicinity of x * is not smaller than U from (2.8): U G ⊆ U . Relation (3.10) suggests using r = 2r G .
For applying the general results of Theorems 2.1 through 2.4, the function ℓ(x) = − m(x)−z 2 /2 has to be sufficiently smooth. We now present some sufficient conditions in terms of the functions m i (x) . These conditions extend (∇m) .
(∇ 4 m) For some C 4 , it holds uniformly over x ∈ U and u ∈ Ê p We are now prepared to state the main result about Laplace's approximation which is a straightforward application of Theorem 2.1 to the measure È f restricted to X 0 .

A Tools and proofs
Here we collect the proofs of the main results and some useful technical statements about the error of Laplace approximation. Below we write x instead of x * . Note that after passing to representation (2.3), many results below apply to any x , not necessarily for x = x * . We only use D 2 G = −∇ 2 f (x) and ω instead of ω(x) . Everywhere we assume the local set U to be fixed by (2.8). We separately study the integrals over U and over its complement. The local error of approximation is measured by As a special case with g(u) ≡ 1 we obtain an approximation of the denominator in (2.3).
In addition, we have to bound the tail integrals The analysis will be split into several steps.

A.1 Overall error of Laplace approximation
First we show how to seam together the error ♦ of local approximation and the bounds for the tail integrals ρ and ρ G ; see (A.2).
Proposition A.1. Suppose that for a function g(u) with g(u) ∈ [0, 1] and some ♦, ♦ g Then with ρ and ρ G from (A.2) Similarly for g(u) ≥ 0 Putting together all these bounds yields (A.5).
The next corollary is straightforward.

A.2 Lower and upper Gaussian measures
This section introduces the lower and upper Gaussian measure which locally sandwich the measure È f using the decomposition from condition (C 0 ) . Denote −∇ 2 f (x) = D 2 G . Definition (2.9) enables us to bound with ω = ω(x) yielding two Gaussian measures which bounds È f locally from above and from below.
The next technical result provides sufficient conditions for their contiguity.

A.3 Local approximation
This section presents the bounds on the error ♦ of local approximation (A.1). The first result only uses ω p G < 1 while the second one also assumes (L 3 ) . These two results allow to bound the total variation distance between È f and N (x, D −2 G ) . The third result is finer and is based on (L 4 ) . It allows to improve the error term of Gaussian approximation over the class of centrally symmetric sets. We also present some extensions for the moments of È f . Proposition A.5. Let ω = ω(x) from (2.9) and p G from (2.7) satisfy Then for any function g(u) with |g(u)| ≤ 1 Moreover, under (L 3 ) , the bound applies with ♦ = ♦ 3 ; see (2.14).

In view of
for a standard normal γ , we derive yielding the second statement.
The result can be extended to the case with a homogeneous function g(u) , e.g.
Proposition A.6. Suppose the conditions of Proposition A.5 and (L 3 ) . Then for m ≥ 1 and any m -homogeneous function g(·) with g(tu) = t m g(u) Proof. Similarly to the proof of Proposition A.5, under (2.12) for R(t) from (A.14) This yields (A.17).
Important special cases correspond to m = 1 .
Proposition A.7. Suppose the conditions of Proposition A.6. Then with γ G ∼ N (0, D −2 G ) , it holds for any linear mapping Q : Ê p → Ê q and any vector a ∈ Ê q Now we state a sharper result based on (L 4 ) .
Proposition A.8. Suppose the conditions of Proposition A.5 and (L 4 ) . Then for any function g(u) with |g(u)| ≤ 1 and g(u) = g(−u) where for a Gaussian element γ G ∼ N (0, D −2 If the function g(·) is not bounded by one but it is symmetric and 2m -homogeneous, i.e. g(tu) = t 2m g(u) , then (A.19) still applies with Proof. Below we write f (3) instead of ∇ 3 f (x) and δ k (u) in place of δ k (x, u) , k = 3, 4 .

A.4 Tail integrals
In this section we also write x in place of x * . Below we evaluate ρ from (A.2) which bounds the integral of e f (x;u) over the complement of the local set U of a special form U = u : Du ≤ ν −1 r G for D from (C 0 ) . Our results help to understand how the radius r G should be fixed to ensure ρ sufficiently small. The main tools for the analysis are deviation probability bounds for Gaussian quadratic forms; see Section B.1.

This yields by r
Integration by parts allows to represent the last integral as This completes the proof of the result (A.22). The second statement (A.23) is about Gaussian probability È T γ ≥ r G for a standard normal element γ , and we derive and (A.23) follows.
The next result extends (A.22).
Proposition A.10. Assume the conditions of Proposition A.9 with for some m ≥ 0 . Then (A.22) can be extended to Proof. The case m > 0 can be proved similarly to m = 0 using m log z ≤ mz .

A.5 Local concentration
Here we show that the measure È f well concentrates on the local set U from (2.8).
Again we fix x = x * .
Proposition A.11. Assume ω ≤ 1/3 . Then Moreover,  A.7 Proof of Theorem 2.4 We make use of the following lemma.
A.8 Proof of Theorem 2.5 and Theorem 2.6 Define C G def = log e − D G u 2 /2 du − log e f (x * ;u) du .
For u = x − x * , it holds as in (2.3) Similarly to (A.15) and (A.16), we can bound On the complement u ∈ U c , we use that The last inequality is based on concavity of ℓ(·) and local approximation ℓ(x * ; u) ≈ − Du 2 /2 within U yielding ℓ(x * ; u) < 0 for u ∈ U c . By Proposition A.10 with We conclude that Similarly This yields e C G ≤ 1 + ♦ 3 + 2e −x and C G ≤ ♦ 3 + 2e −x . Putting all together results in A.10 Proof of Theorem 3.1 The basic idea of the proof is to reduce the statements of the theorem to the case of Theorem 2.1. Conceptually, the most important step of the proof is to show that the function f (x) is strongly concave on the subset X 0 which is a concentration set of Lemma A.13. Suppose (∇m) , (r 0 ) , and (x 0 ) . Then for all x ∈ X 0 Given x ∈ X 0 , denote u = x − x * . The definition of X 0 implies D 0 u ≤ 2r 0 . The first order Taylor expansion yields Hence, by (3.7) and (3.8) The next step is to show that under (x 0 ) , the point x * G is also within X 0 .
Bound (A.31) means that (2.16) is fulfilled with c 3 = 4C 3/2 G C 3 C n . This implies (3.12) similarly to Theorem 2.4. Moreover, (A.32) implies (S 4 ) with c 4 = 8C 2 G C 4 C n , and the second statement of the theorem follows as well.
Now we consider the case of a non-centered quadratic form Bγ, γ /2 + A, γ for a fixed vector A .

B.2 Deviation bounds for Gaussian quadratic forms
The next result explains the concentration effect of Bξ, ξ for a centered Gaussian vector ξ ∼ N (0, V 2 ) and a symmetric trace operator B in Ê p , p ≤ ∞ . We use a version from Laurent and Massart (2000). For completeness, we present a simple proof of the upper bound.