矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解,LDLT(不带平方根乔累斯基)分解等,以及为了满足分解条件又加入行列变换的LPU分解,PLU分解,LUP分解,LDPU分解等。这里矩阵的三角分解系列教程主要是针对在学习三角分解时候的涉及到的一些细节,包括很多方法的来源和证明等,以及其中用到的一些矩阵操作的基础知识,主要包括:
- [矩阵的三角分解系列一] 高斯消元法
- [矩阵的三角分解系列二] LDU基本定理
- [矩阵的三角分解系列三] 杜利特/克劳特分解公式
- [矩阵的三角分解系列四] 乔累斯基(Cholesky)分解公式
- [矩阵的三角分解系列五] 三角分解中的行列变换
- [矩阵的三角分解系列六] Eigen中的三角分解
这个系列后面文章会用到前面文章的理论和技术,所以建议按照顺序查看。
简介
根据[矩阵的三角分解系列二] LDU基本定理文章介绍如果方阵
A
\boldsymbol{A}
A可分解成一个下三角矩阵
L
\boldsymbol{L}
L和一个上三角矩阵
U
\boldsymbol{U}
U的乘积,则称
A
\boldsymbol{A}
A可作三角分解或
L
U
\boldsymbol{LU}
LU分解。如果
L
\boldsymbol{L}
L是单位下三角矩阵,
U
\boldsymbol{U}
U为上三角矩阵,此时的三角分解称为杜利特(Doolittle)分解;如果
L
\boldsymbol{L}
L是下三角矩阵,
U
\boldsymbol{U}
U为单位上三角矩阵,此时的三角分解称为克劳特(Crout)分解。
但之前讲解的都是按照高斯消元过程中介绍的高斯消元过程进行
L
U
\boldsymbol{L U}
LU,
L
D
U
\boldsymbol{LDU}
LDU或者杜利特/克劳特分解的。但是实际使用过程如果矩阵
A
\boldsymbol{A}
A的阶数特别高,使用这种方式分解都会比较麻烦。简单描述下就是需要一步步求出每次的消元矩阵
L
(
k
)
\boldsymbol{L}^{(k)}
L(k),然后和上次计算的
A
(
k
−
1
)
\boldsymbol{A}^{(k-1)}
A(k−1)相乘最后得到
U
\boldsymbol{U}
U,然后所有的
L
(
k
)
\boldsymbol{L}^{(k)}
L(k)的逆相乘得到
L
\boldsymbol{L}
L。不管是求
L
\boldsymbol{L}
L和
U
\boldsymbol{U}
U矩阵都不是那么直接。所以对于具体的杜利特(Doolittle)分解和克劳特(Crout)分解,可以有更直接的方法来进行计算分解矩阵。
杜利特/克劳特分解
A
\boldsymbol{A}
A为一般的
n
n
n阶方阵,有分解式(这里以克劳特分解举例说明,杜利特分解公式类似)
A
=
L
U
,
(1)
\boldsymbol{A} = \boldsymbol{L U}, \tag{1}
A=LU,(1)
即
[
a
11
⋯
a
1
j
⋯
a
1
n
⋮
⋮
⋮
a
i
1
⋯
a
i
j
⋯
a
i
n
⋮
⋮
⋮
a
n
1
⋯
a
n
j
⋯
a
n
n
]
=
[
l
11
⋮
⋱
l
i
1
⋯
l
i
i
⋮
⋱
l
n
1
⋯
l
n
i
⋯
l
n
n
]
[
1
u
12
⋯
u
1
j
⋯
u
1
n
⋱
⋱
⋮
⋮
1
u
j
−
1
,
j
⋯
u
j
−
1
,
n
⋱
⋱
⋮
1
u
n
−
1
,
n
1
]
,
\left[\begin{matrix} a_{11} & \cdots & a_{1 j} & \cdots & a_{1 n} \\ \vdots & & \vdots & & \vdots \\ a_{i 1} & \cdots & a_{i j} & \cdots & a_{i n} \\ \vdots & & \vdots & & \vdots \\ a_{n 1} & \cdots & a_{n j} & \cdots & a_{n n} \end{matrix}\right] = \left[\begin{matrix} l_{11} & & & & \\ \vdots & \ddots & & & \\ l_{i 1} & \cdots & l_{i i} & & \\ \vdots & & & \ddots & \\ l_{n 1} & \cdots & l_{n i} & \cdots & l_{n n} \end{matrix}\right] \left[\begin{matrix} 1 & u_{12} & \cdots & u_{1 j} & \cdots & u_{1 n} \\ & \ddots & \ddots & \vdots & & \vdots \\ & & 1 & u_{j-1, j} & \cdots & u_{j-1, n} \\ & & & \ddots & \ddots & \vdots \\ & & & & 1 & u_{n-1, n} \\ & & & & & 1 \end{matrix}\right],
⎣⎢⎢⎢⎢⎢⎢⎡a11⋮ai1⋮an1⋯⋯⋯a1j⋮aij⋮anj⋯⋯⋯a1n⋮ain⋮ann⎦⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡l11⋮li1⋮ln1⋱⋯⋯liilni⋱⋯lnn⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡1u12⋱⋯⋱1u1j⋮uj−1,j⋱⋯⋯⋱1u1n⋮uj−1,n⋮un−1,n1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤,
目标就是直接求出上面下三角矩阵
L
\boldsymbol{L}
L的每个元素
l
i
j
l_{ij}
lij以及单位上三角矩阵
U
\boldsymbol{U}
U的每个元素
u
i
j
u_{ij}
uij。这里考虑
A
\boldsymbol{A}
A的具体第
i
i
i行第
j
j
j列的元素的计算:
a
i
j
=
(
l
i
1
,
⋯
,
l
i
i
,
0
,
⋯
,
0
)
[
u
1
j
⋮
u
j
−
1
,
j
1
0
⋮
0
]
,
a_{i j}=\left(l_{i 1}, \cdots, l_{i i}, 0, \cdots, 0\right)\left[\begin{matrix} u_{1 j} \\ \vdots \\ u_{j-1, j} \\ 1 \\ 0 \\ \vdots \\ 0 \end{matrix}\right],
aij=(li1,⋯,lii,0,⋯,0)⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡u1j⋮uj−1,j10⋮0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,
当
i
≥
j
i \geq j
i≥j时,也就是计算
n
n
n阶方阵
A
\boldsymbol{A}
A的下三角部分(包括对角线元素),所以
a
i
j
=
∑
k
=
1
j
l
i
k
u
k
j
=
∑
k
=
1
j
−
1
l
i
k
u
k
j
+
l
i
j
,
a_{i j}=\sum_{k=1}^{j} l_{i k} u_{k j}=\sum_{k=1}^{j-1} l_{i k} u_{k j}+l_{i j},
aij=k=1∑jlikukj=k=1∑j−1likukj+lij,
得
l
i
j
=
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
u
k
j
,
i
=
1
,
⋯
,
n
,
j
=
1
,
⋯
,
i
,
(2)
l_{i j}=a_{i j}-\sum_{k=1}^{j-1} l_{i k} u_{k j}, \quad i=1, \cdots, n, \quad j=1, \cdots, i, \tag{2}
lij=aij−k=1∑j−1likukj,i=1,⋯,n,j=1,⋯,i,(2)
又当
i
<
j
i < j
i<j时,也就是计算
n
n
n阶方阵
A
\boldsymbol{A}
A的下三角部分,所以
a
i
j
=
∑
k
=
1
i
l
i
k
u
k
j
=
∑
k
=
1
i
−
1
l
i
k
u
k
j
+
l
i
i
u
i
j
,
a_{i j}=\sum_{k=1}^{i} l_{i k} u_{k j}=\sum_{k=1}^{i-1} l_{i k} u_{k j}+l_{i i} u_{i j},
aij=k=1∑ilikukj=k=1∑i−1likukj+liiuij,
得
u
i
j
=
(
a
i
j
−
∑
k
=
1
l
i
k
u
k
j
)
/
l
i
i
,
i
=
1
,
⋯
,
n
−
1
,
j
=
i
+
1
,
⋯
,
n
,
(3)
u_{i j}=\left(a_{i j}-\sum_{k=1}^{} l_{i k} u_{k j}\right) / l_{i i}, \quad i=1, \cdots, n-1, \quad j=i+1, \cdots, n, \tag{3}
uij=(aij−k=1∑likukj)/lii,i=1,⋯,n−1,j=i+1,⋯,n,(3)
公式
(
9
)
(9)
(9)和公式
(
10
)
(10)
(10)就被称为克劳特分解公式。
看上面分解公式比较复杂,而且是非线性的,看起来用这些公式求解
L
\boldsymbol{L}
L和
U
\boldsymbol{U}
U的元素比较困难。但如果注意计算的顺序“自上而下且从左到右”一行行的去计算的话,当前元素依赖的值已经是已知的了,这样就能顺利的求出
L
\boldsymbol{L}
L和
U
\boldsymbol{U}
U的所有的元素。
所以利用克劳特分解公式进行三角分解的计算步骤为:
第
一
步
,
计
算
l
11
u
12
u
13
⋯
u
1
n
;
第
二
步
,
计
算
l
21
l
22
u
23
⋯
u
2
n
;
⋮
第
n
步
,
计
算
l
n
1
l
n
2
−
l
n
3
⋯
l
n
n
.
\begin{aligned} 第一步,计算 & \quad l_{11} \quad u_{12} \quad u_{13} \quad \cdots \quad u_{1 n}; \\ 第二步,计算 & \quad l_{21} \quad l_{22} \quad u_{23} \quad \cdots \quad u_{2 n}; \\ &\vdots \\ 第n步,计算 & \quad l_{n 1} \quad l_{n 2} \quad-l_{n 3} \quad \cdots \quad l_{n n}. \end{aligned}
第一步,计算第二步,计算第n步,计算l11u12u13⋯u1n;l21l22u23⋯u2n;⋮ln1ln2−ln3⋯lnn.
类似地,可得到
n
n
n阶方阵
A
\boldsymbol{A}
A的杜利特分解的计算公式:
u
i
j
=
a
i
j
−
∑
k
=
1
i
−
1
l
i
k
u
k
j
,
i
=
1
,
⋯
,
n
,
j
=
i
,
⋯
,
n
,
l
i
j
=
(
a
i
j
−
∑
k
=
1
′
l
i
k
u
k
j
)
/
u
j
j
,
i
=
2
,
⋯
,
n
,
j
=
1
,
⋯
,
i
−
1.
(11)
\begin{aligned} u_{i j}&=a_{i j}-\sum_{k=1}^{i-1} l_{i k} u_{k j}, \quad &i=1, \cdots, n, \quad j=i, \cdots, n, \\ l_{i j}&=\left(a_{i j}-\sum_{k=1}^{\prime} l_{i k} u_{k j}\right) / u_{j j}, \quad &i=2, \cdots, n, \quad j=1, \cdots, i-1. \tag{11} \end{aligned}
uijlij=aij−k=1∑i−1likukj,=(aij−k=1∑′likukj)/ujj,i=1,⋯,n,j=i,⋯,n,i=2,⋯,n,j=1,⋯,i−1.(11)
例子
将下列矩阵进行克劳特分解
A
=
[
4
8
4
2
7
2
1
2
3
]
.
\boldsymbol{A} = \left[ \begin{matrix} 4 & 8 & 4 \\ 2 & 7 & 2 \\ 1 & 2 & 3 \end{matrix} \right].
A=⎣⎡421872423⎦⎤.
解:
根据公式
(
9
)
(9)
(9)和公式
(
10
)
(10)
(10)可得
l
11
=
a
11
=
4
u
12
=
a
12
/
l
11
=
2
u
13
=
a
13
/
l
11
=
1
l
21
=
a
21
=
2
l
22
=
a
22
−
l
21
u
12
=
3
u
23
=
(
a
23
−
l
21
u
13
)
/
l
22
=
0
l
31
=
a
31
=
1
l
32
=
a
32
−
l
31
u
12
=
0
l
33
=
a
33
−
l
31
u
13
−
l
32
u
23
=
2
\begin{aligned} l_{11} &=a_{11}=4 \\ u_{12} &=a_{12} / l_{11}=2 \\ u_{13} &=a_{13} / l_{11}=1 \\ l_{21} &=a_{21}=2 \\ l_{22} &=a_{22}-l_{21} u_{12}=3 \\ u_{23} &=\left(a_{23}-l_{21} u_{13}\right) / l_{22}=0 \\ l_{31} &=a_{31}=1 \\ l_{32} &=a_{32}-l_{31} u_{12}=0 \\ l_{33} &=a_{33}-l_{31} u_{13}-l_{32} u_{23}=2 \end{aligned}
l11u12u13l21l22u23l31l32l33=a11=4=a12/l11=2=a13/l11=1=a21=2=a22−l21u12=3=(a23−l21u13)/l22=0=a31=1=a32−l31u12=0=a33−l31u13−l32u23=2
则
A
=
L
U
=
[
4
2
3
1
0
2
]
[
1
2
1
1
0
1
]
\boldsymbol{A} = \boldsymbol{L U} = \left[ \begin{matrix} 4 & \\ 2 & 3 \\ 1 & 0 & 2 \end{matrix} \right] \left[ \begin{matrix} 1 & 2 & 1 \\ & 1 & 0 \\ & & 1 \end{matrix} \right]
A=LU=⎣⎡421302⎦⎤⎣⎡121101⎦⎤
引用
【1】 矩阵论(第二版)