山田的金融日記(7)-LU分解,逆矩陣

05-08

繼續C#Project 記錄

Project的目的是做各種回歸，線性回歸需要逆矩陣，逆矩陣需要LU分解，所以就從LU分解開始。

什麼是LU分解：就是將一個方塊矩陣A分解為一個下三角和上三角矩陣的乘積，自然就是L和U的由來，Lower和Upper，即A=LU;

${displaystyle {egin{bmatrix}a_{11}&a_{12}&a_{13}\a_{21}&a_{22}&a_{23}\a_{31}&a_{32}&a_{33}end{bmatrix}}={egin{bmatrix}l_{11}&0&0\l_{21}&l_{22}&0\l_{31}&l_{32}&l_{33}end{bmatrix}}{egin{bmatrix}u_{11}&u_{12}&u_{13}\0&u_{22}&u_{23}\0&0&u_{33}end{bmatrix}}}$

為什麼要用LU分解，是為了在求解Ax=b的線性方程組的時候分解為LUx=b計算。這樣先算L(Ux)=b得出Ux，然後再算x會降低演算法複雜度。

引用 @PseudoSuffix 的回答

LU 分解的意義在於，將矩陣的「分解」與方程的「求解」分離。有什麼好處？
在不少應用場景中，當需要求解 $Ax=b$ 的時候，左邊的矩陣 $A$ 很多時候是不變的，而右邊的 $b$ 隨著輸入而變化。
做 $LU$ 分解時，只會用到矩陣 $A$ ，所以可以預先準備好 $L$ 與 $U$ ，當有求解 $b$ 的需求時，直接拿來用就好了:
$Ax=LUx=b$
（1）解方程 $Ly=b$ ，得到 $y$

（2）再通過 $Ux=y$ ，得到 $x$
將 $A$ 分解為 $LU$ 的時間複雜度是 $O(N^3)$ ，而求解 $LUx=b$ 只需要 $O(N^2)$ 。

LU分解無非就是將非三角矩陣A轉化為上三角矩陣U的過程的反向操作。

也就是將三角矩陣U還原回A的過程。

${displaystyle {egin{bmatrix}a_{11}&a_{12}\a_{21}&a_{22}\end{bmatrix}}={egin{bmatrix}l_{11}&0\l_{21}&l_{22}\end{bmatrix}}{egin{bmatrix}u_{11}&u_{12}\0&u_{22}\end{bmatrix}}}$

用2x2矩陣來簡單說明

按照我們平時消去的過程，用某一行消去下面的幾行的時候，原來那一行元素我們是不動的，立即就能得出L的對角元素為1，而U的第一行元素和A的第一行元素一樣。

${displaystyle {egin{bmatrix}a_{11}&a_{12}\a_{21}&a_{22}\end{bmatrix}}={egin{bmatrix}1&0\l_{21}&1\end{bmatrix}}{egin{bmatrix}a_{11}&a_{12}\0&a_{22}\end{bmatrix}}}$

而想把0變回原來的 $a_{21}$ 只需要在0上面加上 $u_{11} * (a_{21}/a_{11})$

顯而易見 $l_{21}$ 就是上式的 $(a_{21}/a_{11})$ 。

如果是3x3矩陣或更多 ${egin{bmatrix}a_{11}&a_{12}&a_{13}\a_{21}&a_{22}&a_{23}\a_{31}&a_{32}&a_{33}end{bmatrix}}$ ，可以確定U的第一行，以及L的對角線的1

${displaystyle {egin{bmatrix}a_{11}&a_{12}&a_{13}\a_{21}&a_{22}&a_{23}\a_{31}&a_{32}&a_{33}end{bmatrix}}={egin{bmatrix}1&0&0\l_{21}&1&0\l_{31}&l_{32}&1end{bmatrix}}{egin{bmatrix}a_{11}&a_{12}&a_{13}\0&u_{22}&u_{23}\0&0&u_{33}end{bmatrix}}}$

$l_{21}$ 和 $l_{31}$ 用最上面的方式求出後，對整個U的第二第三行進行操作，變為

${displaystyle {egin{bmatrix}a_{11}&a_{12}&a_{13}\0&a_{22}+a_{12}*l_{21}&a_{23}+a_{13}*l_{21}\0&a_{32}+a_{12}*l_{31}&a_{33}+a_{13}*l_{31}end{bmatrix}}={egin{bmatrix}1&0&0\0&1&0\0&l_{32}&1end{bmatrix}}{egin{bmatrix}a_{11}&a_{12}&a_{13}\0&u_{22}&u_{23}\0&0&u_{33}end{bmatrix}}}$

就變成了解

${displaystyle {egin{bmatrix}a_{22}+a_{12}*l_{21}&a_{23}+a_{13}*l_{21}\a_{32}+a_{12}*l_{31}&a_{33}+a_{13}*l_{31}end{bmatrix}}={egin{bmatrix}1&0\l_{32}&1end{bmatrix}}{egin{bmatrix}u_{22}&u_{23}\0&u_{33}end{bmatrix}}}$

$l_{32} =(a_{32}+ a_{12}*l_{31})/(a_{22}+a_{12}*l_{21})$

以此類推多階的矩陣的LU就可以求了。

在code上面LU放在一起計算就可以簡便很多。

以下C#代碼

double factor = 0; for (int k = 0; k < n-1; k++) { for (int i = k+1; i < n; i++) { factor =data[i, k] / data[k, k]; data[i, k] = factor; for (int j = k+1; j < n; j++) { data[i, j] =data[i, j] - factor * data[k, j]; } } }

裡面最開始data就是A，計算過後的data，上三角部分就是U，下三角部分把對角線換成1就是L。

而解矩陣「相除」就可以用到上面的LU矩陣。

對於Ax=b 我們可以先求LM=b求出M，而M又等於Ux就可以根據Ux=M來求出x

而逆矩陣就是一個特殊的Ax=b，即 $AA^{-1} = I = LUA^{-1}$

我們可以一列一列的來求A的逆矩陣，那麼Ax=b中x就是逆矩陣第n列，b就是單位矩陣的第n列。

int n = Convert.ToInt32(Math.Sqrt(data.Length)); double[,] minverse = new double[n, n]; double[,] UX = new double[n, n]; var U = Lmatrix(data); var L = Umatrix(data); for (int i = 0; i < n; i++) { var temp1 = forward(L, i); var x = backward(U, temp1); for (int j = 0; j < n; j++) { minverse[j, i] = x[j]; } }

//常見的PA=LU中的P是用來給A進行換行的，以免第一行一大堆零無法消去後面的幾行。

參考文獻：

LU Decomposition