拉普拉斯方法

在数学上，以皮埃尔-西蒙·拉普拉斯命名的拉普拉斯方法是用于得出下列积分形式的近似解的方法：

$\int _{a}^{b}\!e^{Mf(x)}\,dx$

其中的 ƒ(x) 是一個二次可微函数， M 是一個很大的數，而積分邊界點 a 與 b 則允許為無限大。此外，函數 ƒ(x) 在此積分範圍內的 全域極大值 所在處必須是唯一的並且不在邊界點上。則它的近似解可以寫為

$\int _{a}^{b}\!e^{Mf(x)}\,dx\approx {\sqrt {\frac {2\pi }{M\left|f''(x_{0})\right|}}}e^{Mf(x_{0})}{\text{ as }}M\to \infty .\,$

其中的 x₀ 為極大值所在處。這方法最早是拉普拉斯在 (1774, pp. 366–367) 所發表。(待考查)

拉普拉斯方法的想法概論

1. 相對誤差

首先，我們得先知道的是，這裡所謂的近似解是以 相對誤差 在講，而不是指 絕對誤差 。因此，令 $s={\sqrt {\frac {2\pi }{M\left|f''(x_{0})\right|}}}$ ， 此 s 很顯然的當 M 很大時為一個微小的數，則上述的積分可以寫為 ${\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}}$ 因此，我們的相對誤差很顯然的為 $\left|\int _{(a-x_{0})/s}^{(b-x_{0})/s}e^{M(f(sy+x_{0})-f(x_{0}))}dy-1\right|$ 現在，讓我們把積分分為 $y\in [-D_{y},D_{y}]$ 的與餘下的部分。

2. 函數 $e^{M\left(f(sy+x_{0})-f(x_{0})\right)}$ 在極大值附近的行為，當 $M$ 很大時會趨近於 $e^{-\pi y^{2}}$

讓我們看一下 $M\left(f(x)-f(x_{0})\right)$ 的 泰勒展開 吧。這裡展開的原點取 x0 ，並將變數由 x 換為 y 因為我們是在對 y 做積分，再加上 x0 為函數極大值所在處，因此， ${\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}}$ 會發現高於二次微分的部分當 M 增大時，會被壓抑掉，而此 指數函數 則會越來越接近 高斯函數 $e^{-\pi y^{2}}$ 且 $\int _{-\infty }^{\infty }e^{-\pi y^{2}}dy=1$ 。如下圖所示： 圖為 $e^{M[f(sy+x_{0})-f(x_{0})]}$ ，其中M為 1, 2 與 3。紅色線為 $e^{-\pi y^{2}}$

3. 當 $M$ 增加，則相對應的 $x$ 範圍變小
$y\in [-D_{y},D_{y}]$ 的部分由上面論述可知當 M 越大則越接近高斯函數；並且有趣的是，相對應的 x 區間則越來越小，因為 $x\in [-sD_{y},sD_{y}]$ 與 ${\sqrt {M}}$ 成反比。

4.若拉普拉斯方法所用的積分會收斂，它的百分誤差相對應的非極大值附近的其餘部分的積分就會隨 $M$ 增大而趨近於0

所以，即便 Dy 取非常大， sDy 最後還是會變得非常小，那麼，我們怎麼保證餘下的部分的積分也會趨近到0？ 簡單的想法是，設法找一個函數 $m(x)$ ，使它在餘下的區域內 $m(x)\geq f(x)$ ，並且隨著 $M$ 的增大，$e^{Mm(x)}$ 的積分會趨近於0 即可。而指數函數只要指數部分為實數，它就一定大於0，並且當指數變大時，指數函數也會跟著變大。基於這兩個理由， $e^{Mf(x)}$ 的積分也會趨近於0。 在此，我們可以選擇此 $m(x)$ 為過 $x=sD_{y}$ 的切線函數，如下圖所示： $m(x)$ 以斜直線代表，為通過位於 $x=sD_{y}+x_{0}$ 的切線，當 $sD_{y}$ 越小，涵蓋的範圍越大 若為有限積分區間，我們可以看出，只要 $M$ 夠大，這切線的斜率就會趨近於0 ，因此，不管 $f(x)$ 在其餘部分是否連續，它們總會在這切線以下，並且很容易證明(後面會證) $e^{Mm(x)}$ 的積分當 $M$ 很大時，就會趨近於0。 然而，若積分範圍是無限大的話，那麼，$m(x)$ 就有可能總是與 $f(x)$ 有相交，若有相交，則表示無法保證 $e^{Mf(x)}$ 的積分最後會趨近於0。比方 $f(x)={\frac {sin(x)}{x}}$ 的情形下，不管 $M$ 為多少，它的 $\int _{0}^{\infty }e^{Mf(x)}dx$ 總是會發散。所以，對於積分範圍為無限大的情形，還得要求 $\int _{d}^{\infty }e^{Mf(x)}dx$ 會收斂才行，那麼，只要 $d$ 越來越大，該積分就會趨近於0。而這 $d$ 可選為$m(x)$ 與 $f(x)$ 相交之處即可，即保證所有餘下的部分的積分會趨近於0。 也許您會問說，為何不是要求 $\int _{d}^{\infty }e^{f(x)}dx$ 會收斂？舉個例子，若 $f(x)$ 在餘下區域內為 $-\ln x$ ，則 $e^{f(x)}={\frac {1}{x}}$ ，積分會發散，然而，當 $M=2$ ， $e^{Mf(x)}={\frac {1}{x^{2}}}$ 的積分卻是收斂的。所以，有的函數在 $M$ 小的時候會發散，然而，當 $M$ 夠大時卻會是收斂的。

