[矩阵的三角分解系列二] LDU基本定理
简介矩阵的三角分解依赖的具体方法是之前讲到的高斯消元过程,所以看本文前建议先把高斯消元过程先看完。矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解,LDLT(不带平方根乔累斯基)分解等。介绍按照高斯消元过程文章中假设矩阵A\mathbf{A}A的前n−1n-1n−1个顺序主子式都不为零,那么矩
矩阵的三角分解是求解线性方程组常用的方法,包括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】 矩阵论(第二版)
更多推荐

所有评论(0)