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

12-28

最近在自學關於數據結構的東西。突然想到如果要大矩陣就必須有特殊的演算法。
於是百度了一下，找到了這種方法，但是不知道怎麼用C語言實現。
問一下有什麼具體的說明或者代碼解釋么？

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

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

假設你有一個線性方程組：
$Ax=b$
其中A是已知矩陣，b是已知向量，x是需要求解的未知向量。

當你有這麼個問題需要解決時，一般的思路是直接求A的逆矩陣，然後x就出來了：
$x=A^{-1}b$

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

$A^{-1}bapprox sum_{i=0}^{m-1}{eta_{i}A^{i}b}=eta_{0}b+ eta_{1}Ab+eta_{2}A^{2}b...+eta_{m-1}A^{m-1}b$
其中 $eta$ 都是未知標量，m是你來假設的一個值，最大不能超過矩陣的維度，比如這裡例子里是1000.
瞧，這麼一處理，我們就不用算 $A^{-1}$ 了。我們只要求出方程里那些 $eta$ 的值，就齊活兒了。

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

(這裡省略了幾步，還要用Arnoldi方法做個循環，先留個空，有同學需要我再補上)

解 $eta$ 值要帶回第一個公式，得到以下方程：
$0=b-A x^{(m)}=b-Asum_{i=0}^{m-1}{eta_{i}A^{i}b}$

有細心的同學一看，說不對勁啊。b的維度是1000，那就是有1000個方程， $eta$ 的數量小於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 eta_{0}} =0, frac{partial r}{partial eta_{1}} =0, frac{partial r}{partial eta_{2}} =0,... frac{partial r}{partial eta_{m-1}} =0$
（註：謝@渭水泱泱提醒，這裡的r指的是 $r^{(m)}$ 的平方和）

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

於是問題轉化為了一個求m個方程m個未知數的方程組的問題，而且m通常不大（當然，m是你自己設定的，設那麼大不是自找麻煩么）

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

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

回顧一下，大概是這樣一個流程：
大型稀疏矩陣求逆--&>Krylov方法--&>線性方程最小二乘問題--&>小矩陣求逆

有意思的問題。
A*Inv(A) = I，假設Inv(A) = {V1, V2, V3, ...}， I = {I1, I2, I3, ...}
然後 A*V1 = I1, ...
在CG/BiCGStab中利用A*{V1, V2, V3, ...} = {I1, I2, I3...}。
這種方法效率不是很好，畢竟A*V1 = I1, A*V2 = I2, ...的收效速度不一樣。
如果A不太大，利用LU/QR分解，倒是個好辦法。