Laplace's method in Bayesian inverse problems
Philipp Wacker

A charged anisotropic well-behaved Adler-Finch-Skea solution Satisfying Karmarkar Condition
    Sumita Banerjee 
    Jan 15 2017
===========================================================================================


In a Bayesian inverse problem setting, the solution consists of a posterior measure obtained by combining prior belief (where we restrict ourselves to Gaussian priors), information about the forward operator, and noisy observational data. This posterior measure is most often given in terms of a density with respect to a reference measure in a high-dimensional (or infinite-dimensional) Banach space. Although Monte Carlo sampling methods provide a way of querying the posterior, the necessity of evaluating the forward operator many times (which will often be a costly PDE solver) prohibits this in practice. For this reason, many practitioners choose a suitable Gaussian approximation of the posterior measure, by a procedure called Laplace's method. Once generated, this Gaussian measure is a easy to sample from and properties like moments are immediately acquired. This paper derives Laplace's approximation of the posterior measure attributed to the inverse problem explicitly by a second-order approximation of the data-misfit functional, specifically and rigorously in the infinite-dimensional setting. By use of a reverse Cauchy-Schwarz inequality we are able to explicitly bound the Hellinger distance between the posterior and its Laplace approximation.  



§ INTRODUCTION

We consider a inverse problem 

    y = G(u) + η

where G:X→ Y is a (possibly) nonlinear mapping between Hilbert spaces X, Y and η∼ N(0, Γ) is additive noise. The challenge consists of inferring the value of u from the noisy (and usually lower-dimensional) observation y. This is an ill-defined problem in general, so some sort of regularization is needed.

In a Bayesian approach (see <cit.>) under Gaussian assumptions, we assume that u∼μ_0 = N(0,C_0), i.e. we have a Gaussian prior on the variable u. For simplicity we assume that the mean is 0, but this assumption can be dropped with slight modifications. The prior acts as a regularization and makes the inverse problem well-defined: Standard theory yields the posterior measure μ on u given an observation y under mild assumptions on the forward operator G:

    μ/μ_0(u) = exp(-Φ(u))/∫exp(-Φ(u)) μ_0(u),

where Φ(u) = y-G(u)_Γ^2/2.
Especially in higher dimensions, the posterior is in practice often approximated by a suitable Gaussian, in order to make computation of moments (or frequentist confidence sets) feasible. For a beautiful example how this is done in practice, see <cit.> in the setting of optimal experimental design of an infinite-dimensional inverse problem, or in a finite-dimensional context in <cit.>. This procedure is called Laplace's method and is the focus of this work. We define the functional I(u) = Φ(u) + u_C_0^2/2 which can also be thought of as a (Tikhonov-)regularized cost functional in the classical sense. The maximum a posteriori point is

    := _u I(u)

and the Laplace approximation is defined as

    ν = N(, HI()^-1)

