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}}$

先將矩陣第一列元素中 a11 以下的所有元素變為 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}}$

再將矩陣第二列元素中 a22 以下的所有元素變為 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_2 級別，較高斯法為優。

### 計算行列式

矩陣 $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.