LU分解

在線性代數與數值分析中，LU分解是矩陣分解的一種，將一個矩陣分解為一個下三角矩陣和一個上三角矩陣的乘積，有時需要再乘上一個置換矩陣。LU分解可以被視為高斯消去法的矩陣形式。在數值計算上，LU分解經常被用來解線性方程組、且在求逆矩陣和計算行列式中都是一個關鍵的步驟。

定義

對於方陣 $A$， $A$ 的 LU 分解是將它分解成一個下三角矩陣 L 與上三角矩陣 U 的乘積，也就是

$A=LU$

如果適當的改變 $A$的行的順序或列的順序，就可以將 $A$做 LU 分解。

舉例來說一個 $3\times 3$ 的矩陣 A ，其 LU 分解會寫成下面的形式：

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

事實上，並不是每個矩陣都有 LU 分解。例如，從上式可知 $a_{11}=l_{11}u_{11}$，若 $a_{11}=0$，則 $l_{11}$或 $u_{11}$等於 0，故 L 或 U 是非可逆矩陣，A 必須也是非可逆矩陣。然而，存在著可逆矩陣 A 滿足 $a_{11}=0$，這些 A 就是沒有 LU 分解的例子。該問題可藉由置換 A 的各列順序來解決，最終會得到一個 A 的 PLU 分解。

PLU 分解

方陣 A 的 PLU 分解是是將它分解成一個置換矩陣 P、一個下三角矩陣 L 與上三角矩陣 U 的乘積，即

$A=PLU$

事實上，所有的方陣都可以寫成 PLU 分解的形式，由於左乘置換矩陣 $P^{-1}$是在交換行的順序，所以由 $P^{-1}A=LU$推得適當的交換 A 的行的順序，即可將 A 做 LU 分解。事實上，PLU 分解有很高的數值穩定性，因此實用上是很好用的工具。

有時為了計算上的方便，會同時間換行與列的順序，此時會將 A 分解成

$A=PLUQ$

其中 P、L、U 同上，Q 是一個置換矩陣。

LDU 分解

方陣 A 的 LDU 分解是是將它分解成一個單位下三角矩陣 L、對角矩陣 D 與單位上三角矩陣 U 的乘積，即

$A=LDU$

其中單位上、下三角矩陣是指對角線上全是 1 的上、下三角矩陣。

事實上，LDU 分解可以推廣到 A 是一般的矩陣，而非方陣。此時，L 和 D 是方陣，並且與 A 有相同的列，U 則有和 A 相同的長寬。注意到現在 U 是上三角的定義改為主對角線的下方都是 0，而主對角線是收集所有 $U_{ij}$ 滿足 i=j。

存在性和唯一性

一個可逆矩陣可以進行LU分解若且唯若它的所有子式都非零。如果要求其中的L矩陣（或U矩陣）為單位三角矩陣，那麼分解是唯一的。同理可知，矩陣的LDU可分解條件也相同，並且總是唯一的。

即使矩陣不可逆，LU仍然可能存在。實際上，如果一個秩為k的矩陣的前k個順序主子式不為零，那麼它就可以進行LU分解，但反之則不然。

目前，在任意域上一個方塊矩陣可進行LU分解的充要條件已經被發現，這些充要條件可以用某些特定子矩陣的秩表示。用高斯消去法來得到LU分解的算法也可以擴張到任意域上。

任意矩陣A（不僅僅是方塊矩陣）都可以進行LUP分解。其中的L和U矩陣是方陣，P矩陣則與A形狀一樣。

正定矩陣

如果矩陣A是埃爾米特矩陣，並且是正定矩陣，那麼可以使，U是L的共軛轉置。也就是說，A可以寫成

$A=LL^{*}\$

這個分解被稱作Cholesky分解。對每一個正定矩陣，Cholesky分解都唯一存在。此外，比起一般的LU分解，計算Cholesky分解更為快捷，並具有更高的數值穩定性。

具體的表達式

由於LDU分解唯一存在，對給定的矩陣，可以給出相應三個矩陣L、D和U的具體的表達式。表達式由A的主子式之比構成（因此要求它們不為零）。設$d_{1},d_{2},\cdots d_{n}$為矩陣D的對角線係數，則有$d_{1}=\mathbf {A} _{1,1}$。對$i=2,\ldots ,n$，$d_{i}$的值等於A的第$i$個順序主子式與第$i-1$個順序主子式之比，其中約定$d_{0}$=1。

算法