This means that the Laplace approximation is the Gaussian measure centered at the MAP point, with covariance operator matching the “local” covariance structure of the posterior measure. In finite dimensions, its density is exactly the normalized exponential of the measure's negative lognormal local quadratic approximation. For a good explanation (and as a general recommendation for a truly enjoyable book) of Laplace's method in finite dimensions, see <cit.>. The paper <cit.> gives a numerical analysis of the Laplace approximation as well as the so-called Bayesian information criterion; <cit.> and the references therein give an overview over the use of Laplace's method in the Machine Learning and imaging community. A treatise about approximation of measures by Gaussian measures can be found in <cit.>, although they employ the Kullback-Leibler divergence (or relative entropy) as a notion of distance between measures. <cit.> is an extensive study of various Gaussian approximation methods (including Laplace's method) in the context of reservoir modelling. Standard results about Laplace's method are recorded in <cit.> and <cit.>, with newer results on Laplace's method in <cit.>. Gaussian approximations in a different context have been treated in <cit.>, <cit.> in the case of diffusion processes and in <cit.> in the context of molecular dynamics. <cit.> presents an approach running “anti-parallely” to Laplace's method: Instead of approximating the posterior measure directly, they approximate the forward operator or the negative log-likelihood. They bound the Hellinger distance between the true posterior and the resulting approximation, which is also the method of attack in this paper. <cit.> discusses possible algorithms for implementing Laplace's method and similar approximation procedures in Bayesian inference problems. The MAP point is here defined as (<ref>), i.e. a minimization point of the functional I. See <cit.> for a characteristation of the MAP point in a broader setting, by dropping the Gaussianity assumption (in particular, they derive an expression for the MAP point in a linear inverse problem with Besov priors). Original accounts of probabilistic methods of tackling inverse problems can be found in e.g. <cit.> for a random processes view in a Hilbert space setting, <cit.> for – among other things – an explanation of how classical regularization can be viewed as application of a Bayesian prior. The difficulty of constructing well-defined linear estimators is addressed in <cit.> and <cit.>. In <cit.>, a convergence result of the Ky-Fan metric (which is another metric between measures) between posterior measure and the minimum-norm-least-squares point (which is related to the maximum likelihood point) is proven. Although inverse problems in practice will always be finite-dimensional, it has proven worthwhile to study inverse problems from an infinite-dimensional viewpoint. This stems from the fact that although many algorithms work nicely for a given finite data set but break in the limit of finer and finer spatial or temporal resolution. One notable example is total variation regularisation, which is found to be edge-preserving, although it loses this property in the infinite-dimensional limit, see <cit.>. Another intuitive point is that algorithms that work in the infinite-dimensional setting should be complexity-invariant to a refinement of spatial resolution (naive algorithms often have high polynomial complexity with respect to the data). For a discussion of this fact see <cit.>. For more details on how to construct well-defined infinite-dimensional inverse problems, see <cit.> and <cit.>.

It can be shown that the Laplace approximation ν coincides with the posterior measure μ if the forward operator G is linear and the prior measure was Gaussian in the first place. Heuristically, the approximation is bad when the posterior measure is multimodal or has different tail properties than a Gaussian. 



We are interested in deriving concrete error bounds for the approximation quality μ≈ν. The Hellinger distance between probability measures lends itself to this cause. Given two measures μ,ν which are absolutely continuous w.r.t. another measure μ_0, the Hellinger distance (which is independent of the choice of μ_0) between μ and ν is
 
    (d_H(μ,ν))^2 = 1/2·∫(√(%s/%s)μμ_0- √(%s/%s)νμ_0)^2μ_0 = 1 - ∫√(%s/%s)μμ_0√(%s/%s)νμ_0μ_0


 The main conclusion of this paper is recorded informally in the following claims; they are correctly stated and proven later.
 
 While the posterior measure is given via its density w.r.t. the prior by
 
    μ/μ_0(u) = exp(-Φ(u))/∫exp(-Φ(u)) μ_0(u),

Laplace's method yields a Gaussian approximation of the form ν = N(, HI()^-1) and its density w.r.t. the prior is

    ν/μ_0(u) = exp(-TΦ(u))/∫exp(-TΦ(u)) μ_0(u).

Here, TΦ(u) = Φ() + DΦ()(u-) + 1/2HΦ()[u-, u-], the second order Taylor approximation of the data-misfit functional Φ at . 
 
   
If there is K∈(0,1) such that

    exp(-Φ(u)/2) - exp(-TΦ(u)/2)_L^2(X, μ_0)≤ K·exp(-I())/√((C_0^1/2· HΦ()· C_0^1/2 + Id)),

