拉普拉斯方法

在數學上，以皮埃爾 - 西蒙・拉普拉斯命名的**拉普拉斯方法**是用於得出下列積分形式的近似解的方法：

　
: $\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_0 為極大值所在處。這方法最早是拉普拉斯在 (1774, pp. 366–367) 所發表。(待考查)

## 拉普拉斯方法的想法概論

基於上述四點，就有辦法證明拉普拉斯方法的可靠性。而 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$

## 其他形式

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

　
: $\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)$ 有關。 \*\*\* 待考查 \*\*\*

至於多維的情形，讓我們令 $\mathbf {x} $ 是一個 $n$ 維向量，而 $f(\mathbf {x} )$ 則是一個純量函數，則此拉普拉斯方法可以寫成

　
: $\int e^{Mf(\mathbf {x} )}\,d^{n}x\approx \left({\frac {2\pi }{M}}\right)^{n/2}|-H(f)(\mathbf {x} _{0})|^{-1/2}e^{Mf(\mathbf {x} _{0})}{\text{ as }}M\to \infty \,$

其中的 $H(f)(\mathbf {x} _{0})$ 是一個在 $\mathbf {x} _{0}$ 取值的 海森矩陣，而 $|\cdot |$ 則是指矩陣的 行列式 ；此外，與單變數的拉普拉斯方法類似，這裡的 海森矩陣 必須為 負定矩陣，即該矩陣的所有本徵值皆小於 0，這樣才會是極大值所在。.

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

此拉普拉斯方法可以被推廣到 複分析 裡頭使用，搭配 留數定理 ，以找出一個過最速下降點的 contour (翻路徑的話，會與 path integral 相衝，所以，這裡還是以英文原字稱呼) 的 曲線積分 ，用來取代原有的複數空間的 contour 積分。因為有時 $f(z)$ 的一階微分為 0 的點未必就在實數軸上，而就算真在實數軸上，也未必二階微分在 $x$ 方向上為小於 0 的實數；此種情況下，就得使用最速下降法了。由於最速下降法中，已經利用另一條通過最速下降的鞍點來取代原有的 contour 積分，經過變數變換後就會變得有如拉普拉斯方法，因此，我們可以透過這新的 contour ，找到原本的積分的漸進近似解，而這將大大的簡化整個計算。就好像原本的路徑像是在蜿蜒的山路開車，而新的路徑就像乾脆繞過這座山開，反正目的只是到達目的地而已，留數定理已經幫我們把中間的差都算好了。請讀 Erdelyi (2012) 與 Arfken & Weber (2005) 的書裡有關 _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_0 就是新的路徑通過的鞍點。
注意，開根號裡的負號是用來指定最速下降的方向，千萬別認為取 $f''(z_{0})$ 的絕對值來取代這個負號，若然，那就大錯特錯了。 
另外要注意的是，如果該被積函數是 亞純函數 ，就有必要將被新舊 contour 包到的極點所貢獻的留數給加入，範例請參考 Okounkov 的文章 arXiv:math/0309074 的第三章。

## 更進一步一般化

最速下降法還可以更進一步的推廣到所謂的 _非線性穩定相位 / 最速下降法_ (nonlinear stationary phase/steepest descent method)。
這方法主要用在解非線性偏微分方程，透過將非線性偏微分方程轉換為求解柯西變換 (Cauchy transform) 的積分形式，就可以藉由最速下降法的想法來得到非線性解的漸進解。

以 艾里方程 (線性) $y''=xy$ 為例，它可以寫成積分形式如下：

　
: $y_{j}(x)={\frac {1}{2\pi i}}\int _{\gamma _{j}}e^{{\frac {1}{3}}t^{3}-xt}dt$

由這條積分式，我們就可以藉由最速下降法 (若 $\gamma _{j}$ 指的是負實數軸，那麼就回到此拉普拉斯方法了) 來得到它的漸進解了。

然而，若方程式如 KdV 方程 $y_{t}+6yy_{x}+y_{xxx}=0$ 是個非線性偏微分方程，想要找到它的解相對應的一個複數 contour 積分的話，就沒那麼簡單，在非線性穩定相位 / 最速下降法裡所用到的概念主要是基於散射逆轉換 (inverse scattering transform) 的處理方式，即先將原本的非線性偏微分方程變成 Lax 對 ，其中一個像是線性的 薛丁格方程式 ，其位能障為我們要找的 $y$ ，本徵值為 $\lambda $ ，波函數為 $\phi _{\lambda }$ (不過，它並非我們所要的 $y$ )；因此可以解它的散射矩陣，若利用解析延拓將原本的波函數由實數 $\lambda $ 延拓到複數空間時，就可以得到黎曼希爾伯特問題 (RHP) 的形式。利用這個黎曼希爾伯特問題 (RHP) ，我們可以解得 $\phi $ 的柯西變換的積分形式，再利用此線性薛丁格方程的特性，就可以反推出 $y$ 的複數 contour 積分 形式了。