基於上述四點，就有辦法證明拉普拉斯方法的可靠性。而 Fog(2008) 又將此方法推廣到任意精確。 ***待考查***

此方法的正式表述與證明：

假設 $f(x)$ 是一個在 $x_{0}\in (a,b)$ 這點滿足 (1) $f(x_{0})=\max _{[a,b]}f(x)$ ，(2)唯一全域最大，(3) $x_{0}$ 附近為二階可微且 $f''(x_{0})<0$ (4) 當拉普拉斯方法的積分範圍為無限大時，此積分會收斂，

則，

$\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$

證明部分請參閱 Laplace's method 。 ***待翻譯***

其他形式

有時拉普拉斯方法也會被寫成其他形式，如：

$\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 \,$

其中 $h$為正 (好像不必要)。

重要的是，這方法精確度與函數 $g(x)$ 和 $h(x)$ 有關。 ***待考查***

百分誤差的推導

首先，由於極大值所在設為 $x_{0}=0$ 並不影響計算結果，所以以下的推導皆如此假設。此外，我們想要的是百分誤差 $\left|R\right|$，所以 $\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}$ 其中 $s\equiv {\sqrt {\frac {2\pi }{M\left|g''(0)\right|}}}$ ，所以，若我們令 $A\equiv {\frac {h(sy)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}$ 且 $A_{0}\equiv e^{-\pi y^{2}}$ ，則 $\left|R\right|=\left|\int _{a/s}^{b/s}A\,dy-\int _{-\infty }^{\infty }A_{0}\,dy\right|$ 所以，只要我們求得上式的上限為何，即可得百分誤差。注意，這裡的推導還可以再優化，不過，由於此處我只想顯示它會收斂到0，因此，不再細推。 由於 $\left|A+B\right|\leq |A|+|B|$ ，因此 $|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})$ 與 $(a_{2})$ 雷同，所以只算 $(a_{1})$ ，經過 $z\equiv \pi y^{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}}}$ 也就是說，只要 $D_{y}$ 取得夠大，它就會趨近於0。 而 $(b_{1})$ 與 $(b_{2})$ 也雷同，所以也只算一個： $(b_{1})\leq \left|\int _{D_{y}}^{b/s}\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{Mm(sy)}dy\right|$ 其中 $m(x)\geq {\frac {1}{M}}\ln {\frac {h(x)}{h(0)}}+g(x)-g(0)\,\,{\text{as}}\,\,x\in [sD_{y},b]$ 並且 $h(x)$ 在此範圍內要與 $h(0)$ 同號。 這裡讓我們選過 $x=sD_{y}$ 的切線為 $m(x)$ 吧！即 $m(sy)=g(sD_{y})-g(0)+g'(sD_{y})\left(sy-sD_{y}\right)$ ，如圖所示： $m(x)$ 以斜直線代表，為通過位於 $x=sD_{y}$ 的切線 由圖可以看出，當 $s$ 或者 $D_{y}$ 變小時，滿足上列不等式的區間會變大，因此，為了涵蓋整個區間， $D_{y}$ 有了上限。此外，$e^{-\alpha x}$ 的積分也相對容易，因此，很適合用來預估誤差。 由泰勒展開我們得知， ${\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}}$ 且 ${\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}}$ 然後將它們代回 $(b_{1})$ 的計算裡，不過，您可以看到，上述的兩個餘項皆與 $M$ 的開根號成反比，為了式子的漂亮，讓我省略它們吧！不省略只是看起來較醜陋而已，不過，那樣子較嚴謹便是。 ${\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}}$ 所以，原則上也是 $D_{y}$ 越大則越趨近於0。只不過，在決定 $m(x)$ 的過程，設下了 $D_{y}$ 的上限。 至於靠近極大值的點的積分，我們一樣可以利用泰勒展開來求，當 $h(x)$ 在此點的一階微分不為0時， ${\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}}$ 會看到它原則上與 $M$ 的開根號成反比，其實，當 $h(x)$ 為常數時，積分出來的也會有如此的行為。 所以，概括的講，靠近極大點的附近的積分原則上會隨著 ${\sqrt {M}}$ 的增大而變小，而其餘的部分，只要 $D_{y}$ 夠大，也會趨近於0，只是這$D_{y}$ 是有上限的，取決於我們所找的上限函數 $m(x)$ 是否在這區間內總是大於 $g(x)-g(0)$ ；不過，一旦有一個滿足條件的 $m(x)$ 被找到，由於切線的起點為 $x=sD_{y}$ ，因此，$D_{y}$ 可以取正比於 ${\sqrt {M}}$ 即可，所以， $M$ 越大，$D_{y}$ 也可以設越大。