then the Hellinger distance between the posterior and its Laplace's approximation can be bounded:

    d_H(μ,ν) ≤K/√(1 + (1-K)^2)

 
 While approximation results for Laplace's method have been achieved before (e.g. in the references made above), it is to the best of our knowledge that the claims as stated have not been formulated in an Inverse problem setting where Laplace's method is often used.
 
 The paper is organized as follows: First the more intuitive one-dimensional case is presented and representation (<ref>) is derived. The approach taken can not be generalized to the infinite dimensional case, as we use densities with respect to a Lebesgue measure. Then the infinite-dimensional equivalent is proven, where we show the equality (<ref>) directly by means of characteristic functions. The following section shows that the problem of bounding the Hellinger distance can be reduced to a reverse Cauchy-Schwarz inequality. After recording a few elementary results about reverse CS inequalities, we immediately obtain claim <ref> and also a more practical (but less tight) version of it.
 
 
 A caveat regarding practicality: Although proposition <ref> and corollary <ref> yield expressions that seem to offer readily-checked conditions for explicit bounds on the Hellinger distance, a non-artificially constructed example where claim 2 easily yields a non-trivial bound on d_H seems hard do obtain and the applicability of the second claim should thus not be overstated. This is explained in more detail in the section about application. Rather, the main point of this paper is the derivation and proof of characterization of expression (<ref>).
 


§ THE LAPLACE APPROXIMATION IN ONE DIMENSION

We recall 

    μ/μ_0(u) = exp(-Φ(u))/∫exp(-Φ(u)) μ_0(u),

and with for the Lebesgue measure λ in one dimension, the Laplace approximation of μ, which is a Gaussian centered at the maximum a posteriori point, with variance equal to the inverse  of the second derivative of I at this point. In finite dimensions, any measure is most easily recorded by its density w.r.t. λ.

    ν/λ(u) = √(%s/%s)I”()2πexp( -I”()/2· (u-)^2),

hence

    ν/μ_0(u) = ν/λ(u)/μ_0/λ(u) = √(I”())·σ·exp( -I”()/2· (u-)^2 + u^2/2σ^2).

As we mentioned, second-order Taylor approximations will play a leading role, hence we define 

    R(u) = I(u) - I() - I”()/2(u-)^2,

