# Laplace's method

In mathematics, Laplace's method, named after Pierre-Simon Laplace, is a technique used to approximate integrals of the form

$\int _{a}^{b}\!e^{Mf(x)}\,dx$ where ƒ(x) is some twice-differentiable function, M is a large number, and the integral endpoints a and b could possibly be infinite. This technique was originally presented in Laplace (1774, pp. 366–367).

## The idea of Laplace's method The function e(x), in blue, is shown on top for M = 0.5, and at the bottom for M = 3. Here, ƒ(x) = sin x/x, with a global maximum at x0 = 0. It is seen that as M grows larger, the approximation of this function by a Gaussian function (shown in red) is getting better. This observation underlies Laplace's method.

Assume that the function ƒ(x) has a unique global maximum at x0. Then, the value ƒ(x0) will be larger than other values ƒ(x). If we multiply this function by a large number M, the ratio between (x0) and (x) will stay the same (since (x0)/(x) = ƒ(x0)/ƒ(x)), but it will grow exponentially in the function (see figure)

$e^{Mf(x)}.\,$ Thus, significant contributions to the integral of this function will come only from points x in a neighborhood of x0, which can then be estimated.

## General theory of Laplace's method

To state and motivate the method, we need several assumptions. We will assume that x0 is not an endpoint of the interval of integration, that the values ƒ(x) cannot be very close to ƒ(x0) unless x is close to x0, and that the second derivative $f''(x_{0})<0$ .

We can expand ƒ(x) around x0 by Taylor's theorem,

$f(x)=f(x_{0})+f'(x_{0})(x-x_{0})+{\frac {1}{2}}f''(x_{0})(x-x_{0})^{2}+R$ where $R=O\left((x-x_{0})^{3}\right).$ Since ƒ has a global maximum at x0, and since x0 is not an endpoint, it is a stationary point, so the derivative of ƒ vanishes at x0. Therefore, the function ƒ(x) may be approximated to quadratic order

$f(x)\approx f(x_{0})-{\frac {1}{2}}|f''(x_{0})|(x-x_{0})^{2}$ for x close to x0 (recall that the second derivative is negative at the global maximum ƒ(x0)). The assumptions made ensure the accuracy of the approximation

$\int _{a}^{b}\!e^{Mf(x)}\,dx\approx e^{Mf(x_{0})}\int _{a}^{b}e^{-M|f''(x_{0})|(x-x_{0})^{2}/2}\,dx$ (see the picture on the right). This latter integral is a Gaussian integral if the limits of integration go from −∞ to +∞ (which can be assumed because the exponential decays very fast away from x0), and thus it can be calculated. We find

$\int _{a}^{b}\!e^{Mf(x)}\,dx\approx {\sqrt {\frac {2\pi }{M|f''(x_{0})|}}}e^{Mf(x_{0})}{\text{ as }}M\to \infty .\,$ A generalization of this method and extension to arbitrary precision is provided by Fog (2008).

Formal statement and proof:

Then,

$\lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)=1$ Template:Hidden begin Lower bound:

Then we have the following lower bound:

$\int _{a}^{b}e^{nf(x)}\,dx\geq \int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx\geq e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})-\varepsilon )(x-x_{0})^{2}}\,dx=e^{nf(x_{0})}{\sqrt {\frac {1}{n(-f''(x_{0})+\varepsilon )}}}\int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy$ where the last equality was obtained by a change of variables $y={\sqrt {n(-f''(x_{0})+\varepsilon )}}(x-x_{0})$ . Remember that $f''(x_{0})<0$ so that is why we can take the square root of its negation.

If we divide both sides of the above inequality by $e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}$ and take the limit we get:

$\lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\geq \lim _{n\to +\infty }{\frac {1}{\sqrt {2\pi }}}\int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy{\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}$ since this is true for arbitrary $\varepsilon$ we get the lower bound:

$\lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\geq 1$ Upper bound:

Then we can calculate the following upper bound:

$\int _{a}^{b}e^{nf(x)}\,dx\leq \int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx$ $\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{-\infty }^{+\infty }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx$ $\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0})-\varepsilon )}}}$ If we divide both sides of the above inequality by $e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}$ and take the limit we get:

$\lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\leq \lim _{n\to +\infty }\left((b-a)e^{-\eta n}{\sqrt {\frac {n(-f''(x_{0}))}{2\pi }}}+{\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}\right)={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}$ $\lim _{n\to +\infty }\left({\frac {\int _{a}^{b}e^{nf(x)}\,dx}{\left(e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}\right)}}\right)\leq 1$ And combining this with the lower bound gives the result.