而 Lax 對 的另一個偏微分方程則是決定每個 $\phi _{\lambda }$ 隨時間變化的行為，由於 $y$ 在 $x\rightarrow \pm \infty $ 時被要求為 0 ，會發現整個偏微分方程會變得十分簡單，並且只決定 $c(\lambda ,t)\phi _{\lambda }$ 裡的 $c(\lambda ,t)$ 的值，不過，條件是 $\phi _{\lambda }$ 必須是指由正負無限遠入射或反射波的解。這樣，我們所得到的這個只與時間與本徵值有關的係數 $c(\lambda ,t)$ 就可以直接被應用在上述的黎曼希爾伯特問題 (RHP) 裡的躍變矩陣裡了。

接著就是非線性穩定相位 / 最速下降法所要做的工作，即找出 鞍點 來，在該點附近基於最速下降法的精神做近似。不過，這近似因著考慮到收斂性，需要將原本的 contour 變形，與將原本的黎曼希爾伯特問題 (RHP) 作轉換，所以有再比原本的最速下降法多出一些步驟來。

這整個方法最早由 Deift 與 Zhou 在 1993 基於 Its 之前的工作所提出的，後續又有許多人加以改進，主要的應用則有 孤波 理論，可積模型等。

## 複數積分

以下的積分常用在 拉普拉斯變換 #拉普拉斯逆變換 裡，

　
: ${\frac {1}{2\pi i}}\int _{c-i\infty }^{c+i\infty }g(s)e^{st}\,ds$

假定我們想要得到該積分在 $t>>1$ 時的結果 (若 $t$ 為時間，通常就是在找經過夠久時間後達穩定的結果)，我們可以透過 解析延拓 的概念，先將這時間換成虛數，如 _t_ = _iu_ 並且一併做 $s=c+ix$ 的變換，則我們可以將上式轉換為如下的 拉普拉斯變換 #雙邊拉普拉斯變換

　
: ${\frac {1}{2\pi }}\int _{-\infty }^{\infty }g(c+ix)e^{-ux}e^{icu}\,dx.$

這裡就可以套用此拉普拉斯方法求漸進解，最後，再利用 _u_ = _t_ / _i_ 把 _t_ 換回來，就可以得到該拉普拉斯逆變換的漸進解了。

## 例子 1：斯特靈公式

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

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

 **證明：** 

由 Γ 函數 的積分定義，我們可以得到

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

接著讓我們做變數變換，

　
: : $x=Nz\,$

因此

　
: : $dx=N\,dz.$

將這些代回 Γ 函數 的積分定義裡，我們可以得到

　
: ${\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_0= 1 而且在該點的二次微分為 $-1$ 。因此，我們得到

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

## 例子 2：貝氏網路 ，參數估計與概率推理

關於概率推理的簡介，請參考 http://doc.baidu.com/view/e9c1c086b9d528ea81c77989.html 。
而在 Azevedo-Filho & Shachter 1994 的文章裡，則回顧了如何應用此拉普拉斯方法 (無論是單變數或者多變數) 如何應用在 _概率推理_ 上，以加速得到系統的 後驗矩 (數學) (posterior moment) ， 貝氏參數 等，並舉醫學診斷上的應用為例。

## 相關維基百科文章

* 平穩相近似

## 參考文獻

* Deift, P.; Zhou, X., A steepest descent method for oscillatory Riemann–Hilbert problems. Asymptotics for the MKdV equation, Ann. of Math. **137** (2), 1993, **137** (2): 295–368, doi:10.2307/2946540 .
* Azevedo-Filho, A.; Shachter, R., Laplace's Method Approximations for Probabilistic Inference in Belief Networks with Continuous Variables, Mantaras, R.; Poole, D. (編), Uncertainty in Artificial Intelligence, San Francisco, CA: Morgan Kauffman, 1994, CiteSeerX: 10.1.1.91.2064 .
* Erdelyi, A., Asymptotic Expansions, Dover Publications, 2012, ISBN 9780486155050 .
* Arfken, G.B.; Weber, H.J., Mathematical Methods For Physicists International Student Edition, Elsevier Science, 2005, ISBN 9780080470696 .