當前位置: 華文問答 > 科學

如何使用Krylov方法求解矩陣的運算尤其是逆?

2014-04-06科學

不太懂數據結構上的東西,我就從數學上講講Krylov方法,希望能給你一點啟發:

Krylov方法是一種 「降維打擊」 手段,有利有弊。其特點一是犧牲了精度換取了速度,二是在沒有辦法求解大型稀疏矩陣時,他給出了一種辦法,雖然不精確。

假設你有一個線性方程式組:

Ax=b

其中A是已知矩陣,b是已知向量,x是需要求解的未知向量。

當你有這麽個問題需要解決時,一般的思路是直接求A的逆矩陣,然後x就出來了:

x=A^{-1}b

但是,如果A的維度很高,比方說1000*1000的矩陣,那麽A就是一個大型矩陣,大型矩陣是很難求逆的,如果A還是一個稀疏矩陣,那就更難求了。這時聰明的Krylov想到了一種方法來替換A^{-1}

A^{-1}b\approx \sum_{i=0}^{m-1}{\beta_{i}A^{i}b}=\beta_{0}b+ \beta_{1}Ab+\beta_{2}A^{2}b...+\beta_{m-1}A^{m-1}b

其中\beta 都是未知純量,m是你來假設的一個值,最大不能超過矩陣的維度,比如這裏例子裏是1000.

瞧,這麽一處理,我們就不用算A^{-1} 了。我們只要求出方程式裏那些\beta 的值,就齊活兒了。

(Krylov透過數學上的推導證明了,當m趨近於矩陣維度時(這裏是1000),算出來的值就是精確解了。當然很少有人會真的把m提到那個數量級來算,那樣就等於新構建了一個大型線形方程式組,計算量還是很大。不過這麽轉換一下也不是沒有好處,畢竟從稀疏矩陣變為了非稀疏矩陣,好求一點,沒準就能直接求逆了。)



(這裏省略了幾步,還要用Arnoldi方法做個迴圈,先留個空,有同學需要我再補上)


解\beta 值要帶回第一個公式,得到以下方程式:

0=b-A x^{(m)}=b-A\sum_{i=0}^{m-1}{\beta_{i}A^{i}b}

有細心的同學一看,說不對勁啊。b的維度是1000,那就是有1000個方程式,\beta 的數量小於1000. 那不是方程式數大於未知數了嗎?這種情況應該沒法兒求解啊。

對的,這種情況確實沒法兒精確求解,只能求近似解。

方程式數大於未知數時常用的方法之一是最小平方法。那麽這裏可不可以用最小平方法呢?

一般來說,最小平方法套用的最重要的條件之一,就是方程式須是線性的,最小平方法一般只用來解線性方程式,解非線性的就非常困難,需要進行一些「魔改」,比如基於最小平方法的Levenberg-Marquardt and trust-region methods,就是matlab裏的fsolve函式呼叫的演算法,這裏我就不鋪開講了,免得讀者分心。我們觀察了一下這個方程式,正好就是線性的,那麽就可以用。

(岔個話,非線性方程式組的求解一直是個「老大難」的問題,一般可用的方法只有Newton(牛頓)法,對就是三百年前英國那個牛頓,這麽些年一直沒啥進步。我們研究Krylov方法,其最重要最廣泛的套用,就是可以跟Newton法結合起來,把牛頓法裏一般需要手動求解的一個非常復雜的Jacobian矩陣給省去了。創造這一天才結合的科學家將這種耦合演算法稱作JFNK,就是Jacobian-free Newton Krylov的縮寫,意圖一目了然,從此科學家們省去了手推Jacobian矩陣的煩惱,人人用了都說好,所以學Krylov演算法的同學不順便學一學JFNK就是「入寶山而空手歸」了。)

令r^{(m)}=b-A x^{(m)}

這裏r^{(m)} 是指當m為m時的殘量,所謂殘量,就是error,就是我們不想要他存在的一個量。從上面的第一個公式就可以看出來,如果我們最終得出的x完全精確,那麽r應該等於0. 於是現在這個問題轉變為求一個含有多個自變量的運算式的最小值問題。

含有多個自變量的運算式的最小值問題,可以用最小平方法來解決。最小平方法的核心就是以下這些個公式:

\frac{\partial r}{\partial \beta_{0}} =0, \frac{\partial r}{\partial \beta_{1}} =0, \frac{\partial r}{\partial \beta_{2}} =0,... \frac{\partial r}{\partial \beta_{m-1}} =0

(註:謝@

渭水泱泱

提醒,這裏的r指的是r^{(m)} 的平方和)

意思就是在r為最小值的時候,r關於所有變量的偏導都應當為0,這是毫無疑問的。

於是問題轉化為了一個求m個方程式m個未知數的方程式組的問題,而且m通常不大(當然,m是你自己設定的,設那麽大不是自找麻煩麽)

這種問題就很好解了,一般用前面的x=A^{-1}b 方法就可以搞定了。

然後問題解決,戰鬥結束。

回顧一下,大概是這樣一個流程:

大型稀疏矩陣求逆-->Krylov方法-->線性方程式最小平方問題-->小矩陣求逆