LU分解
LU分解是矩阵分解的一种，将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，有时需要再乘上一个置换矩阵。
LU分解可以被视为高斯消元法的矩阵形式。在数值计算上，LU分解经常被用来解线性方程组、且在求逆矩阵和计算行列式中都是一个关键的步骤。
一、定义
对于方阵
A
A
A，
A
A
A 的LU分解是将它分解成一个下三角矩阵 L 与上三角矩阵 U 的乘积，也就是
A
=
L
U
A=LU
A=LU。
举例来说一个
3
×
3
{\displaystyle 3\times 3}
3×3的矩阵
A
A
A ，其 LU 分解会写成下面的形式：
A
=
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
=
[
l
11
0
0
l
21
l
22
0
l
31
l
32
l
33
]
[
u
11
u
12
u
13
0
u
22
u
23
0
0
u
33
]
{\displaystyle 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}}}
A=⎣⎡a11a21a31a12a22a32a13a23a33⎦⎤=⎣⎡l11l21l310l22l3200l33⎦⎤⎣⎡u1100u12u220u13u23u33⎦⎤
LDU 分解
方阵 A 的LDU分解是是将它分解成一个单位下三角矩阵 L、对角矩阵 D 与单位上三角矩阵 U 的乘积，即
A
=
L
D
U
{\displaystyle A=LDU}
A=LDU
其中单位上、下三角矩阵是指对角线上全是 1 的上、下三角矩阵。事实上，LDU 分解可以推广到 A 是一般的矩阵，而非方阵。
存在性和唯一性
- 一个可逆矩阵可以进行LU分解当且仅当它的所有子式都非零。
- 如果要求其中的L矩阵（或U矩阵）为单位三角矩阵，那么分解是唯一的。
- 即使矩阵不可逆，LU仍然可能存在。实际上，如果一个秩为k的矩阵的前k个顺序主子式不为零，那么它就可以进行LU分解，但反之则不然。
二、算法
LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法（Doolittle algorithm），从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。
杜尔里特算法
对给定的N × N矩阵
A
=
(
a
n
,
n
)
{\displaystyle A=(a_{n,n})}
A=(an,n)有
A
(
0
)
:
=
A
A^{{(0)}}:=A
A(0):=A
然后定义对于
n
=
1
,
.
.
.
,
N
−
1
n = 1,...,N-1
n=1,...,N−1 的情况如下：
在第n步，消去矩阵
A
(
n
−
1
)
A^{(n-1)}
A(n−1)的第n列主对角线下的元素：将
A
(
n
−
1
)
A^{(n-1)}
A(n−1)的第n行乘以
l
i
,
n
=
−
a
i
,
n
(
n
−
1
)
a
n
,
n
(
n
−
1
)
{\displaystyle l_{i,n}=-{\frac {a_{i,n}^{(n-1)}}{a_{n,n}^{(n-1)}}}}
li,n=−an,n(n−1)ai,n(n−1)之后加到第
i
i
i行上去。其中
i
=
n
+
1
,
…
,
N
{\displaystyle i=n+1,\ldots ,N}
i=n+1,…,N。这相当于在
A
(
n
−
1
)
A^{(n-1)}
A(n−1)的左边乘上一个单位下三角矩阵：
L
n
=
[
1
0
⋱
1
l
n
+
1
,
n
⋱
⋮
⋱
0
l
N
,
n
1
]
{\displaystyle L_{n}={\begin{bmatrix}1&&&&&0\\&\ddots &&&&\\&&1&&&\\&&l_{n+1,n}&\ddots &&\\&&\vdots &&\ddots &\\0&&l_{N,n}&&&1\\\end{bmatrix}}}
Ln=⎣⎢⎢⎢⎢⎢⎢⎢⎡10⋱1ln+1,n⋮lN,n⋱⋱01⎦⎥⎥⎥⎥⎥⎥⎥⎤
于是，设
A
(
n
)
=
L
n
A
(
n
−
1
)
{\displaystyle A^{(n)}=L_{n}A^{(n-1)}}
A(n)=LnA(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
)
=
…
=
L
1
−
1
…
L
N
−
1
−
1
A
(
N
−
1
)
{\displaystyle 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=L1−1L1A(0)=L1−1A(1)=L1−1L2−1L2A(1)=L1−1L2−1A(2)=…=L1−1…LN−1−1A(N−1)
这时，矩阵
A
(
N
−
1
)
A^{(N-1)}
A(N−1) 就是U，
L
=
L
1
−
1
…
L
N
−
1
−
1
{\displaystyle L=L_{1}^{-1}\ldots L_{N-1}^{-1}}
L=L1−1…LN−1−1。 下三角矩阵
L
k
{\displaystyle L_{k}}
Lk的逆依然是下三角矩阵，而且下三角矩阵的乘积仍是下三角矩阵，所以
L
{\displaystyle L}
L是下三角矩阵。 于是我们得到分解：
A
=
L
U
{\displaystyle A=LU}
A=LU。
显然，要是算法成立，在每步操作时必须有
a
n
,
n
(
n
−
1
)
≠
0
{\displaystyle a_{n,n}^{(n-1)}\neq 0}
an,n(n−1)=0。如果这一条件不成立，就要将第n行和另一行交换，由此就会出现一个置换矩阵P。这就是为什么一般来说LU分解里会带有一个置换矩阵的原因。
例子：
将一个简单的3×3矩阵A进行LU分解：
A
=
[
1
2
3
2
5
7
3
5
3
]
{\displaystyle A={\begin{bmatrix}1&2&3\\2&5&7\\3&5&3\\\end{bmatrix}}}
A=⎣⎡123255373⎦⎤
先将矩阵第一列元素中
a
11
a_{11}
a11以下的所有元素变为0，即
L
1
A
=
[
1
0
0
−
2
1
0
−
3
0
1
]
×
[
1
2
3
2
5
7
3
5
3
]
=
[
1
2
3
0
1
1
0
−
1
−
6
]
{\displaystyle 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}}}
L1A=⎣⎡1−2−3010001⎦⎤×⎣⎡123255373⎦⎤=⎣⎡10021−131−6⎦⎤
再将矩阵第二列元素中
a
22
a_{22}
a22以下的所有元素变为0，即
L
2
(
L
1
A
)
=
[
1
0
0
0
1
0
0
1
1
]
×
[
1
2
3
0
1
1
0
−
1
−
6
]
=
[
1
2
3
0
1
1
0
0
−
5
]
=
U
{\displaystyle 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}
L2(L1A)=⎣⎡100011001⎦⎤×⎣⎡10021−131−6⎦⎤=⎣⎡10021031−5⎦⎤=U
L
=
L
1
−
1
L
2
−
1
=
[
1
0
0
2
1
0
3
0
1
]
×
[
1
0
0
0
1
0
0
−
1
1
]
=
[
1
0
0
2
1
0
3
−
1
1
]
{\displaystyle 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}}}
L=L1−1L2−1=⎣⎡123010001⎦⎤×⎣⎡10001−1001⎦⎤=⎣⎡12301−1001⎦⎤
