矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解,LDLT(不带平方根乔累斯基)分解等,以及为了满足分解条件又加入行列变换的LPU分解,PLU分解,LUP分解,LDPU分解等。这里矩阵的三角分解系列教程主要是针对在学习三角分解时候的涉及到的一些细节,包括很多方法的来源和证明等,以及其中用到的一些矩阵操作的基础知识,主要包括:
- [矩阵的三角分解系列一] 高斯消元法
- [矩阵的三角分解系列二] LDU基本定理
- [矩阵的三角分解系列三] 杜利特/克劳特分解公式
- [矩阵的三角分解系列四] 乔累斯基(Cholesky)分解公式
- [矩阵的三角分解系列五] 三角分解中的行列变换
- [矩阵的三角分解系列六] Eigen中的三角分解
这个系列后面文章会用到前面文章的理论和技术,所以建议按照顺序查看。
简介
矩阵的三角分解依赖的具体方法是之前讲到的高斯消元过程,所以看本文前建议先把高斯消元过程先看完。
按照高斯消元过程文章中假设矩阵
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1个顺序主子式都不为零,那么矩阵
A
\boldsymbol{A}
A高斯消元过程能够进行到底。根据高斯消元过程令
U
=
A
(
n
)
=
L
(
n
−
1
)
⋯
L
(
2
)
L
(
1
)
A
.
(1)
\boldsymbol{U} = \boldsymbol{A^{(n)}} = \boldsymbol{L^{(n-1)}}\cdots\boldsymbol{L^{(2)}}\boldsymbol{L^{(1)}}\boldsymbol{A}. \tag{1}
U=A(n)=L(n−1)⋯L(2)L(1)A.(1)
又有
d
e
t
L
(
k
)
=
1
≠
0
(
k
=
1
,
2
,
⋯
,
n
−
1
)
det\boldsymbol{L^{(k)}} = 1 \not = 0(k=1,2,\cdots,n-1)
detL(k)=1=0(k=1,2,⋯,n−1),所以
n
−
1
n-1
n−1个消元矩阵
L
(
k
)
\boldsymbol{L^{(k)}}
L(k)是可逆的,于是有
A
=
(
L
(
1
)
)
−
1
(
L
(
2
)
)
−
1
⋯
(
L
(
n
−
1
)
)
−
1
U
\boldsymbol{A} = (\boldsymbol{L^{(1)}})^{-1}(\boldsymbol{L^{(2)}})^{-1}\cdots(\boldsymbol{L^{(n-1)}})^{-1}\boldsymbol{U}
A=(L(1))−1(L(2))−1⋯(L(n−1))−1U
根据逆矩阵的定义知
L
(
k
)
=
[
1
⋱
1
l
k
+
1
,
k
1
⋮
⋱
l
n
k
1
]
,
k
=
1
,
2
,
⋯
,
n
−
1.
\boldsymbol{L^{(k)}} = \left[ \begin{matrix} 1 & \\ & \ddots \\ & & 1 \\ & & l_{k+1, k} & 1 \\ & & \vdots & & \ddots \\ & & l_{nk} & & & 1 \end{matrix} \right], k = 1, 2,\cdots,n-1.
L(k)=⎣⎢⎢⎢⎢⎢⎢⎢⎡1⋱1lk+1,k⋮lnk1⋱1⎦⎥⎥⎥⎥⎥⎥⎥⎤,k=1,2,⋯,n−1.
相当于负号去掉。
而且都是下三角矩阵,连乘仍是下三角矩阵,则
L
=
(
L
(
1
)
)
−
1
(
L
(
2
)
)
−
1
⋯
(
L
(
n
−
1
)
)
−
1
=
[
1
l
21
1
l
31
l
32
1
l
41
l
42
l
43
⋱
⋮
⋮
⋮
⋱
1
l
n
1
l
n
2
l
n
3
⋯
l
n
,
n
−
1
1
]
.
(2)
\boldsymbol{L} = (\boldsymbol{L^{(1)}})^{-1}(\boldsymbol{L^{(2)}})^{-1}\cdots(\boldsymbol{L^{(n-1)}})^{-1} = \left[ \begin{matrix} 1 \\ l_{21} & 1 \\ l_{31} & l_{32} & 1 \\ l_{41} & l_{42} & l_{43} & \ddots \\ \vdots & \vdots & \vdots & \ddots & 1 \\ l_{n1} & l_{n2} & l_{n3} & \cdots & l_{n, n-1} & 1 \end{matrix} \right]. \tag{2}
L=(L(1))−1(L(2))−1⋯(L(n−1))−1=⎣⎢⎢⎢⎢⎢⎢⎢⎡1l21l31l41⋮ln11l32l42⋮ln21l43⋮ln3⋱⋱⋯1ln,n−11⎦⎥⎥⎥⎥⎥⎥⎥⎤.(2)
该矩阵是对角线都是1的下三角矩阵,称为单位下三角矩阵,则可得
A
=
L
U
.
(3)
\boldsymbol{A} = \boldsymbol{L}\boldsymbol{U}. \tag{3}
A=LU.(3)
这样就把矩阵
A
\boldsymbol{A}
A分解为一个单位下三角矩阵
L
\boldsymbol{L}
L与一个上三角矩阵
U
\boldsymbol{U}
U的乘积。
一般有如下的定义:
如果方阵
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{LU}
LU分解肯定不是唯一的,因为上面提到的杜利特分解和克劳特分解本身就是两种特殊的三角分解方式。其实方阵的三角分解的方式有无穷多个,因为假设
D
\boldsymbol{D}
D是行列式不为零的任意对角矩阵,那么
A
=
L
U
=
L
D
D
−
1
U
=
(
L
D
)
(
D
−
1
U
)
=
L
~
U
~
\boldsymbol{A=L U=L D D^{-1} U=(L D)(D^{-1}U)=\widetilde{L} \widetilde{U}}
A=LU=LDD−1U=(LD)(D−1U)=L
U
其中
L
~
,
U
~
\widetilde{L}, \widetilde{U}
L
,U
也分别是下三角矩阵和上三角矩阵,然后
A
=
L
~
U
~
\boldsymbol{A}=\widetilde{L} \widetilde{U}
A=L
U
也是
A
\boldsymbol{A}
A 的一个三角分解。又因为
D
\boldsymbol{D}
D的任意性,所以三角分解不惟一。
但是可以判断满足什么条件一个方阵存在三角分解?什么形式的三角分解才是惟一的?
LDU基本定理
定理内容为
A
\boldsymbol{A}
A为
n
n
n阶方阵,则
A
\boldsymbol{A}
A可以惟一地分解为
A
=
L
D
U
(4)
\boldsymbol{A=LDU} \tag{4}
A=LDU(4)
的充分必要条件是
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1个顺序主子式
Δ
k
≠
0
(
k
=
1
,
2
,
⋯
,
n
−
1
)
\Delta_{k} \not = 0(k=1,2,\cdots,n-1)
Δk=0(k=1,2,⋯,n−1)(也就是在高斯消元过程文中说的高斯消元过程能完整进行下去的条件)。
其中
L
,
U
\boldsymbol{L, U}
L,U分别是单位下下三角矩阵和单位上三角矩阵,
D
\boldsymbol{D}
D是对角矩阵,形式如下
D
=
diag
(
d
1
,
d
2
,
⋯
,
d
n
)
d
k
=
Δ
k
Δ
k
−
1
=
a
k
k
(
k
)
,
k
=
1
,
2
,
⋯
,
n
D=\operatorname{diag}(d_{1}, d_{2}, \cdots, d_{n}) \\ d_{k}=\frac{\Delta_{k}}{\Delta_{k-1}}=a_{kk}^{(k)}, \quad k=1,2, \cdots, n
D=diag(d1,d2,⋯,dn)dk=Δk−1Δk=akk(k),k=1,2,⋯,n
下面会分别证明这个条件的充分性和必要性。
充分性证明
高斯消元过程文中定理2可知如果
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1个顺序主子式
Δ
k
≠
0
(
k
=
1
,
2
,
⋯
,
n
−
1
)
\Delta_{k} \not = 0(k=1,2,\cdots,n-1)
Δk=0(k=1,2,⋯,n−1),高斯消元过程得以完成,把矩阵
A
\boldsymbol{A}
A分解成一个单位下三角矩阵和上三角矩阵乘积,即实现了一个杜利特分解,
A
=
L
U
~
\boldsymbol{A=L\widetilde{U}}
A=LU
,其中
L
\boldsymbol{L}
L为单位下三角矩阵,
U
~
\boldsymbol{\widetilde{U}}
U
为上三角矩阵。
U
~
=
[
u
11
u
12
⋯
u
1
n
u
22
⋯
u
2
n
⋱
⋮
u
n
n
]
=
[
a
11
(
1
)
a
12
(
1
)
⋯
a
1
n
(
1
)
a
22
(
2
)
⋯
a
2
n
(
2
)
⋱
⋮
a
n
n
(
n
)
]
=
A
(
n
)
,
\boldsymbol{\widetilde{U}}=\left[ \begin{matrix} u_{11} & u_{12} & \cdots & u_{1 n} \\ & u_{22} & \cdots & u_{2 n} \\ & & \ddots & \vdots \\ & & & u_{n n} \end{matrix} \right] = \left[ \begin{matrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1 n}^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2 n}^{(2)} \\ & & \ddots & \vdots \\ & & & a_{n n}^{(n)} \end{matrix} \right] = \boldsymbol{A^{(n)}},
U
=⎣⎢⎢⎢⎡u11u12u22⋯⋯⋱u1nu2n⋮unn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡a11(1)a12(1)a22(2)⋯⋯⋱a1n(1)a2n(2)⋮ann(n)⎦⎥⎥⎥⎥⎤=A(n),
同时
u
i
i
≡
a
i
i
(
i
)
≠
0
(
i
=
1
,
2
,
⋯
,
n
−
1
)
u_{i i} \equiv a_{i i}^{(i)} \not = 0(i=1,2, \cdots, n-1)
uii≡aii(i)=0(i=1,2,⋯,n−1)。
存在性证明
分两种情况讨论:
- 当
A
\boldsymbol{A}
A非奇异,那么
Δ
n
=
a
11
(
1
)
a
22
(
2
)
⋯
a
n
n
(
n
)
=
∣
A
∣
≠
0
\Delta_{n}=a_{11}^{(1)} a_{22}^{(2)} \cdots a_{n n}^{(n)}=|A| \not = 0
Δn=a11(1)a22(2)⋯ann(n)=∣A∣=0,所以
a
n
n
(
n
)
=
u
n
n
≠
0
a_{n n}^{(n)} = u_{n n } \not = 0
ann(n)=unn=0。然后令
D
=
diag
(
a
11
(
1
)
,
a
22
(
2
)
,
⋯
,
a
n
n
(
n
)
)
\boldsymbol{D}=\operatorname{diag}(a_{11}^{(1)}, a_{22}^{(2)}, \cdots, a_{n n}^{(n)})
D=diag(a11(1),a22(2),⋯,ann(n)),则
D − 1 = diag ( 1 a 11 ( 1 ) , 1 a 22 ( 2 ) , ⋯ , 1 a n n ( n ) ) , D^{-1}=\operatorname{diag}(\frac{1}{a_{11}^{(1)}}, \frac{1}{a_{22}^{(2)}}, \cdots, \frac{1}{a_{n n}^{(n)}}), D−1=diag(a11(1)1,a22(2)1,⋯,ann(n)1),
所以
A = L U ~ = L D ( D − 1 U ~ ) = L D U (5) \boldsymbol{A=L \widetilde{U}=L D\left(D^{-1} \widetilde{U}\right)=L D U} \tag{5} A=LU =LD(D−1U )=LDU(5)
是 A \boldsymbol{A} A的一个 L D U \boldsymbol{L D U} LDU分解。 - 当
A
\boldsymbol{A}
A奇异, 那么
a
n
n
(
n
)
≡
u
n
n
=
0
a_{n n}^{(n)} \equiv u_{n n } = 0
ann(n)≡unn=0,令
D
=
diag
(
a
11
(
1
)
,
⋯
,
a
n
−
1
,
n
−
1
(
n
−
1
)
,
0
)
\boldsymbol{D}=\operatorname{diag}(a_{11}^{(1)}, \cdots, a_{n-1, n-1}^{(n-1)}, 0)
D=diag(a11(1),⋯,an−1,n−1(n−1),0),
D
n
−
1
=
diag
(
a
11
(
1
)
,
⋯
,
a
n
−
1
,
n
−
1
(
n
−
1
)
)
\boldsymbol{D_{n-1}}=\operatorname{diag}(a_{11}^{(1)}, \cdots, a_{n-1, n-1}^{(n-1)})
Dn−1=diag(a11(1),⋯,an−1,n−1(n−1)),
α
=
(
u
1
n
,
⋯
,
u
n
−
1
,
n
)
T
\boldsymbol{\alpha}=\left(u_{1 n}, \cdots, u_{n-1, n}\right)^\mathrm{T}
α=(u1n,⋯,un−1,n)T,则
U ~ ≡ [ U ~ n − 1 α 0 T 0 ] = [ D n − 1 0 0 T 0 ] [ D n − 1 − 1 U ~ n − 1 D n − 1 − 1 α 0 T 1 ] = D U \widetilde{\boldsymbol{U}} \equiv \left[\begin{matrix} \widetilde{\boldsymbol{U}}_{n-1} & \boldsymbol{\alpha} \\ \boldsymbol{0}^{\mathrm{T}} & 0 \end{matrix}\right]= \left[\begin{matrix} \boldsymbol{D}_{n-1} & \boldsymbol{0} \\ \boldsymbol{0}^{\mathrm{T}} & 0 \end{matrix}\right] \left[\begin{matrix} \boldsymbol{D}_{n-1}^{-1} \widetilde{\boldsymbol{U}}_{n-1} & \boldsymbol{D}_{n-1}^{-1} \boldsymbol{\alpha} \\ \boldsymbol{0}^{\mathrm{T}} & 1 \end{matrix}\right]= \boldsymbol{D} \boldsymbol{U} U ≡[U n−10Tα0]=[Dn−10T00][Dn−1−1U n−10TDn−1−1α1]=DU
所以不论在何种情况下,只要 Δ k ≠ 0 ( k = 1 , 2 , ⋯ , n − 1 ) \Delta_{k} \not = 0(k=1,2,\cdots,n-1) Δk=0(k=1,2,⋯,n−1),总是存在一个 L D U \boldsymbol{L D U} LDU分解,且 d k = a k k ( k ) = Δ k Δ k − 1 ( k = 1 , 2 , ⋯ , n ) , Δ 0 = 1 d_{k}=a_{kk}^{(k)}=\frac{\Delta_{k}}{\Delta_{k-1}}(k=1,2, \cdots, n),\Delta_0=1 dk=akk(k)=Δk−1Δk(k=1,2,⋯,n),Δ0=1。
唯一性证明
也分两种情况讨论:
- 当
A
\boldsymbol{A}
A非奇异时,
∣
A
∣
=
∣
L
∣
∣
D
∣
∣
U
∣
≠
0
\boldsymbol{|A|=|L||D||U|} \neq 0
∣A∣=∣L∣∣D∣∣U∣=0,所以
L
,
D
,
U
\boldsymbol{L, D, U}
L,D,U都是非奇异。如果还矩阵
A
\boldsymbol{A}
A还存在另一个
L
,
D
,
U
\boldsymbol{L, D, U}
L,D,U分解
A
=
L
1
D
1
U
1
\boldsymbol{A=L_{1} D_{1} U_{1}}
A=L1D1U1,这里
L
1
,
D
1
,
U
1
\boldsymbol{L_{1}, D_{1}, U_{1}}
L1,D1,U1也非奇异,于是有
L D U = L 1 D 1 U 1 (6) \boldsymbol{L D U}=\boldsymbol{L_{1} D_{1} U_{1}} \tag{6} LDU=L1D1U1(6)
上式两端左乘以 L 1 − 1 \boldsymbol{L}_{1}^{-1} L1−1以及右乘以 U − 1 \boldsymbol{U}^{-1} U−1和 D − 1 \boldsymbol{D}^{-1} D−1,得
L 1 − 1 L = D 1 U 1 U − 1 D − 1 \boldsymbol{L}_{1}^{-1} \boldsymbol{L}=\boldsymbol{D}_{1} \boldsymbol{U}_{1} \boldsymbol{U}^{-1} \boldsymbol{D}^{-1} L1−1L=D1U1U−1D−1
但上式左端是单位下三角阵,右端是单位上三角阵,所以都应该是单位阵。因此
L 1 − 1 L = I , D 1 U 1 U − 1 D − 1 = I \boldsymbol{L}_{1}^{-1} \boldsymbol{L}=\boldsymbol{I}, \quad \boldsymbol{D}_{1} \boldsymbol{U}_{1} \boldsymbol{U}^{-1} \boldsymbol{D}^{-1}=\boldsymbol{I} L1−1L=I,D1U1U−1D−1=I
即
L 1 = L , U 1 U − 1 = D 1 − 1 D \boldsymbol{L}_{1}=\boldsymbol{L}, \quad \boldsymbol{U}_{1} \boldsymbol{U}^{-1}=\boldsymbol{D}_{1}^{-1} \boldsymbol{D} L1=L,U1U−1=D1−1D
由后一个等式类似地可得
U 1 U 1 = I , D 1 − 1 D = I \boldsymbol{U}_{1} \boldsymbol{U}^{1}=\boldsymbol{I}, \quad \boldsymbol{D}_{1}^{-1} \boldsymbol{D}=\boldsymbol{I} U1U1=I,D1−1D=I
即有
U 1 = U , D 1 = D \boldsymbol{U}_{1}=\boldsymbol{U}, \quad \boldsymbol{D}_{1}=\boldsymbol{D} U1=U,D1=D - 当
A
A
A奇异,则公式
(
6
)
(6)
(6)可写成分块形式
[ L ~ 1 0 β 1 T 1 ] [ D ~ 1 0 0 T 0 ] [ U ~ 1 α 1 0 T 1 ] = [ L ~ 0 β T 1 ] [ D ~ 0 0 T 0 ] [ U ~ α 0 T 1 ] , \left[\begin{matrix} \widetilde{\boldsymbol{L}}_{1} & \boldsymbol{0} \\ \boldsymbol{\beta}_{1}^{\mathrm{T}} & 1 \end{matrix}\right] \left[\begin{matrix} \widetilde{\boldsymbol{D}}_{1} & \boldsymbol{0} \\ \boldsymbol{0}^{\mathrm{T}} & 0 \end{matrix}\right] \left[\begin{matrix} \widetilde{\boldsymbol{U}}_{1} & \boldsymbol{\alpha}_{1} \\ \boldsymbol{0}^{\mathrm{T}} & 1 \end{matrix}\right]= \left[\begin{matrix} \widetilde{\boldsymbol{L}} & \boldsymbol{0} \\ \boldsymbol{\beta}^{\mathrm{T}} & 1 \end{matrix}\right] \left[\begin{matrix} \widetilde{\boldsymbol{D}} & \boldsymbol{0} \\ \boldsymbol{0}^{\mathrm{T}} & 0 \end{matrix}\right] \left[\begin{matrix} \widetilde{\boldsymbol{U}} & \boldsymbol{\alpha} \\ \boldsymbol{0}^{\mathrm{T}} & 1 \end{matrix}\right], [L 1β1T01][D 10T00][U 10Tα11]=[L βT01][D 0T00][U 0Tα1],
其中 L ~ , L ~ 1 \widetilde{\boldsymbol{L}}, \widetilde{\boldsymbol{L}}_{1} L ,L 1是 n − 1 n-1 n−1阶单位下三角阵, U ~ , U ~ 1 \widetilde{\boldsymbol{U}}, \widetilde{\boldsymbol{U}}_{1} U ,U 1是 n − 1 n-1 n−1阶上三角阵, D ~ , D ~ 1 \widetilde{\boldsymbol{D}}, \widetilde{\boldsymbol{D}}_{1} D ,D 1是 n − 1 n-1 n−1阶对角阵, α , α 1 , β , β 1 \boldsymbol{\alpha}, \boldsymbol{\alpha}_{1}, \boldsymbol{\beta}, \boldsymbol{\beta}_{1} α,α1,β,β1是 n − 1 n-1 n−1维列向量。
由此得出
[ L ~ 1 D ~ 1 U ~ 1 L ~ 1 D ~ 1 α 1 β 1 T D ~ 1 U ~ 1 β 1 T D ~ 1 α 1 ] = [ L ~ D ~ U ~ L ~ D ~ α β T D ~ U ~ β T D ~ α ] , \left[\begin{matrix} \widetilde{\boldsymbol{L}}_{1}\widetilde{\boldsymbol{D}}_{1}\widetilde{\boldsymbol{U}}_{1} & \widetilde{\boldsymbol{L}}_{1}\widetilde{\boldsymbol{D}}_{1}\boldsymbol{\alpha}_{1} \\ \boldsymbol{\beta}_{1}^{\mathrm{T}}\widetilde{\boldsymbol{D}}_{1}\widetilde{\boldsymbol{U}}_{1} & \boldsymbol{\beta}_{1}^{\mathrm{T}}\widetilde{\boldsymbol{D}}_{1}\boldsymbol{\alpha}_{1} \end{matrix}\right] = \left[\begin{matrix} \widetilde{\boldsymbol{L}}\widetilde{\boldsymbol{D}}\widetilde{\boldsymbol{U}} & \widetilde{\boldsymbol{L}}\widetilde{\boldsymbol{D}}\boldsymbol{\alpha} \\ \boldsymbol{\beta}^{\mathrm{T}}\widetilde{\boldsymbol{D}}\widetilde{\boldsymbol{U}} & \boldsymbol{\beta}^{\mathrm{T}}\widetilde{\boldsymbol{D}}\boldsymbol{\alpha} \end{matrix}\right], [L 1D 1U 1β1TD 1U 1L 1D 1α1β1TD 1α1]=[L D U βTD U L D αβTD α],
其中 L ~ 1 , D ~ 1 , U ~ 1 \widetilde{\boldsymbol{L}}_{1}, \widetilde{\boldsymbol{D}}_{1}, \widetilde{\boldsymbol{U}}_{1} L 1,D 1,U 1和 L ~ , D ~ , U ‾ \widetilde{\boldsymbol{L}}, \widetilde{\boldsymbol{D}}, \overline{\boldsymbol{U}} L ,D ,U皆非奇昇,类似于前面的推理,可得
L ~ 1 = L ~ , D ~ 1 = D ~ , U ~ 1 = U ~ α 1 = α , β 1 T = β T \begin{matrix} \widetilde{\boldsymbol{L}}_{1}=\widetilde{\boldsymbol{L}}, & \widetilde{\boldsymbol{D}}_{1}=\widetilde{\boldsymbol{D}}, & \widetilde{\boldsymbol{U}}_{1}=\widetilde{\boldsymbol{U}} \\ \boldsymbol{\alpha}_{1}=\boldsymbol{\alpha}, & \boldsymbol{\beta}_{1}^{\mathrm{T}}=\boldsymbol{\beta}^{\mathrm{T}} \end{matrix} L 1=L ,α1=α,D 1=D ,β1T=βTU 1=U
充分性证明完毕!!!!
必要性证明
假定
A
\boldsymbol{A}
A有一个惟一的
L
D
U
\boldsymbol{L D U}
LDU分解,写成分块的形式便是
[
A
n
−
1
x
y
T
a
n
n
]
=
[
L
n
−
1
0
β
T
1
]
[
D
n
−
1
0
0
d
n
]
[
U
n
−
1
α
0
T
1
]
,
(7)
\left[\begin{matrix} \boldsymbol{A}_{n-1} & \boldsymbol{x} \\ \boldsymbol{y}^{\mathrm{T}} & a_{n n} \end{matrix}\right]= \left[\begin{matrix} \boldsymbol{L}_{n-1} & \boldsymbol{0} \\ \boldsymbol{\beta}^{\mathrm{T}} & 1 \end{matrix}\right] \left[\begin{matrix} \boldsymbol{D}_{n-1} & \boldsymbol{0} \\ \boldsymbol{0} & d_{n} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{U}_{n-1} & \boldsymbol{\alpha} \\ \boldsymbol{0}^{\mathrm{T}} & 1 \end{matrix}\right], \tag{7}
[An−1yTxann]=[Ln−1βT01][Dn−100dn][Un−10Tα1],(7)
其中
L
n
−
1
,
D
n
−
1
,
U
n
−
1
,
A
n
−
1
\boldsymbol{L_{n-1}, D_{n-1}, U_{n-1}, A_{n-1}}
Ln−1,Dn−1,Un−1,An−1分别是
L
,
D
,
U
,
A
\boldsymbol{L, D, U, A}
L,D,U,A的
n
−
1
n-1
n−1 阶顺序主子矩阵,
x
,
y
,
α
,
β
\boldsymbol{x, y, \alpha, \beta}
x,y,α,β为
n
−
1
n-1
n−1。
所以可知它们关系为
A
n
−
1
=
L
n
−
1
D
n
−
1
U
n
−
1
y
T
=
β
T
D
n
−
1
U
n
−
1
x
=
L
n
−
1
D
n
−
1
α
a
n
n
=
β
T
D
n
−
1
α
+
d
n
(8)
\begin{aligned} \boldsymbol{A}_{n-1}&=\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1} \boldsymbol{U}_{n-1} \\ \boldsymbol{y}^{\mathrm{T}}&=\boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{D}_{n-1} \boldsymbol{U}_{n-1} \\ \boldsymbol{x}&=\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1} \boldsymbol{\alpha} \\ a_{n n}&=\boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{D}_{n-1} \boldsymbol{\alpha}+d_{n} \end{aligned} \tag{8}
An−1yTxann=Ln−1Dn−1Un−1=βTDn−1Un−1=Ln−1Dn−1α=βTDn−1α+dn(8)
开始证明:
假设
Δ
n
−
1
=
∣
A
n
−
1
∣
=
0
\Delta_{n-1}=\left|\boldsymbol{A}_{n-1}\right|=0
Δn−1=∣An−1∣=0, 则
∣
A
n
−
1
∣
=
∣
L
n
−
1
∣
∣
D
n
−
1
∣
∣
U
n
−
1
∣
=
∣
D
n
−
1
∣
=
0
\left|\boldsymbol{A}_{n-1}\right|=\left|\boldsymbol{L}_{n-1}\right|\left|\boldsymbol{D}_{n-1}\right|\left|\boldsymbol{U}_{n-1}\right|=\left|\boldsymbol{D}_{n-1}\right|=0
∣An−1∣=∣Ln−1∣∣Dn−1∣∣Un−1∣=∣Dn−1∣=0
于是有
∣
L
n
−
1
D
n
−
1
∣
=
∣
D
n
−
1
∣
=
0
\left|\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1}\right|=\left|\boldsymbol{D}_{n-1}\right|=0
∣Ln−1Dn−1∣=∣Dn−1∣=0,即
L
n
−
1
D
n
−
1
\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1}
Ln−1Dn−1奇异,那么对于非齐次线性方程组
L
n
−
1
D
n
−
1
α
=
x
\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1} \boldsymbol{\alpha}=\boldsymbol{x}
Ln−1Dn−1α=x存在无穷多非零解。不妨设有
α
′
\boldsymbol{\alpha}^{\prime}
α′, 使
L
n
−
1
D
n
−
1
α
′
=
x
\boldsymbol{L}_{n-1} \boldsymbol{D}_{n-1} \boldsymbol{\alpha}^{\prime}=\boldsymbol{x}
Ln−1Dn−1α′=x,而
α
′
≠
α
\boldsymbol{\alpha}^{\prime} \neq \boldsymbol{\alpha}
α′=α。同理,因
D
n
−
1
U
n
−
1
\boldsymbol{D}_{n-1} \boldsymbol{U}_{n-1}
Dn−1Un−1奇异,故
(
D
n
−
1
U
n
−
1
)
T
=
D
n
−
1
T
U
n
−
1
T
\left(\boldsymbol{D}_{n-1} \boldsymbol{U}_{n-1}\right)^{\mathrm{T}}=\boldsymbol{D}_{n-1}^{\mathrm{T}} \boldsymbol{U}_{n-1}^{\mathrm{T}}
(Dn−1Un−1)T=Dn−1TUn−1T也奇异,就有
β
′
≠
β
\boldsymbol{\beta}^{\prime} \neq \boldsymbol{\beta}
β′=β,使
U
n
−
1
T
D
n
−
1
T
β
=
y
\boldsymbol{U}_{n-1}^{\mathrm{T}} \boldsymbol{D}_{n-1}^{\mathrm{T}} \boldsymbol{\beta}=\boldsymbol{y}
Un−1TDn−1Tβ=y和
U
n
−
1
T
D
n
−
1
T
β
′
=
y
\boldsymbol{U}_{n-1}^{\mathrm{T}} \boldsymbol{D}_{n-1}^{\mathrm{T}} \boldsymbol{\beta^{\prime}}=\boldsymbol{y}
Un−1TDn−1Tβ′=y。取
d
n
′
=
a
n
n
−
β
′
T
D
n
−
1
α
′
d_{n}^{\prime}=a_{n n}-\boldsymbol{\beta}^{\prime \mathrm{T}} \boldsymbol{D}_{n-1} \boldsymbol{\alpha}^{\prime}
dn′=ann−β′TDn−1α′,则有
[
A
n
−
1
x
y
T
a
n
n
]
=
[
L
n
−
1
0
β
′
T
1
]
[
D
n
−
1
0
0
d
n
′
]
[
U
n
−
1
α
′
0
T
1
]
.
\left[\begin{matrix} \boldsymbol{A}_{n-1} & \boldsymbol{x} \\ \boldsymbol{y}^{\mathrm{T}} & a_{n n} \end{matrix}\right]= \left[\begin{matrix} \boldsymbol{L}_{n-1} & \boldsymbol{0} \\ \boldsymbol{\beta}^{\prime \mathrm{T}} & 1 \end{matrix}\right] \left[\begin{matrix} \boldsymbol{D}_{n-1} & \boldsymbol{0} \\ \boldsymbol{0} & d_{n}^{\prime} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{U}_{n-1} & \boldsymbol{\alpha}^{\prime}\\ \boldsymbol{0}^{\mathrm{T}} & 1 \end{matrix}\right].
[An−1yTxann]=[Ln−1β′T01][Dn−100dn′][Un−10Tα′1].
这与
A
\boldsymbol{A}
A的
L
D
U
\boldsymbol{L D U}
LDU分解的惟一性矛盾,因此
Δ
n
−
1
≠
0
\Delta_{n-1} \neq 0
Δn−1=0。
同理可以把
n
−
1
n-1
n−1阶顺序主矩阵
A
n
−
1
\boldsymbol{A}_{n-1}
An−1也按照公式
(
7
)
(7)
(7)写成分块形式,同样有
A
n
−
2
=
L
n
−
2
D
n
−
2
U
n
−
2
\boldsymbol{A}_{n-2}=\boldsymbol{L}_{n-2} \boldsymbol{D}_{n-2} \boldsymbol{U}_{n-2}
An−2=Ln−2Dn−2Un−2,由于
∣
D
n
−
1
∣
≠
0
\left|\boldsymbol{D}_{n-1}\right| \neq 0
∣Dn−1∣=0,所以
∣
D
n
−
2
∣
≠
0
\left|\boldsymbol{D}_{n-2}\right| \neq 0
∣Dn−2∣=0,然后
∣
A
n
−
2
∣
=
∣
L
n
−
2
∣
∣
D
n
−
2
∣
∣
U
n
−
2
∣
=
∣
D
n
−
2
∣
≠
0
\left|\boldsymbol{A}_{n-2}\right|=\left|\boldsymbol{L}_{n-2}\right|\left|\boldsymbol{D}_{n-2}\right|\left|\boldsymbol{U}_{n-2}\right|=\left|\boldsymbol{D}_{n-2}\right| \neq 0
∣An−2∣=∣Ln−2∣∣Dn−2∣∣Un−2∣=∣Dn−2∣=0
从而
Δ
n
−
2
≠
0
\Delta_{n-2} \neq 0
Δn−2=0,以此类推可得
Δ
k
≠
0
(
k
=
1
,
2
,
⋯
,
n
−
1
)
\Delta_k \neq 0(k=1,2,\cdots,n-1)
Δk=0(k=1,2,⋯,n−1)。
必要性证明完毕!!!!
推论
在
A
\boldsymbol{A}
A的三角分解中,只要有一个三角矩阵是单位三角矩阵,则分解总是唯一的,上面说到的克劳特/杜利特分解是唯一的,而且可以进行杜利特/克劳特唯一分解的充分必要条件是
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1阶顺序主子式不为零。如果
A
\boldsymbol{A}
A是非奇异矩阵,那么充要条件则是
A
\boldsymbol{A}
A矩阵的各阶顺序主子式都不为零。如果
A
\boldsymbol{A}
A是奇异矩阵,杜利特分解的上三角矩阵
U
\boldsymbol{U}
U的最后一个对角线元素
u
n
n
=
0
u_{nn} = 0
unn=0,同理克劳特分解的下三角矩阵
L
\boldsymbol{L}
L的最后一个对角线元素
l
n
n
=
0
l_{nn} = 0
lnn=0。
否则,如果
L
\boldsymbol{L}
L和
U
\boldsymbol{U}
U都不是单位三角矩阵,那么分解不唯一。
稳定性
前面提到只要
n
n
n阶方阵
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1个顺序主子式
Δ
k
≠
0
(
k
=
1
,
2
,
⋯
,
n
−
1
)
\Delta_{k} \not = 0(k=1,2,\cdots,n-1)
Δk=0(k=1,2,⋯,n−1),就存在
L
U
\boldsymbol{L U}
LU分解,同时存在唯一的杜利特/克劳特分解以及
L
D
U
\boldsymbol{L D U}
LDU分解。
上面的条件主要是为了确保在高斯消元中作为被除数的
a
k
k
(
k
)
a_{kk}^{(k)}
akk(k)不能为零,但是如果
a
k
k
(
k
)
a_{kk}^{(k)}
akk(k)特别小的话,会导致最终分解矩阵的元素中出现特别大的数,导致分解的不稳定。
这里举个例子来看一下,比如矩阵
A
=
[
ϵ
1
1
1
]
\boldsymbol{A}=\left[\begin{matrix} \epsilon & 1 \\ 1 & 1 \end{matrix}\right]
A=[ϵ111]
L
D
U
\boldsymbol{LDU}
LDU分解为
A
=
[
1
1
ϵ
1
]
[
ϵ
1
−
1
ϵ
]
[
1
1
ϵ
1
]
=
L
D
U
\boldsymbol{A}= \left[\begin{matrix} 1 \\ \frac{1}{\epsilon} & 1 \end{matrix}\right] \left[\begin{matrix} \epsilon \\ & 1 - \frac{1}{\epsilon} \end{matrix}\right] \left[\begin{matrix} 1 & \frac{1}{\epsilon} \\ & 1 \end{matrix}\right] = \boldsymbol{LDU}
A=[1ϵ11][ϵ1−ϵ1][1ϵ11]=LDU
杜利特分解
A
=
[
1
1
ϵ
1
]
[
ϵ
1
1
−
1
ϵ
]
=
L
U
~
\boldsymbol{A}= \left[\begin{matrix} 1 \\ \frac{1}{\epsilon} & 1 \end{matrix}\right] \left[\begin{matrix} \epsilon & 1 \\ & 1 - \frac{1}{\epsilon} \end{matrix}\right] = \boldsymbol{L \widetilde{U}}
A=[1ϵ11][ϵ11−ϵ1]=LU
克劳特分解
A
=
[
ϵ
1
1
−
1
ϵ
]
[
1
1
ϵ
1
]
=
L
~
U
\boldsymbol{A}= \left[\begin{matrix} \epsilon \\ 1 & 1 - \frac{1}{\epsilon} \end{matrix}\right] \left[\begin{matrix} 1 & \frac{1}{\epsilon} \\ & 1 \end{matrix}\right] = \boldsymbol{\widetilde{L} U}
A=[ϵ11−ϵ1][1ϵ11]=L
U
上面三角分解中元素中出现
1
ϵ
\frac{1}{\epsilon}
ϵ1,是非常大的数,所以说
L
U
\boldsymbol{L U}
LU是不稳定的。
对于上面的例子来说,要解决大数不稳定问题比较简单,就是把矩阵
A
\boldsymbol{A}
A的第一行和第二行交换,然后在进行三角分解,基本的思路就是把小数下移大数上移,消除高斯消元过程中会出现的大数除以小数导致的不稳定问题。在这里就不展开讲了,上面的例子可以自己尝试行交换后的分解结果,发现不会出现
1
ϵ
\frac{1}{\epsilon}
ϵ1大数问题。
在后面的[矩阵的三角分解系列五] 三角分解中的行列变换里面会在详细的讲解这个方法。
例子
求矩阵
A
=
[
2
−
1
3
1
2
1
2
4
2
]
\boldsymbol{A}=\left[\begin{matrix} 2 & -1 & 3 \\ 1 & 2 & 1 \\ 2 & 4 & 2 \end{matrix}\right]
A=⎣⎡212−124312⎦⎤
的
L
D
U
\boldsymbol{LDU}
LDU分解。
解:
矩阵
A
\boldsymbol{A}
A的前
n
−
1
n-1
n−1阶顺序主子式为
Δ
1
=
2
Δ
2
=
5
\Delta_{1}=2 \\ \Delta_{2}=5
Δ1=2Δ2=5
所以矩阵
A
\boldsymbol{A}
A有惟一的
L
D
U
\boldsymbol{LDU}
LDU分解。
下面我们仿照高斯消元过程的计算步骤来计算矩阵
A
\boldsymbol{A}
A的
L
D
U
\boldsymbol{LDU}
LDU分解。
按照高斯消元过程中介绍的过程得到矩阵
A
\boldsymbol{A}
A消元阵
L
(
1
)
=
[
1
0
0
−
1
2
1
0
−
1
0
1
]
,
(
L
(
1
)
)
−
1
=
[
1
0
0
1
2
1
0
1
0
1
]
.
\boldsymbol{L}^{(1)}=\left[ \begin{matrix} 1 & 0 & 0 \\ -\frac{1}{2} & 1 & 0 \\ -1 & 0 & 1 \end{matrix} \right],\quad (\boldsymbol{L}^{(1)})^{-1}=\left[\begin{matrix} 1 & 0 & 0 \\ \frac{1}{2} & 1 & 0 \\ 1 & 0 & 1 \end{matrix}\right].
L(1)=⎣⎡1−21−1010001⎦⎤,(L(1))−1=⎣⎡1211010001⎦⎤.
所以得
A
(
2
)
=
L
(
1
)
A
=
[
2
−
1
3
0
5
2
−
1
2
0
5
−
1
]
.
\boldsymbol{A}^{(2)}= \boldsymbol{L}^{(1)} \boldsymbol{A}=\left[\begin{matrix} 2 & -1 & 3 \\ 0 & \frac{5}{2} & -\frac{1}{2} \\ 0 & 5 & -1 \end{matrix}\right].
A(2)=L(1)A=⎣⎡200−12553−21−1⎦⎤.
再由
A
(
2
)
\boldsymbol{A}^{(2)}
A(2)计算消元阵
L
(
2
)
=
[
1
0
0
0
1
0
0
−
2
1
]
,
(
L
(
2
)
)
−
1
=
[
1
0
0
0
1
0
0
2
1
]
.
\boldsymbol{L}^{(2)}=\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -2 & 1 \end{matrix}\right], \quad (\boldsymbol{L}^{(2)})^{-1}=\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 2 & 1 \end{matrix}\right].
L(2)=⎣⎡10001−2001⎦⎤,(L(2))−1=⎣⎡100012001⎦⎤.
得
A
(
3
)
=
L
(
2
)
A
(
2
)
=
[
2
−
1
3
0
5
2
−
1
2
0
0
0
]
=
U
~
=
D
U
\boldsymbol{A}^{(3)} = \boldsymbol{L}^{(2)} \boldsymbol{A}^{(2)}= \left[\begin{matrix} 2 & -1 & 3 \\ 0 & \frac{5}{2} & -\frac{1}{2} \\ 0 & 0 & 0 \end{matrix}\right]= \widetilde{\boldsymbol{U}} = \boldsymbol{D} \boldsymbol{U}
A(3)=L(2)A(2)=⎣⎡200−12503−210⎦⎤=U
=DU
又
D
=
[
a
11
(
1
)
0
0
0
a
22
(
2
)
0
0
0
a
33
(
3
)
]
=
[
2
0
0
0
5
2
0
0
0
0
]
,
(
D
)
−
1
=
[
1
2
0
0
0
2
5
0
0
0
0
]
.
\boldsymbol{D} = \left[\begin{matrix} a_{11}^{(1)} & 0 & 0 \\ 0 & a_{22}^{(2)} & 0 \\ 0 & 0 & a_{33}^{(3)} \end{matrix}\right] = \left[\begin{matrix} 2 & 0 & 0 \\ 0 & \frac{5}{2} & 0 \\ 0 & 0 & 0 \end{matrix}\right], \quad (\boldsymbol{D})^{-1}=\left[\begin{matrix} \frac{1}{2} & 0 & 0 \\ 0 & \frac{2}{5} & 0 \\ 0 & 0 & 0 \end{matrix}\right].
D=⎣⎢⎡a11(1)000a22(2)000a33(3)⎦⎥⎤=⎣⎡2000250000⎦⎤,(D)−1=⎣⎡21000520000⎦⎤.
所以
U
=
(
D
)
−
1
U
~
=
[
1
−
1
2
3
2
0
1
−
1
5
0
0
1
]
.
U = (\boldsymbol{D})^{-1} \widetilde{\boldsymbol{U}} = \left[\begin{matrix} 1 & -\frac{1}{2} & \frac{3}{2} \\ 0 & 1 & -\frac{1}{5} \\ 0 & 0 & 1 \end{matrix}\right].
U=(D)−1U
=⎣⎡100−211023−511⎦⎤.
同时
L
=
(
L
(
1
)
)
−
1
(
L
(
2
)
)
−
1
=
[
1
0
0
1
2
1
0
1
2
1
]
.
\boldsymbol{L} = (\boldsymbol{L}^{(1)})^{-1}(\boldsymbol{L}^{(2)})^{-1} = \left[\begin{matrix} 1 & 0 & 0 \\ \frac{1}{2} & 1 & 0 \\ 1 & 2 & 1 \end{matrix}\right].
L=(L(1))−1(L(2))−1=⎣⎡1211012001⎦⎤.
即
D
U
=
U
~
=
A
(
3
)
=
L
(
2
)
L
(
1
)
A
=
L
−
1
A
.
\boldsymbol{D} \boldsymbol{U} = \widetilde{\boldsymbol{U}} = \boldsymbol{A}^{(3)} = \boldsymbol{L}^{(2)}\boldsymbol{L}^{(1)}\boldsymbol{A} = \boldsymbol{L}^{-1}\boldsymbol{A}.
DU=U
=A(3)=L(2)L(1)A=L−1A.
故
A
=
[
1
0
0
1
2
1
0
1
2
1
]
[
2
0
0
0
5
2
0
0
0
0
]
[
1
−
1
2
3
2
0
1
−
1
5
0
0
1
]
=
L
D
U
\boldsymbol{A} = \left[\begin{matrix} 1 & 0 & 0 \\ \frac{1}{2} & 1 & 0 \\ 1 & 2 & 1 \end{matrix}\right] \left[\begin{matrix} 2 & 0 & 0 \\ 0 & \frac{5}{2} & 0 \\ 0 & 0 & 0 \end{matrix}\right] \left[\begin{matrix} 1 & -\frac{1}{2} & \frac{3}{2} \\ 0 & 1 & -\frac{1}{5} \\ 0 & 0 & 1 \end{matrix}\right] = \boldsymbol{L D U}
A=⎣⎡1211012001⎦⎤⎣⎡2000250000⎦⎤⎣⎡100−211023−511⎦⎤=LDU
引用
【1】 矩阵论(第二版)