LU分解在本質上是高斯消去法的一種表達形式。實質上是將A通過初等行轉換變成一個上三角矩陣，其轉換矩陣就是一個單位下三角矩陣。這正是所謂的杜爾里特算法（Doolittle algorithm）：從下至上地對矩陣A做初等行轉換，將對角線左下方的元素變成零，然後再證明這些行轉換的效果等同於左乘一系列單位下三角矩陣，這一系列單位下三角矩陣的乘積的逆就是L矩陣，它也是一個單位下三角矩陣。

這類算法的複雜度一般在${\frac {2n^{3}}{3}}$左右，對充分消去的分解則不然。

杜爾里特算法

對給定的N × N矩陣

$A=(a_{n,n})$

有

$A^{(0)}:=A$

然後定義對於n = 1,...,N-1的情況如下：

在第n步，消去矩陣A^(n-1)的第n列主對角線下的元素：將A^(n-1)的第n行乘以$l_{i,n}:=-{\frac {a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}}$之後加到第i行上去。其中$i=n+1,\ldots ,N$。

這相當於在A^(n-1)的左邊乘上一個單位下三角矩陣：

$L_{n}={\begin{bmatrix}1&&&&&0\\&\ddots &&&&\\&&1&&&\\&&l_{n+1,n}&\ddots &&\\&&\vdots &&\ddots &\\0&&l_{N,n}&&&1\\\end{bmatrix}}$

於是，定義為：設

$A^{(n)}:=L_{n}A^{(n-1)}$

經過N-1輪操作後，所有在主對角線下的係數都為0了，於是我們得到了一個上三角矩陣：A^(N-1)，這時就有：

$A=L_{1}^{-1}L_{1}A^{(0)}=L_{1}^{-1}A^{(1)}=L_{1}^{-1}L_{2}^{-1}L_{2}A^{(1)}=L_{1}^{-1}L_{2}^{-1}A^{(2)}=\ldots =L_{1}^{-1}\ldots L_{N-1}^{-1}A^{(N-1)}$

這時，矩陣A^(N-1) 就是U，$L=L_{1}^{-1}\ldots L_{N-1}^{-1}$。 下三角矩陣$L_{k}$的逆依然是下三角矩陣，而且下三角矩陣的乘積仍是下三角矩陣，所以$L$是下三角矩陣。 於是我們得到分解：$A=LU$。

顯然，要是算法成立，在每步操作時必須有$a_{n,n}^{(n-1)}\neq 0$。如果這一條件不成立，就要將第n行和另一行交換，由此就會出現一個置換矩陣P。這就是為什麼一般來說LU分解里會帶有一個置換矩陣的原因。

例子

將一個簡單的3×3矩陣A進行LU分解：

$A={\begin{bmatrix}1&2&3\\2&5&7\\3&5&3\\\end{bmatrix}}$

先將矩陣第一列元素中a₁₁以下的所有元素變為0，即

$L_{1}A={\begin{bmatrix}1&0&0\\-2&1&0\\-3&0&1\\\end{bmatrix}}\times {\begin{bmatrix}1&2&3\\2&5&7\\3&5&3\\\end{bmatrix}}={\begin{bmatrix}1&2&3\\0&1&1\\0&-1&-6\\\end{bmatrix}}$

再將矩陣第二列元素中a₂₂以下的所有元素變為0，即

$L_{2}(L_{1}A)={\begin{bmatrix}1&0&0\\0&1&0\\0&1&1\\\end{bmatrix}}\times {\begin{bmatrix}1&2&3\\0&1&1\\0&-1&-6\\\end{bmatrix}}={\begin{bmatrix}1&2&3\\0&1&1\\0&0&-5\\\end{bmatrix}}=U$

$L=L_{1}^{-1}L_{2}^{-1}={\begin{bmatrix}1&0&0\\2&1&0\\3&0&1\\\end{bmatrix}}\times {\begin{bmatrix}1&0&0\\0&1&0\\0&-1&1\\\end{bmatrix}}={\begin{bmatrix}1&0&0\\2&1&0\\3&-1&1\\\end{bmatrix}}$

還有一種方法是通過方程式求解，如下所示，將以下矩陣進行LU分解：

${\begin{bmatrix}4&3\\6&3\\\end{bmatrix}}={\begin{bmatrix}l_{11}&0\\l_{21}&l_{22}\\\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}\\0&u_{22}\\\end{bmatrix}}$

由於矩陣階數只是2，可以直接列方程式解：