the error term of the second order Taylor approximation of I in . Note that the first-order term vanishes because of  being a minimum of I, which we assume to be C^2. Interestingly, R(u) is also an error term for the second order approximation of Φ, albeit in a slightly different way:

    R(u)    = Φ(u) + u^2/2σ^2 - Φ() - ^2/2σ^2-Φ”()/2(u-)^2 - 1/σ^2(u-)^2
       = · (u-)/σ^2 + Φ'()(u-) 
       +[ Φ(u) - Φ() - Φ'()(u-) - Φ”()/2(u-)^2]
       = Φ(u) - T^(2)_Φ(u).

Note that the terms in the second line vanish because of 0 = I'() = Φ'() + /σ^2 and the second term is exactly the error between Φ(u) and its second order Taylor polynom developed in . It holds that I'() = 0 but in general Φ'()≠ 0.
With this definition of R, and especially

    -I”()/2· (u-)^2 + u^2/2σ^2 = -Φ(u) + I() + R(u),

we obtain for the density of ν w.r.t. μ_0:

    ν/μ_0(u) = √(I”())·σ·exp(I()) ·exp(-Φ(u) + R(u)).

An easy calculation shows

    exp(-I())/σ√(I”())   = ∫exp(u^2/2σ^2-I()-I”()/2(u-)^2)μ_0( u)
       = ∫exp(-Φ(u)+R(u))μ_0( u) = ∫exp(-T^(2)_Φ(u))μ_0(u).

This also follows immediately from the normalization of ν/μ_0 and thus

    ν/μ_0(u) = exp(-Φ(u)+R(u))/∫exp(-Φ(u)+R(u))μ_0( u) = exp(-TΦ(u))/∫exp(-TΦ(u))μ_0( u).





















































































§ THE LAPLACE APPROXIMATION IN GENERAL HILBERT SPACES

In this section, we will show claim <ref> in the setting of a general (possibly infinite-dimensional) Hilbert space. In order to do this, we need a slight adaptation of a useful Gaussian integral calculation.

In Proposition 1.2.8 in <cit.>, the authors prove the following: Let μ_0 = N(0, Q) be a Gaussian measure on a real Hilbert space H. Assume that M is a symmetric operator such that ⟨ Q^1/2MQ^1/2 u, u⟩ < ⟨ u,u⟩ for all 0≠ u∈ H. Then for b∈ H

    ∫_H exp(1/2⟨ M u,u⟩ + ⟨ b, u⟩) μ_0(u) = exp(1/2|(1-Q^1/2MQ^1/2)^-1/2· Q^1/2b|^2 )/√( (1-Q^1/2MQ^1/2)).

We will need a generalization of this formula, which follows from analytical continuation of the (real) Hilbert space's inner product to its complex extension. Recall that this continuation will not be positively definite anymore: ⟨λ_1 a + b, λ_2 c + d⟩ = λ_1 λ_2 ⟨ a, c⟩ + λ_2 ⟨ b, c⟩ + λ_1⟨ a, d⟩ + ⟨ b,d⟩, without any complex conjugations on any λ_i. The following lemma is stated without proof as it follows immediately from the bilinearity of the analytical continuation stated above.


Let μ_0 = N(0, Q) be a Gaussian measure on a real Hilbert space H. Assume that M is a symmetric operator such that ⟨ Q^1/2MQ^1/2 u, u⟩ < ⟨ u,u⟩ for all 0≠ u∈ H. Then, with L:=Q^1/2(1-Q^1/2MQ^1/2)^-1Q^1/2 and for b_1, b_2∈ H

    ∫_H    exp(1/2⟨ M u,u⟩ + ⟨ b_1+ib_2, u⟩) μ_0(u) 
       = exp(1/2⟨ Lb_1,b_2⟩ + i⟨ Lb_1, b_2⟩ - 1/2⟨ Lb_2, b_2⟩)/√( (1-Q^1/2MQ^1/2)).



Now we can prove our main result: 

Consider the inverse problem <ref> with prior μ_0 and posterior μ given by 
    μ/μ_0(u) = exp(-Φ(u))/∫exp(-Φ(u))μ_0.
 
The functional I(u) = Φ(u) + 1/2u_C_0^2 is assumed to be C^2 in a neighborhood of =  I(u).
Then the Laplace approximation of μ given by

    ν = N(, I̋()^-1)

is equivalently defined by

    ν/μ_0(u) = exp(-TΦ(u))/∫exp(-TΦ(u))μ_0,
 
where TΦ(u) = Φ() + DΦ()(u-) + 1/2Φ̋()[u-, u-] is the second order Taylor approximation of Φ generated in 




This is done by comparing the Fourier transform (or characteristic function) of both representations for ν. We calculate


    ∫exp(-TΦ(u))μ_0(u) = e^-I()/√( C^1/2I̋()C^1/2)

and

    ∫exp(i⟨λ, u⟩ -TΦ(u))μ_0(u) = e^-I()·exp( i⟨, λ⟩ - 1/2·I̋()^-1[λ, λ])/√( C^1/2I̋()C^1/2)

We show (<ref>), as (<ref>) follows from setting λ = 0. We use 

    -TΦ(u) = R(u) - Φ(u) = u_C^2/2 - I() - 1/2I̋()[u-, u-].

and write J = I̋() and v= for brevity. Note that for a bilinear operator K we will identify K(w, z) = ⟨ K w, z⟩. Then 
    ∫exp(i⟨λ, u⟩   -TΦ(u))μ_0(u)
       = ∫exp(i⟨λ, u ⟩  + u_C^2/2 - I(v) - 1/2J[u-v, u-v] )μ_0(u)
       =e^-I(v)·∫exp( ⟨ iλ, u⟩ + 1/2⟨ C^-1 u,u⟩  -1/2⟨ J(u-v),u-v⟩)μ_0(u)
       = e^-I(v) - 1/2⟨ Jv,v⟩·∫exp(⟨ Jv + iλ, u⟩ + 1/2⟨ (C^-1-J)u,u⟩)μ_0(u).
    This is formula (<ref>) with M = C^-1-J and b_1 = Jv, b_2=λ. In this case, 1 - C^1/2MC^1/2 = 1 - C^1/2(C^-1-J)C^1/2 = C^1/2JC^1/2 and thus (1 - C^1/2MC^1/2)^-1/2 = J^-1/2C^-1/2. Continuing,   = e^-I(v) - 1/2⟨ Jv,v⟩·exp(1/2|J^1/2v|^2 + i ⟨ J^1/2v,J^-1/2λ⟩- 1/2|J^-1/2λ|^2)/√( C^1/2JC^1/2)
       = e^-I(v)·exp(i ⟨ v,λ⟩ - 1/2⟨ J^-1λ, λ⟩)/√( C^1/2JC^1/2)

This proves equation (<ref>) and it follows that the characteristic function of the measure ν̃ defined by ν̃/μ_0 = 1/Z·exp(-TΦ) fulfills

    ν̂̃̂(λ)    = ∫exp(i⟨λ, u⟩) ν̃(u)/Z = ∫exp(i⟨λ, u⟩)·exp(-TΦ(u)) μ_0(u)/∫exp(i⟨λ, u⟩) μ_0(u)
       = exp{i⟨, λ⟩ + 1/2⟨I̋()^-1(u),u⟩},

i.e. ν̃= N(, I̋()^-1) as claimed.

As in the one-dimensional setting, we conclude

    μ/μ_0 = exp(-Φ(u))/∫exp(-Φ(u))μ_0(u)

and

    ν/μ_0 = exp(-TΦ(u))/∫exp(-TΦ(u))μ_0(u).


With these expressions, we derive a bound on the Hellinger distance between  μ and ν in the next section.


§ HELLINGER DISTANCE


There are many notions of metrics and semi-metrics between measures, notably total variation, Hellinger, Wasserstein, Prokhorov and Kullback-Leibler. A survey of these and more probability metrics, including a detailed exposition of their relations is <cit.>. 

We choose the Hellinger distance, mainly because of its good analytic properties, its consistency with the total variation metric and because of the following lemma which allows us to bound the difference between expectations under the different measures in question:


Let μ, μ' be two probability measures which are absolutely continuous w.r.t. another measure ν on a Banach space (X,·_X). Assume that f: X→ E, where (E,·) has second moments with respect to both μ and μ'. Then

    ^μ f - ^μ' f≤ 2 √(^μf^2 + ^μ'f^2)· d_H(μ,μ').



Recall that


    d_H(μ,ν)^2 = 1 - ∫√(μ/μ_0(u))√(ν/μ_0(u))μ_0(u).

Thus, with the results of the preceding sections,

    d_H(μ,ν)^2    = 1 - ∫exp(-Φ(u) + TΦ(u)/2)μ_0(u)/√(∫exp(-Φ(u))μ_0(u))√(∫exp(-TΦ(u))μ_0(u))
       = 1 - ⟨exp(-1/2Φ), exp(-1/2TΦ)⟩_L^2(X,μ_0)/exp(-1/2Φ)_L^2(X,μ_0)·exp(-1/2TΦ)_L^2(X,μ_0)



From positivity of the exponential and the Cauchy-Schwarz-Bunyakowsky inequality it can easily be seen that d_H(μ,ν)∈ [0,1].

We would like to bound the Hellinger distance, i.e. optimally we would like to prove something like d_H(μ,ν) ≤ϵ for some ϵ > 0. This amounts to a reverse Cauchy-Schwarz inequality, or a statement of the kind

    ⟨ e^-Φ/2, e^-(TΦ)/2⟩_L^2(X,μ_0) > (1-^2) ·e^-Φ/2_L^2(X,μ_0)·e^-(TΦ)/2_L^2(X,μ_0).

In the next section, we present a few elementary results about this kind of inequality.


§ A REVERSE CAUCHY-SCHWARZ INEQUALITY

In this section, H will always be a Hilbert space with inner product ⟨·⟩ and norm ·. The study of of reverse Cauchy-Schwarz-Bunyakowsky inequalities is an active field of research by itself. We refer to <cit.> for further reading, although we will need only a very basic form of a reverse CSB inequality.




















Let f,g∈ H and D>0 with f-g^2 ≤ D· (f^2 + g^2). Then

    ⟨ f, g ⟩≥1-D/2· (f^2+g^2)≥ (1-D)·f·g.
 



    ⟨ f, g ⟩   - 1-D/2·f^2 - 1-D/2·g^2 
       = 1/2f^2 + 1/2g^2 - 1/2f-g^2- 1-D/2·f^2 - 1-D/2·g^2
       = 1/2[D· (f^2+g^2)-f-g^2] ≥ 0.

The last inequality in (<ref>) is just Young's inequality in .

 Lemma <ref> can be thought of as a reverse Young's inequality (which gives for all f,g∈ H that ⟨ f, g⟩≤1/2f^2 + 1/2f^2 and thus if the conditions on f,g as in lemma <ref> are fulfilled, we can bound

    1-D/2· (f^2+g^2) ≤⟨ f, g⟩≤1/2· (f^2+g^2).
 


Let f,g∈ H and K ∈ (0,1) with f-g≤ K·f. Then

    ⟨ f,g⟩≥[1-K/1 + (1-K)^2]·(f^2+g^2) ≥ 2 ·[1-K/1 + (1-K)^2] ·f·g.




    f-g^2/f^2+g^2   ≤f-g^2/f^2 + (f-f-g)^2= f-g^2/2f^2 - 2ff-g + f-g^2
       ≤f-g^2/2f^2 - 2Kf^2 + f-g^2= 1 - (2-2K)f^2/(2-2K)f^2 + f-g^2
       ≤ 1 - 2(1-K)/1 + (1-K)^2

And using lemma <ref> (with D = 1 - 2(1-K)/1 + (1-K)^2), we obtain the result.


This is another form of a reverse Young's inequality, but with a different prerequisite: If f-g≤ K·f, we have

    1/2·(1 - K^2/1 + (1-K)^2) · (f^2 + g^2) ≤⟨ f,g⟩≤1/2· (f^2+g^2)



















§ CONDITIONS FOR GOOD APPROXIMATION

From lemma <ref> and the expression for _H(μ, ν) in (<ref>) (the Hellinger distance between the posterior measure μ and its Laplace approximation ν) we immediately obtain the following:
























































With the assumptions of proposition <ref>, if for some K∈ (0,1)

    exp(-Φ(u)/2) - exp(-TΦ(u)/2)_L^2(X, μ_0)   ≤ K·exp(-I())/√((C_0^1/2·I̋()· C_0^1/2))
       = K·exp(-I())/√((Id + C_0^1/2· HΦ ()· C_0^1/2))

then 

    d_H(μ,ν) ≤K/√(1 + (1-K)^2)



We use lemma <ref>, where H = L^2(X,μ_0), f=exp(-TΦ/2) and g = exp(-Φ/2).  Note that f_H^2 = exp(-I())/√((C_0^1/2·I̋() · C_0^1/2)) due to (<ref>). Then we can set 1-^2 = 2  [1-K/1 + (1-K)^2], hence = K/√(1 + (1-K)^2) in equation (<ref>).

The equality in condition  (<ref>) follows from I(u) = Φ(u) + 1/2u_C_0.






Note that 

    C_0^1/2HΦ(u)C_0^1/2[h_1,h_2]    = ⟨ DG(u)(C_0^1/2h_1), Γ^-1DG(u)(C_0^1/2h_2)⟩
       - ⟨ HG(u)[C_0^1/2h_1, C_0^1/2h_2], Γ^-1(y-G(u))⟩

where DG is the Frechet derivative of the forward operator of G and HG is its second Frechet derivative (the Hessian).


The following corollary uses assumptions which can be checked more easily but is less strict.

Define 

    K^2:= √((C_0^1/2I̋()C_0^1/2))/exp(-I())·∫exp (-min{Φ(u),TΦ(u)})·min{|Φ(u)-TΦ(u)|^2/4,1}μ_0( u).

Then we have

    d_H(μ,ν) ≤K/√(1 + (1-K)^2).



This is due to the elementary inequality |e^-x-e^-y|≤ e^-(x∧ y)· (|x-y|∧ 1) from which we obtain

    ∫ |e^-Φ(u)/2-e^-TΦ(u)/2|^2μ_0    ≤∫ e^-(Φ(u)∧ TΦ(u))· (|Φ(u)-TΦ(u)|^2/4 ∧ 1)μ_0
       =K^2 ·∫ e^-TΦ(u)μ_0.




§ A CAVEAT ON APPLICATION


We check on a one-dimensional example whether the bounds on the Hellinger distance obtained are tight. We consider a one-dimensional example by setting G:→, with G(u) = exp(u). This means that Φ(u) = |y-e^u|^2/2.

We consider two cases: y=-2 and y=2. In the first case, the Laplace approximation is not very good (d_H(μ,ν) ≈ 0.3260) whereas in the second case, there is only a slight misfit between posterior measure and its Laplace approximation (d_H(μ, ν)≈ 0.0958).















  
The case y=-2.
Recall that the actual value for the Hellinger distance is d_H(μ,ν) ≈ 0.32595




  * The K from proposition <ref> is K≈ 0.46621, hence d_H(μ,ν)≤ 0.41128.

  * Explicit calculation of K^2 in corollary <ref> yields K≈ 0.55328, so d_H(μ, ν)≤ 0.50517.


















 


  
The case y=2.
Recall that the actual value for the Hellinger distance is d_H(μ,ν) ≈ 0.095810




  * The K from proposition <ref> is K≈ 0.13648, hence d_H(μ,ν)≤ 0.10330.

  * Explicit calculation of K^2 in corollary <ref> yields K≈ 0.17422, so d_H(μ, ν)≤ 0.13434.
















 
In conclusion, the bounds are not strict even in the one-dimensional setting (which is obvious, given the nature of inequalities we have used in lemmata <ref>, <ref>, but they give a good idea of the magnitude of the Hellinger distance nonetheless.




The calculation of the necessary quantities was done by explicit integration in order to check for strictness of the bound. This is unfeasible in higher dimensions which are of much higher interest (in one dimension there is no reason to even use a Laplace approximation instead of the fully non-Gaussian posterior). So we turn our attention to ways of obtaining results by other methods. We will see that this is difficult even in one dimension: Corollary <ref> demands calculation of the integral 

    K^2 = √((C_0^1/2I̋()C_0^1/2))/exp(-I())·∫exp (-min{Φ(u),TΦ(u)})·min{|Φ(u)-TΦ(u)|^2/4,1}μ_0( u).

As a simpler example, consider the integral

    ∫_-∞^∞ x^2 exp(-|x|) N(0,1)( x).

An intuitive approach is to split the real line in subsets A = [-,], B = [-log(1/), -] ∪ [ ,log(1/)], C = {|u|> log(1/)}. Then the quadratic term is small on A and the exponential term is small on C. The set B makes a lot more trouble as neither the quadratic nor the exponential part is consistently small and taking the maximum x^2 exp(-|x|)|_x=2 as a bound for the function on B leads to a suboptimal result. It is the same with the formula for K: For u near , the error between Φ and TΦ is small, while for large u, both the misfit Φ and the quadratic TΦ become large, so the exponential term gets small very quickly. The interface between those two domains though is hard to handle and taking an analytical bound ony any term in this interface yields a highly suboptimal bound. This makes an explicit calculation of the integral almost impossible (save taking a large number of subsets tesselating the domain, but this amounts to explicit numerical integration in the limit). This is even more true in higher (and infinite dimensions), so the practical use of this formula is very limited.

This problem makes it necessary to apply MCMC sampling to obtain the value K, but in this case we could have used to definition of d_H in the first place to get the exact Hellinger distance. 



§ ACKNOWLEDGMENTS
 The author thanks Serge Kräutle for help with reverse CSB inequalities as well as Carlo Beenakker for a helpful discussion on MathOverflow.

siam


