目录
- Gauss 消去法的矩阵意义
- 矩阵的 LU 分解
- 其它的三角分解
-
- LDU 分解
- 对称正定矩阵
- 带状矩阵的分解
-
- 三对角阵的分解
- 矩阵分解的应用
-
- 1.求解方程组 $Ax=b$ 的解
- 2.求矩阵行列式
- 3.求矩阵的逆
- 4.求解成组方程组
Gauss 消去法的矩阵意义
注意：这里面讨论的不是列主元高斯消去法，所以算法可能不稳定。
根据《线性代数》我们知道，对于矩阵 A 行初等变换相当于左乘相应初等矩阵，对于矩阵 A 列初等变换相当于右乘相应初等矩阵。
由此可知：
A = L U A=LU A=LU
即 A 可以表示成为一个单位下三角阵 L 与一个上三角阵 U 的乘积
这就是矩阵的 LU 分解,也称为三角分解或 Doolittle 分解
矩阵的 LU 分解
那是不是所有的方阵都能够进行 LU 分解呢，如果不是，那又是哪些能够进行 LU 分解呢？
联系前面我们知道 Gauss 消元法顺利进行的条件是各个消元的主元不为 0 （即 a k k ( k − 1 ) ≠ 0 a_{kk}^{(k-1)} \neq 0 akk(k−1)̸=0 ）。我们知道，能够进行 Gauss 消元与能够进行 LU 分解是等价的， Gauss 消元法也可以用来求行列式。
若将矩阵 A 的左上方各阶顺序主子矩阵记为 A 1 , A 2 , . . . , A n A_1,A_2,...,A_n A1,A2,...,An 。使用高斯消元法使得第 k 个消元的主元如果为 0 ，那么 d e t ( A k ) = 0 det(A_k)=0 det(Ak)=0 ，所以有：
若 n 阶方阵的各阶顺序主子行列式不为 0 ，则存在唯一的单位下三角矩阵 L 和上三角矩阵 U ，满足 A=LU
由此定理可以得到 LU 分解的计算公式：
对于 A 中的每一项，有
因此得到 LU 分解的每项的计算公式：
{ μ 1 j = α 1 j , j = 1 , 2 , . . . , n l i 1 = α i 1 μ 11 , i = 2 , 3 , . . . , n μ i j = α i j − ∑ k = 1 i − 1 l i k u k j , j = i , i + 1 , . . . , n , i = 2 , 3 , . . , n l k i = 1 μ i i ( α k i − ∑ t = 1 i − 1 l k t μ t i ) , k = i + 1 , . . . , n , i = 2 , 3 , . . , n \left\{ \begin{array}{l} {\mu _{1j}} = {\alpha _{1j}}{\rm{ }},j = 1,2,...,n\\ {l_{i1}} = \frac{{{\alpha _{i1}}}}{{{\mu _{11}}}}{\rm{ }},i = 2,3,...,n\\ {\mu _{ij}} = {\alpha _{ij}} - \sum\limits_{k = 1}^{i - 1} {{l_{ik}}{u_{kj}}} {\rm{ }},j = i,i + 1,...,n,i = 2,3,..,n\\ {l_{ki}} = \frac{1}{{{\mu _{ii}}}}({\alpha _{ki}} - \sum\limits_{t = 1}^{i - 1} {{l_{kt}}{\mu _{ti}}} ){\rm{ }},k = i + 1,...,n,i = 2,3,..,n \end{array} \right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧μ1j=α1j,j=1,2,...,nli1=μ11αi1,i=2,3,...,nμij=αij−k=1∑i−1likukj,j=i,i+1,...,n,i=2,3,..,nlki=μii1(αki−t=1∑i−1lktμti),k=i+1,...,n,i=2,3,..,n
所以要先计算 U 中的第 i 行元素，再计算 L 中的第 i 列元素（因为计算 L 中的第 i 列元素需要知道 u i i u_{ii} uii ）。总的时间复杂度为 O ( n 3 ) O(n^3) O(n3) 。
我们在存储时，可以将 LU 存储在一个矩阵中，减小空间消耗：
( a 11 ( 0 ) a 12 ( 0 ) a 13 ( 0 ) ⋯ a 1 n ( 0 ) β 1 ( 0 ) l 21 a 22 ( 1 ) a 23 ( 1 ) ⋯ a 2 n ( 1 ) β 2 ( 1 ) l 31 l 32 a 33 ( 2 ) ⋯ a 3 n ( 2 ) β 3 ( 2 ) ⋱ ⋮ ⋮ l n 1 l n 2 l n 3 ⋯ a n n ( n − 1 ) β n ( n − 1 ) ) {\begin{pmatrix} {a_{11}^{(0)}}&{a_{12}^{(0)}}&{a_{13}^{(0)}}& \cdots &{a_{1n}^{(0)}}&{\beta _1^{(0)}}\\ {{l_{21}}}&{a_{22}^{(1)}}&{a_{23}^{(1)}}& \cdots &{a_{2n}^{(1)}}&{\beta _2^{(1)}}\\ {{l_{31}}}&{{l_{32}}}&{a_{33}^{(2)}}& \cdots &{a_{3n}^{(2)}}&{\beta _3^{(2)}}\\ {}&{}&{}& \ddots & \vdots & \vdots \\ {{l_{n1}}}&{{l_{n2}}}&{{l_{n3}}}& \cdots &{a_{nn}^{(n - 1)}}&{\beta _n^{(n - 1)}} \end{pmatrix}} ⎝⎜⎜⎜⎜⎜⎜⎛a11(0)l21l31ln1a12(0)a22(1)l32ln2a13(0)a23(1)a33(2)ln3⋯⋯⋯⋱⋯a1n(0)a2n(1)a3n(2)⋮ann(n−1)β1(0)β2(1)β3(2)⋮βn(n−1)⎠⎟⎟⎟⎟⎟⎟⎞
故而空间复杂度为 O ( n 2 ) O(n^2) O(n2) 。
其它的三角分解
LDU 分解
对于上三角矩阵 U ，如果 U 的对角元存在 0 ，则 d e t ( U ) = 0 det(U)=0 det(U)=0 。所以当 A 的三角分解存在时，上三角矩阵 U 的全体对角元 u k k ≠ 0 u_{kk}\neq 0 ukk̸=0,若令
D = d i a g ( u 11 , u 22 , . . . , u n n ) D=diag(u_{11},u_{22},...,u_{nn}) D=diag(u11,u22,...,unn)
则可以将 LU 分解转化为：
A = L U = L D D − 1 U = L D M T A=LU=LDD^{-1}U=LDM^T A=LU=LDD−1U=LDMT
我们知道 D − 1 = d i a g ( 1 u 11 , 1 u 22 , . . . , 1 u n n ) D^{-1}=diag(\frac 1 {u_{11}},\frac 1 {u_{22}},...,\frac 1 {u_{nn}}) D−1=diag(u111,u221,...,unn1), D − 1 U D^{-1}U D−1U 为单位上三角矩阵。
M T M^T MT 为单位上三角矩阵
M = ( 1 m 21 1 m 31 m 32 1 ⋮ ⋮ ⋮ ⋱ m n 1 m n 2 m n 3 ⋯ 1 ) , m j i = μ i j μ i i , ( i = 1 , 2 , . . . , n − 1 ; j = n + 1 , . . . , n ) M = {\begin{pmatrix} 1&{}&{}&{}&{}\\ {{m_{21}}}&1&{}&{}&{}\\ {{m_{31}}}&{{m_{32}}}&1&{}&{}\\ \vdots & \vdots & \vdots & \ddots &{}\\ {{m_{n1}}}&{{m_{n2}}}&{{m_{n3}}}& \cdots &1 \end{pmatrix}} ,{m_{ji}} = \frac{{{\mu _{ij}}}}{{{\mu _{ii}}}},\begin{array}{l} (i = 1,2,...,n - 1;\\ j = n + 1,...,n) \end{array} M=⎝⎜⎜⎜⎜⎜⎛1m21m31⋮mn11m32⋮mn21⋮mn3⋱⋯1⎠⎟⎟⎟⎟⎟⎞,mji=μiiμij,(i=1,2,...,n−1;j=n+1,...,n)
其中 LDU 分解和 LU 分解的条件是一样的。
我们知道，矩阵的 LU 分解能够通过计算式计算出来，那么 L 、 D 、 M 的元的公式也可以通过 LU 的元的公式计算出来。
{ d k = α k k − ∑ p = 1 k − 1 l k p ( d p m k p ) l i k = 1 d k [ α i k − ∑ p = 1 k − 1 l i p ( d p m k p ) ] m i k = 1 d k [ α k i − ∑ p = 1 k − 1 ( l k p d p ) m i p ) ] k = 1 , 2 , . . . , n , i = k + 1 , k + 2 , . . . , n \left\{ \begin{matrix} {d_k} = {\alpha _{kk}} - \sum\limits_{p = 1}^{k - 1} {{l_{kp}}({d_p}{m_{kp}})} {\rm{ }}\\ {l_{ik}} = \frac{1}{{{d_k}}}[{\alpha _{ik}} - \sum\limits_{p = 1}^{k - 1} {{l_{ip}}({d_p}{m_{kp}})]} \\ {m_{ik}} = \frac{1}{{{d_k}}}[{\alpha _{ki}} - \sum\limits_{p = 1}^{k - 1} {({l_{kp}}{d_p}){m_{ip}})]} {\rm{ }}\\ {\rm{ }}k = 1,2,...,n,i = k + 1,k + 2,...,n \end{matrix} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧dk=αkk−p=1∑k−1lkp(dpmkp)lik=dk1[αik−p=1∑k−1lip(dpmkp)]mik=dk1[αki−p=1∑k−1(lkpdp)mip)]k=1,2,...,n,i=k+1,k+2,...,n
时间复杂度为 O ( n 3 ) O(n^3) O(n3) 。
当 A 为对称矩阵时（即 A T = A A^T=A AT=A)，有 A = L D M T A=LDM^T A=LDMT ，则：
A = A T = M D L T A=A^T=MDL^T A=AT=MDLT,又因为 A 的 LDU 分解具有唯一性，所以：
M = L M=L M=L
所以，当 A 对称时，有分解式：
A = L D L T A=LDL^T A=LDLT
此时，只需要计算 L 和 D ，计算量减少约一半。
对称正定矩阵
任何对称正定矩阵 A 都存在唯一的下三角矩阵 G,满足 A = G G T A=GG^T A=GGT
理由如下：
在实际应用中，对称正定矩阵非常重要和常见。
在线性代数中，正定矩阵等价于顺序主子式全大于 0 ，即可以进行 LU 分解或者 LDU 分解。
只要 A 对称且可以分解，那么 A = L D L T A=LDL^T A=LDLT 成立，令 D 1 2 = d i a g ( d 1 , d 2 , . . . , d n , ) D^{\frac 1 2}=diag(\sqrt {d_1},\sqrt {d_2},...,\sqrt {d_n},) D21=diag(d1
,d2
,...,dn
,)
令 G = L D 1 2 G=LD^{\frac 1 2} G=LD21
则 A = L D L T = L D 1 2 D 1 2 L T = G G T A=LDL^T=LD^{\frac 1 2}D^{\frac 1 2}L^T=GG^T A=LDLT=LD21D21LT=GGT
称为 G 的Cholesky 分解
带状矩阵的分解
带状矩阵在有限元分析等实际应用场景非常普遍，所以非常重要。
定义：设 A = ( α i j ) A=(\alpha_{ij}) A=(αij) 是 n 阶方阵，对于小于 n 的正整数 p 和 q ，当 j>i+p 或 i>j+p 时，有 α i j = 0 \alpha_{ij}=0 αij=0 ，则称 A 是具有上带宽 q 和下带宽 p 的带状矩阵。
分解结果：
设 n 阶矩阵 A 有 LU 分解， A=LU,如果 A 是上下带宽分别为 q 和 p 的带状矩阵，则 L 是带宽为 p 的下三角矩阵， U 是带宽为 q 的上三角矩阵。
计算公式：
{ μ i j = α i j − ∑ t = max { 1 , i − p , j − q } i − 1 l i t μ t j , j = i , i + 1 , . . . , min { n , i + q } l k i = 1 μ i i ( α k i − ∑ t = max { 1 , k − p , i − q } i − 1 l i t μ t j ) , k = i + 1 , . . . , min { n , i + p } \left\{ \begin{array}{l} {\mu _{{\rm{ij}}}} = {\alpha _{ij}} - \sum\limits_{t = \max \{ 1,i - p,j - q\} }^{i - 1} {{l_{it}}{\mu _{tj}},j = i,i + 1,...,\min \{ n,i + q\} } \\ {l_{ki}} = \frac{1}{{{\mu _{ii}}}}({\alpha _{ki}} - \sum\limits_{t = \max \{ 1,k - p,i - q\} }^{i - 1} {{l_{it}}{\mu _{tj}}),k = i + 1,...,\min \{ n,i + p\} } \end{array} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧μij=αij−t=max{1,i−p,j−q}∑i−1litμtj,j=i,i+1,...,min{n,i+q}lki=μii1(αki−t=max{1,k−p,i−q}∑i−1litμtj),k=i+1,...,min{n,i+p}
只需要计算带状矩阵中不为 0 的元素，所以运算量比 O ( n 3 ) O(n^3) O(n3) 小得多， p 和 q 固定时，运算量是 O ( n ) O(n) O(n)
三对角阵的分解
当 p=q=1 时，此时的带状矩阵称为三对角阵
T = ( b 1 c 1 a 2 b 2 c 2 a 3 b 3 c 3 ⋱ ⋱ ⋱ a n − 1 b n − 1 c n − 1 a n b n ) {\rm{ }}T = {\begin{pmatrix} {{b_1}}&{{c_1}}&{}&{}&{}&{}\\ {{a_2}}&{{b_2}}&{{c_2}}&{}&{}&{}\\ {}&{{a_3}}&{{b_3}}&{{c_3}}&{}&{}\\ {}&{}& \ddots & \ddots & \ddots &{}\\ {}&{}&{}&{{a_{n - 1}}}&{{b_{n - 1}}}&{{c_{n - 1}}}\\ {}&{}&{}&{}&{{a_n}}&{{b_n}} \end{pmatrix}} {\rm{ }} T=⎝⎜⎜⎜⎜⎜⎜⎛b1a2c1b2a3c2b3⋱c3⋱an−1⋱bn−1ancn−1bn⎠⎟⎟⎟⎟⎟⎟⎞
根据带状矩阵分解的定理,上面带状矩阵可以被分解为：
L = ( 1 l 2 1 l 3 1 ⋱ ⋱ l n 1 ) U = ( μ 1 r 1 μ 2 r 2 ⋱ ⋱ μ n − 1 r n − 1 μ n ) L = \left( {\begin{matrix} 1&{}&{}&{}&{}\\ {{l_2}}&1&{}&{}&{}\\ {}&{{l_3}}&1&{}&{}\\ {}&{}& \ddots & \ddots &{}\\ {}&{}&{}&{{l_n}}&1 \end{matrix}} \right) U = \left( {\begin{matrix} {{\mu _1}}&{{r_1}}&{}&{}&{}\\ {}&{{\mu _2}}&{{r_2}}&{}&{}\\ {}&{}& \ddots & \ddots &{}\\ {}&{}&{}&{{\mu _{n - 1}}}&{{r_{n - 1}}}\\ {}&{}&{}&{}&{{\mu _n}} \end{matrix}} \right) L=⎝⎜⎜⎜⎜⎛1l21l31⋱⋱ln1⎠⎟⎟⎟⎟⎞U=⎝⎜⎜⎜⎜⎛μ1r1μ2r2⋱⋱μn−1rn−1μn⎠⎟⎟⎟⎟⎞
因此得到计算公式：
{ μ 1 = b 1 , r 1 = c 1 l k = α k μ k − 1 , μ k = b k − l k r k − 1 , r k = c k . k = 2 , 3 , . . . , n \left\{ \begin{matrix} {\mu _1} = {b_1},{r_1} = {c_1}\\ {l_k} = \frac{{{\alpha _k}}}{{{\mu _{k - 1}}}},{\mu _k} = {b_k} - {l_k}{r_{k - 1}},{r_k} = {c_k}.k = 2,3,...,n \end{matrix} \right. {μ1=b1,r1=c1lk=μk−1αk,μk=bk−lkrk−1,rk=ck.k=2,3,...,n
求解三对角方程组： T x = d Tx=d Tx=d,其中 d = ( d 1 , d 2 , . . . , d n ) T d={(d_1,d_2,...,d_n)}^T d=(d1,d2,...,dn)T ，将 T 进行三角分解 T=LU ，则 LUx=d,令 y=Ux ，则 Ly=d,得到计算 x 和 y 的公式：
可以推导出：
追的过程（从前往后）：
{ y 1 = d 1 y k = d k − l k y k − 1 , k = 2 , 3 , . . . , n \left\{ \begin{array}{l} {y_1} = {d_1}\\ {y_k} = {d_k} - {l_k}{y_{k - 1}},k = 2,3,...,n \end{array} \right. {y1=d1yk=dk−lkyk−1,k=2,3,...,n
赶的过程（从后往前）：
{ x n = y n μ n x k = 1 μ k ( y k − r k x k + 1 ) , k = n − 1 , n − 2 , . . . , 2 , 1 \left\{ \begin{array}{l} {x_n} = \frac {y_n} {\mu _n}\\ {x_k} = \frac{1}{{{\mu _k}}}({y_k} - {r_k}{x_{k + 1}}),k = n - 1,n - 2,...,2,1 \end{array} \right. {xn=μnynxk=μk1(yk−rkxk+1),k=n−1,n−2,...,2,1
所以该方法被称为追赶法。
矩阵分解的应用
若对 A 进行 LU 分解
1.求解方程组 A x = b Ax=b Ax=b 的解
原方程等价于：
{ L y = b U x = y \left\{ \begin{array}{l} Ly = b\\ Ux = y \end{array} \right. {Ly=bUx=y
则可以通过类似三对角方程组追赶法求解。 LU 分解的时间复杂度为 O ( n 3 ) O(n^3) O(n3) ，之后的求解过程的时间复杂度为 O ( n 2 ) O(n^2) O(n2)
2.求矩阵行列式
det ( A ) = det ( L ) det ( U ) \det (A) = \det (L)\det (U) det(A)=det(L)det(U)
因为 L 是单位下三角矩阵，所以
det ( A ) = det ( L ) det ( u ) = det ( u ) = μ 11 μ 22 ⋯ μ n n \det (A) = \det (L)\det (u) = \det (u) = {\mu _{11}}{\mu _{22}} \cdots {\mu _{nn}} det(A)=det(L)det(u)=det(u)=μ11μ22⋯μnn
3.求矩阵的逆
对称矩阵 A 进行 LDU 分解有 A = L D L T A=LDL^T A=LDLT ， A − 1 = L − T D − 1 L − 1 {A^{ - 1}} = {L^{ - T}}{D^{ - 1}}{L^{ - 1}} A−1=L−TD−1L−1
D − 1 = d i a g ( d 1 − 1 , d 2 − 1 , . . . , d n − 1 ) {D^{ - 1}} = diag(d_1^{ - 1},d_2^{ - 1},...,d_n^{ - 1}) D−1=diag(d1−1,d2−1,...,dn−1)
L − 1 L^{-1} L−1 是单位下三角矩阵，由追赶法可以解出 L − 1 L^{-1} L−1 ，进而求出 A − 1 A^{-1} A−1
4.求解成组方程组
若有一系列方程组
A x = B k ( k = 1 , . . . , m , m ≥ 1 ) Ax = {B_k}{\rm{ (}}k = 1,...,m,m \ge 1) Ax=Bk(k=1,...,m,m≥1)
若 A 存在 LU 分解， A=LU ，
{ L y = B k U x = y \left\{ \begin{array}{l} Ly = {B_k}\\ Ux = y \end{array} \right. {Ly=BkUx=y
对每个 B k {B_k} Bk 分别按照求解方程组的方法就可以求解成组方程组。