$l_{11}*u_{11}+0*0=4$

$l_{11}*u_{12}+0*u_{22}=3$

$l_{21}*u_{11}+l_{22}*0=6$

$l_{21}*u_{12}+l_{22}*u_{22}=3$

這個線性方程組有無數多組解。因此，可以假設其中一個是單位三角矩陣，比如說L，也就是說其對角線上的兩個係數都是1。這時可以解出：

$l_{21}=1.5$

$u_{11}=4$

$u_{12}=3$

$u_{22}=-1.5$。

也就是說

${\begin{bmatrix}4&3\\6&3\\\end{bmatrix}}={\begin{bmatrix}1&0\\1.5&1\\\end{bmatrix}}{\begin{bmatrix}4&3\\0&-1.5\\\end{bmatrix}}$

稀疏矩陣分解

對於階數很大的稀疏矩陣，有特別的簡便算法來獲得其LU分解：這時的L和U也是稀疏矩陣。理論上來說，算法的複雜度約等於非零係數的個數，而不是矩陣的大小階數。這些算法通過運用行和列的交換，使得過程中零係數因為操作而變成非零係數的次數減到最少。

一般的將零係數因為操作而變成非零係數的次數減到最少的方法是運用圖論。

應用

求解線性方程式

對於給定的線性方程組

$Ax=LUx=b\,$

要解出x，可以進行以下步驟：

1. 首先，解方程式$Ly=b$得到$y$；
2. 然後解方程式$Ux=y$ 得到$x$。

在兩次的求解中，我們遇到的都是三角矩陣，因此運用向前（向後）替代法就可以簡潔地求解（參見三角矩陣），而不需要用到高斯消去法。然而，在將A進行LU分解時，仍然要用到高斯消去法。因此，這個方法適合在要對許多個不同的b求解時用。

求逆矩陣

求矩陣A的逆時，可以直接求L和U的逆矩陣，然後代入：$A^{-1}=U^{-1}L^{-1}$ 。也可以將單位矩陣分解成n個列向量，然後用上面求解線性方程式的方法解出逆矩陣的列向量，然後拼起來。後者的複雜度在n²級別，較高斯法為優。

計算行列式

矩陣$L$和$U$可以用來快速地計算矩陣$A$的行列式，因為det(A) = det(L) det(U)，而三角矩陣的行列式就是對角線元素的乘積。如果要求L 是單位三角矩陣，那麼$\det(A)=\det(L)\det(U)=\prod _{i=1}^{n}u_{ii}$

同樣的方法也可以應用於LUP分解，只需乘上P的行列式，即相應置換的符號差。

參見

• 分塊LU分解
• Cholesky分解
• 矩陣分解
• LU約簡

參考來源

• Trefethen, Lloyd; Bau, David, Numerical Linear Algebra 
• Cormen, T.H.; Leisserson, C.E; Rivest, R.L., Introduction to Algorithms 
• Golub, Gene H.; Van Loan, Charles F., Matrix Computations 3rd, Baltimore: Johns Hopkins, 1996, ISBN 978-0-8018-5414-9 .
• Horn, Roger A.; Johnson, Charles R., Matrix Analysis, Cambridge University Press, 1985, ISBN 0-521-38632-2 . See Section 3.5.
• Okunev, Pavel; Johnson, Charles, Necessary And Sufficient Conditions For Existence of the LU Factorization of an Arbitrary Matrix, 1997, arXiv:math.NA/0506382 .
• Householder, Alston, The Theory of Matrices in Numerical Analysis, 1975 .
• LU decomposition （頁面存檔備份，存於網際網路檔案館） on MathWorld.
• LU decomposition （頁面存檔備份，存於網際網路檔案館） on Math-Linux.
• LU decomposition （頁面存檔備份，存於網際網路檔案館）

外部連結

• LAPACK （頁面存檔備份，存於網際網路檔案館） is a collection of FORTRAN subroutines for solving dense linear algebra problems
• ALGLIB （頁面存檔備份，存於網際網路檔案館） includes a partial port of the LAPACK to C++, C#, Delphi, etc.
• Online Matrix Calculator performs LU decomposition
• LU decomposition （頁面存檔備份，存於網際網路檔案館） at Holistic Numerical Methods Institute
• Module for LU Factorization with Pivoting
• LU Decomposition （頁面存檔備份，存於網際網路檔案館） by Ed Pegg, Jr.，The Wolfram Demonstrations Project，2007.