Note that the above proof obviously fails when $a=-\infty$ or $b=\infty$ (or both). To deal with these cases, we need some extra assumptions. A sufficient (not necessary) assumption is that for $n=1$ , the integral $\int _{a}^{b}e^{nf(x)}\,dx$ is finite, and that the number $\eta$ as above exists (note that this must be an assumption in the case when the interval $[a,b]$ is infinite). The proof proceeds otherwise as above, but the integrals

$\int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx$ must be approximated by

$\int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq \int _{a}^{b}e^{f(x)}e^{(n-1)(f(x_{0})-\eta )}\,dx=e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx$ ${\frac {e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}=e^{-(n-1)\eta }{\sqrt {n}}e^{-f(x_{0})}\int _{a}^{b}e^{f(x)}\,dx{\sqrt {\frac {-f''(x_{0})}{2\pi }}}$ whose limit as $n\rightarrow \infty$ is $0$ . The rest of the proof (the analysis of the interesting term) proceeds as above.

The given condition in the infinite interval case is, as said above, sufficient but not necessary. However, the condition is fulfilled in many, if not in most, applications: the condition simply says that the integral we are studying must be well-defined (not infinite) and that the maximum of the function at $x_{0}$ must be a "true" maximum (the number $\eta >0$ must exist). There is no need to demand that the integral is finite for $n=1$ but it is enough to demand that the integral is finite for some $n=N$ .

This method relies on 4 basic concepts such as Template:Hidden begin

1. Relative error

First of all, we need to have an understanding about the so-called “approximation” in this method is related to the relative error instead of the absolute error. Therefore, if we set

, this integration can be written as

{\begin{aligned}\int _{a}^{b}\!e^{Mf(x)}\,dx&=se^{Mf(x_{0})}{\frac {1}{s}}\int _{a}^{b}\!e^{M(f(x)-f(x_{0}))}\,dx\\&=se^{Mf(x_{0})}\int _{(a-x_{0})/s}^{(b-x_{0})/s}\!e^{M(f(sy+x_{0})-f(x_{0}))}\,dy\end{aligned}} , where $s$ is a small number when $M$ is a large number obviously and the relative error will be

$\left|\int _{(a-x_{0})/s}^{(b-x_{0})/s}e^{M(f(sy+x_{0})-f(x_{0}))}dy-1\right|.$ Now, let us separate this integration into two parts: $y\in [-D_{y},D_{y}]$ region and the rest part.

2. function $e^{M\left(f(sy+x_{0})-f(x_{0})\right)}$ will tend to $e^{-\pi y^{2}}$ around the stationary point when $M$ is large enough

Let’s look at the Taylor expansion of $M\left(f(x)-f(x_{0})\right)$ around x0 and translate x to y because we do the comparison in y-space, we will get

{\begin{aligned}M\left(f(x)-f(x_{0})\right)&={\frac {Mf''(x_{0})}{2}}s^{2}y^{2}+{\frac {Mf'''(x_{0})}{6}}s^{3}y^{3}+\cdots \\&=-\pi y^{2}+O\left({\frac {1}{\sqrt {M}}}\right).\end{aligned}} Note that $f'(x_{0})=0$ because $x_{0}$ is a stationary point. From this equation you will find that the terms higher than second derivative in this Taylor expansion is suppressed as the order of ${\frac {1}{\sqrt {M}}}$ so that $\exp \left(M\left(f(x)-f(x_{0})\right)\right)$ will get closer to the Gaussian function as shown in figure. Besides,

3. The bigger $M$ is, the smaller range of $x$ is related
4. If the integration used by Laplace’s method is converged, the contribution of the region which is not around the stationary point $x_{0}$ of the integration of its relative error will tend to zero when $M$ is increased.

Relying on the 3rd concept, even if we choose a very large Dy , sDy will finally be a very small number when $M$ is increased to a huge number. Then, how can we guarantee the integration of the rest part will tend to 0 when $M$ is large enough?

If the interval of the integration of this method is finite, we will find that no matter $f(x)$ is continue in the rest region, it will be always smaller than $m(x)$ shown above when $M$ is large enough. By the way, it will be proved later that the integration of $e^{Mm(x)}$ will tend to zero when $M$ is large enough.

Based on these four concepts, we can derive the relative error of this Laplace's method.

## Other formulations

Laplace's approximation is sometimes written as

$\int _{a}^{b}\!h(x)e^{Mg(x)}\,dx\approx {\sqrt {\frac {2\pi }{M|g''(x_{0})|}}}h(x_{0})e^{Mg(x_{0})}{\text{ as }}M\to \infty \,$ Importantly, the accuracy of the approximation depends on the variable of integration, that is, on what stays in $g(x)$ and what goes into $h(x)$ .

Template:Hidden begin First of all, let me set the global maximum is located at $x_{0}=0$ which can simplify the derivation and does not lost any important information; therefore, all the derivation inside this sub-section is under this assumption. Besides, what we want is the relative error $\left|R\right|$ as shown below

$\int _{a}^{b}\!h(x)e^{Mg(x)}\,dx=h(0)e^{Mg(0)}s\underbrace {\int _{a/s}^{b/s}{\frac {h(x)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}dy} _{1+R},$ $\left|R\right|=\left|\int _{a/s}^{b/s}A\,dy-\int _{-\infty }^{\infty }A_{0}\,dy\right|$ Owing to $\left|A+B\right|\leq |A|+|B|$ , we can separate this integration into 5 parts with 3 different types (a), (b) and (c), respectively. Therefore,

$|R|<\underbrace {\left|\int _{D_{y}}^{\infty }A_{0}dy\right|} _{(a_{1})}+\underbrace {\left|\int _{D_{y}}^{b/s}Ady\right|} _{(b_{1})}+\underbrace {\left|\int _{-D_{y}}^{D_{y}}\left(A-A_{0}\right)dy\right|} _{(c)}+\underbrace {\left|\int _{a/s}^{-D_{y}}Ady\right|} _{(b_{2})}+\underbrace {\left|\int _{-\infty }^{-D_{y}}A_{0}dy\right|} _{(a_{2})}$ $(a_{1})=\left|{\frac {1}{2{\sqrt {\pi }}}}\int _{\pi D_{y}^{2}}^{\infty }e^{-z}z^{-1/2}dz\right|<{\frac {e^{-\pi D_{y}^{2}}}{2\pi D_{y}}}.$ This means that as long as $D_{y}$ is large enough, it will tend to zero.

$(b_{1})\leq \left|\int _{D_{y}}^{b/s}\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{Mm(sy)}dy\right|$ where

$m(x)\geq {\frac {1}{M}}\ln {\frac {h(x)}{h(0)}}+g(x)-g(0)\,\,{\text{as}}\,\,x\in [sD_{y},b]$ Based on Taylor expansion, we can get

{\begin{aligned}M\left[g(sD_{y})-g(0)\right]&=M\left[{\frac {g''(0)}{2}}s^{2}D_{y}^{2}+{\frac {g'''(\xi )}{6}}s^{3}D_{y}^{3}\right]\,\,{\text{as}}\,\,\xi \in [0,sD_{y}]\\&=-\pi D_{y}^{2}+{\frac {(2\pi )^{3/2}g'''(\xi )D_{y}^{3}}{6{\sqrt {M}}|g''(0)|^{3/2}}},\end{aligned}} and

{\begin{aligned}Msg'(sD_{y})&=Ms\left(g''(0)sD_{y}+{\frac {g'''(\zeta )}{2}}s^{2}D_{y}^{2}\right),\,\,{\text{as}}\,\,\zeta \in [0,sD_{y}]\\&=-2\pi D_{y}+{\sqrt {\frac {2}{M}}}\left({\frac {\pi }{|g''(0)|}}\right)^{3/2}g'''(\zeta )D_{y}^{2},\end{aligned}} and then substitute them back into the calculation of $(b_{1})$ ; however, you can find that the remainders of these two expansions are both inversely proportional to the square root of $M$ , let me drop them out to beautify the calculation. Keeping them is better, but it will make the formula uglier.

{\begin{aligned}(b_{1})&\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}\int _{0}^{b/s-D_{y}}e^{-2\pi D_{y}y}dy\right|\\&\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}{\frac {1}{2\pi D_{y}}}\right|.\end{aligned}} Therefore, it will tend to zero when $D_{y}$ gets larger, but don't forget that the upper bound of $D_{y}$ should be considered during this calculation.

{\begin{aligned}(c)&\leq \int _{-D_{y}}^{D_{y}}e^{-\pi y^{2}}\left|{\frac {sh'(\xi )}{h(0)}}y\right|\,dy\\&<{\sqrt {\frac {2}{\pi M|g''(0)|}}}\left|{\frac {h'(\xi )}{h(0)}}y\right|_{\text{max}}\left(1-e^{-\pi D_{y}^{2}}\right)\end{aligned}} $\int e^{Mf(\mathbf {x} )}\,d\mathbf {x} \approx \left({\frac {2\pi }{M}}\right)^{d/2}|-H(f)(\mathbf {x} _{0})|^{-1/2}e^{Mf(\mathbf {x} _{0})}{\text{ as }}M\to \infty \,$ ## Laplace's method extension: Steepest descent

{{#invoke:main|main}}

In extensions of Laplace's method, complex analysis, and in particular Cauchy's integral formula, is used to find a contour of steepest descent for an (asymptotically with large M) equivalent integral, expressed as a line integral. In particular, if no point x0 where the derivative of ƒ vanishes exists on the real line, it may be necessary to deform the integration contour to an optimal one, where the above analysis will be possible. Again the main idea is to reduce, at least asymptotically, the calculation of the given integral to that of a simpler integral that can be explicitly evaluated. See the book of Erdelyi (1956) for a simple discussion (where the method is termed steepest descents).

The appropriate formulation for the complex z-plane is

$\int _{a}^{b}\!e^{Mf(z)}\,dz\approx {\sqrt {\frac {2\pi }{-Mf''(z_{0})}}}e^{Mf(z_{0})}{\text{ as }}M\to \infty .\,$ for a path passing through the saddle point at z0. Note the explicit appearance of a minus sign to indicate the direction of the second derivative: one must not take the modulus. Also note that if the integrand is meromorphic, one may have to add residues corresponding to poles traversed while deforming the contour (see for example section 3 of Okounkov's paper Symmetric functions and random partitions).

## Further generalizations

An extension of the steepest descent method is the so-called nonlinear stationary phase/steepest descent method. Here, instead of integrals, one needs to evaluate asymptotically solutions of Riemann–Hilbert factorization problems.

Given a contour C in the complex sphere, a function ƒ defined on that contour and a special point, say infinity, one seeks a function M holomorphic away from the contour C, with prescribed jump across C, and with a given normalization at infinity. If ƒ and hence M are matrices rather than scalars this is a problem that in general does not admit an explicit solution.

An asymptotic evaluation is then possible along the lines of the linear stationary phase/steepest descent method. The idea is to reduce asymptotically the solution of the given Riemann–Hilbert problem to that of a simpler, explicitly solvable, Riemann–Hilbert problem. Cauchy's theorem is used to justify deformations of the jump contour.

The nonlinear stationary phase was introduced by Deift and Zhou in 1993, based on earlier work of Its. A (properly speaking) nonlinear steepest descent method was introduced by Kamvissis, K. McLaughlin and P. Miller in 2003, based on previous work of Lax, Levermore, Deift, Venakides and Zhou.

The nonlinear stationary phase/steepest descent method has applications to the theory of soliton equations and integrable models, random matrices and combinatorics.

## Complex integrals

For complex integrals in the form:

${\frac {1}{2\pi i}}\int _{c-i\infty }^{c+i\infty }g(s)e^{st}\,ds$ with t >> 1, we make the substitution t = iu and the change of variable s = c + ix to get the Laplace bilateral transform:

${\frac {1}{2\pi }}\int _{-\infty }^{\infty }g(c+ix)e^{-ux}e^{icu}\,dx.$ We then split g(c+ix) in its real and complex part, after which we recover u = t / i. This is useful for inverse Laplace transforms, the Perron formula and complex integration.

## Example 1: Stirling's approximation

Laplace's method can be used to derive Stirling's approximation

$N!\approx {\sqrt {2\pi N}}N^{N}e^{-N}\,$ for a large integer N.

From the definition of the Gamma function, we have

$N!=\Gamma (N+1)=\int _{0}^{\infty }e^{-x}x^{N}\,dx.$ Now we change variables, letting

$x=Nz\,$ so that

$dx=N\,dz.$ Plug these values back in to obtain

{\begin{aligned}N!&=\int _{0}^{\infty }e^{-Nz}\left(Nz\right)^{N}N\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{-Nz}z^{N}\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{-Nz}e^{N\ln z}\,dz\\&=N^{N+1}\int _{0}^{\infty }e^{N(\ln z-z)}\,dz.\end{aligned}} This integral has the form necessary for Laplace's method with

$f\left(z\right)=\ln {z}-z$ which is twice-differentiable:

$f'(z)={\frac {1}{z}}-1,\,$ $f''(z)=-{\frac {1}{z^{2}}}.\,$ The maximum of ƒ(z) lies at z0 = 1, and the second derivative of ƒ(z) has the value −1 at this point. Therefore, we obtain

$N!\approx N^{N+1}{\sqrt {\frac {2\pi }{N}}}e^{-N}={\sqrt {2\pi N}}N^{N}e^{-N}.\,$ ## Example 2: parameter estimation and probabilistic inference

Template:Harvnb reviews Laplace's method results (univariate and multivariate) and presents a detailed example showing the method used in parameter estimation and probabilistic inference under a Bayesian perspective. Laplace's method is applied to a meta-analysis problem from the medical domain, involving experimental data, and compared to other techniques.