LU分解

在线性代数中， LU分解是矩阵分解的一种，可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积（有时是它们和一个置换矩阵的乘积）。LU分解主要应用在数值分析中，用来解线性方程、求反矩陣或计算行列式。

定义

设A是一个方块矩阵。A的LU分解是将它分解成如下形式：

$A = LU, \,$

其中L和U分别是下三角矩阵和上三角矩阵。

例如对于一个$3 \times 3$的矩阵，就有

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

一个LDU分解是一个如下形式的分解：

$A = LDU$

其中D是对角矩阵，L和U是单位三角矩阵（对角线上全是1的三角矩阵）。

一个LUP分解是一个如下形式的分解：

$A = LUP$

其中L和U仍是三角矩阵，P是一个置换矩阵。

一个充分消元的LU分解为如下形式：

$PAQ = LU$

存在性和唯一性

一个可逆矩阵可以进行LU分解当且仅当它的所有子式都非零。如果要求其中的L矩阵（或U矩阵）为单位三角矩阵，那么分解是唯一的。同理可知，矩阵的LDU可分解条件也相同，并且总是唯一的。

即使矩阵不可逆，LU仍然可能存在。实际上，如果一个秩为k的矩阵的前k个顺序主子式不为零，那么它就可以进行LU分解，但反之则不然。

目前，在任意域上一个方块矩阵可进行LU分解的充要条件已经被发现，这些充要条件可以用某些特定子矩阵的秩表示。用高斯消元法来得到LU分解的算法也可以扩张到任意域上。

任意矩阵A（不仅仅是方块矩阵）都可以进行LUP分解。其中的L和P矩阵是方阵，U矩阵则与A形状一样。

正定矩阵

如果矩阵A是埃尔米特矩阵，并且是正定矩阵，那么可以使，U是L的共轭转置。也就是说，A可以写成

$A = L L^{*}. \,$

这个分解被称作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{pmatrix} 1 & & & & & 0 \\ & \ddots & & & & \\ & & 1 & & & \\ & & l_{n+1,n} & \ddots & & \\ & & \vdots & & \ddots & \\ 0 & & l_{N,n} & & & 1 \\ \end{pmatrix}.$

于是，定义为：设

$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}$。于是我们得到分解：$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也是稀疏矩阵。理论上来说，算法的复杂度约等于非零系数的个数，而不是矩阵的大小阶数。这些算法通过运用行和列的交换，使得过程中零系数因为操作而变成非零系数的次数减到最少。

一般的将零系数因为操作而变成非零系数的次数减到最少的方法是运用图论。

应用

求解线性方程

对于给定的线性方程组

$A x = L U x = 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.