至於多維的情形，由於我對於 Laplace's method 裡所用的 $d\mathbf {x}$ 的定義有疑慮，暫不翻譯。

拉普拉斯方法的推廣：最速下降法

主条目：method of steepest descent

此拉普拉斯方法可以被推廣到 複分析 裡頭使用，搭配 留數定理 ，以找出一個過最速下降點的 contour (翻路徑的話，會與path integral 相衝，所以，這裡還是以英文原字稱呼) 的 曲線積分 ，用來取代原有的複數空間的 contour積分。因為有時 $f(z)$ 的一階微分為0的點未必就在實數軸上，而就算真在實數軸上，也未必二階微分在 $x$ 方向上為小於0 的實數；此種情況下，就得使用最速下降法了。由於最速下降法中，已經利用另一條通過最速下降的鞍點來取代原有的 contour 積分，經過變數變換後就會變得有如拉普拉斯方法，因此，我們可以透過這新的 contour ，找到原本的積分的漸進近似解，而這將大大的簡化整個計算。就好像原本的路徑像是在蜿蜒的山路開車，而新的路徑就像乾脆繞過這座山開，反正目的只是到達目的地而已，留數定理已經幫我們把中間的差都算好了。請讀 Erdelyi (1956) 的書裡有關 steepest descents 的章節 ***未考證***。

以下就是該方法在z 平面下的形式：

$\int _{a}^{b}\!e^{Mf(z)}\,dz\approx {\sqrt {\frac {2\pi }{-Mf''(z_{0})}}}e^{Mf(z_{0})}{\text{ as }}M\to \infty .\,$

其中 z₀ 就是新的路徑通過的鞍點。 注意，開根號裡的負號是用來指定最速下降的方向，千萬別認為取 $f''(z_{0})$ 的絕對值來取代這個負號，若然，那就大錯特錯了。 另外要注意的是，如果該被積函數是 meromorphic ，就有必要將被新舊 contour 包到的極點所貢獻的留數給加入， (請看example section 3 of Okounkov's paper Symmetric functions and random partitions). ***待考查***

更進一步一般化

***待翻譯***

複數積分

***待翻譯***

例子1：斯特靈公式

拉普拉斯方法可以用在推導 斯特靈公式 上；當 N 很大時，

$N!\approx {\sqrt {2\pi N}}N^{N}e^{-N}\,$

證明：

由 Gamma function 的積分定義，我們可以得到

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

接著讓我們做變數變換，

$x=Nz\,$

因此

$dx=N\,dz.$

將這些代回 Gamma function 的積分定義裡，我們可以得到

${\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}}$

經由此變數變換後，我們有了拉普拉斯方法所需要的$f(x)$ 為

$f\left(z\right)=\ln {z}-z$

而它乃為二次可微函數，且

$f'(z)={\frac {1}{z}}-1,\,$

$f''(z)=-{\frac {1}{z^{2}}}.\,$

因此， ƒ(z) 的極大值出現在 z₀ = 1 而且在該點的二次微分為 $-1$ 。因此，我們得到

$N!\approx N^{N+1}{\sqrt {\frac {2\pi }{N}}}e^{-N}={\sqrt {2\pi N}}N^{N}e^{-N}.\,$

例子2：參數估計與概率推論

***待翻譯***

相關維基百科文章

• Method of stationary phase

參考文獻

請參閱 Laplace's method 的 References。