线性方程组的直接法

三角分解

张瑞
中国科学技术大学数学科学学院

三角分解

由矩阵理论，增广矩阵的一次行操作，与左乘一个初等矩阵是等价的。 在 Gauss 消元的第k步，相当与左乘了如下的初等矩阵

\[\begin{pmatrix} 1 & & & & \\ & \ddots & & & \\ & & 1 & \\ & & l_{k+1,k} & 1 & \\ & & & & \ddots & \\ & & & & & 1 \\ \end{pmatrix} \cdots \begin{pmatrix} 1 & & & & \\ & \ddots & & & \\ & & 1 & \\ & & & 1 & \\ & & & & \ddots & \\ & & l_{n,k} & & & 1 \\ \end{pmatrix} \]

其中$l_{i,k}=\frac{-a_{ik}^{(k)}}{a_{kk}^{(k)}}$。是一系列的单位下三角阵。

因此，整个Gauss消元过程，相当与左乘了一系列的单位三角阵，最后将矩阵变为了上三角阵。即有$LA=U$。利用单位下三角阵的逆矩阵仍然是单位三角阵， 可以证明如下结论

定理 1.
若矩阵$A$的所有顺序主子式都不为$0$，则存在单位下三角阵$L$，和上三角阵$U$，使得$A=LU$。

若有$A=LU$，则线性方程组$Ax=b$，可以按如下方式求得

\[\begin{cases} Ly=b \\ Ux=y \end{cases} \]

Dollitte 分解

将矩阵$A$分解为单位下三角阵$L$与上三角阵$U$的乘积的过程，称为Dollitte分解。即

\[A=\begin{pmatrix} 1 & \\ l_{21} & 1 & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots & \\ & & & u_{nn} \end{pmatrix} \]
  • 将$L$与$U$相乘后，与$A$的元素作比较，即可以顺序解得$L$与$U$的每个元素。
\[\begin{pmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} & \cdots & l_{21}u_{1n}+u_{2n} \\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33} & \cdots \\ \vdots & \vdots \\ l_{n1}u_{11} & l_{n1}u_{12}+l_{n2}u_{22} & l_{n1}u_{13}+l_{n2}u_{23}+l_{n3}u_{33} & \end{pmatrix} \]

比较第1行与第1列，有

\[\begin{aligned} u_{1j}=a_{1j}, & j=1,2,\cdots,n \\ l_{i1}=\frac{a_{i1}}{u_{11}}, & i=2,\cdots,n \end{aligned} \]

计算量为： $n-1$

比较第2行与第2列，有

\[\begin{aligned} u_{2j}=a_{2j}-l_{21}u_{1j}, & j=2,\cdots,n \\ l_{i2}=\frac{a_{i2}-l_{i1}u_{12}}{u_{22}}, & i=3,\cdots,n \end{aligned} \]

计算量为： $n-1+(n-2)2$

比较第k行与第k列，有

\[\begin{aligned} u_{kj}=a_{kj}-\sum_{m=1}^{k-1}l_{km}u_{mj}, & j=k,\cdots,n \\ l_{ik}=\frac{a_{ik}-\displaystyle\sum_{m=1}^{k-1}l_{im}u_{mk}}{u_{kk}}, & i=k+1,\cdots,n \end{aligned} \]

计算量为： $(n-k+1)(k-1)+(n-k)k$

整个LU分解的计算量为

\[\sum_{k=1}^n [(n-k+1)(k-1)+(n-k)k] \]

分解过程中，不需要新的存贮，可以直接覆盖原来的矩阵

\[\begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ l_{21} & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots & \\ l_{n1} & l_{n2} & \cdots & u_{nn} \end{pmatrix} \]
  • 分解能够进行的条件与Gauss消元是一样的，需要矩阵的所有顺序主子式均不为0
  • 分解只与矩阵有关，而Gauss消元法则是对增广矩阵的操作
  • 分解完成后，解线性方程组就只需要$n^2$的计算量

例 1. 将下面矩阵分解为单位下三角阵与上三角阵的乘积

\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ ? & 1 & 0 \\ ? & ? & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & ? & ? \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ ? & ? & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & ? & ? \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ ? & ? & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & 0.5 & 3 \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ 2 & ? & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & 0.5 & 3 \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ 2 & 3 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & 0.5 & 3 \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ 2 & 3 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -3 \\ 0 & 0.5 & 3 \\ 0 & 0 & 2.5 \end{pmatrix} \]

列主元的LU分解

  • 分解过程中，$u_{kk}$是除数。为了算法能够进行，需要保证$u_{kk}\neq 0$。
  • 为了减少传入误差的影响，$|u_{kk}|$应该尽可能的大。
  • 类似与列主元的Gauss消元，三角分解法也可以做列主元的操作
  • 列主元分解法会使用行的交换，在解方程的时候，需要把右端项先做同样的行交换，然后才能解出正确的解

分解第k行与第k列前，矩阵内的数据为

\[\begin{pmatrix} u_{11} & \cdots & u_{1,k-1} & u_{1,k} & \cdots & u_{1n} \\ \vdots & \ddots & \vdots & \vdots & & \vdots \\ l_{k-1,1} & \cdots & u_{k-1,k-1} & u_{k-1,k} & \cdots & u_{k-1,n} \\ l_{k,1} & \cdots & l_{k,k-1} & a_{k,k} & \cdots & a_{k,n} \\ \vdots & & \vdots & \vdots & \ddots & \vdots & \\ l_{n1} & \cdots & l_{n,k-1} & a_{nk} & \cdots & a_{nn} \end{pmatrix} \]

取

\[\tilde u_{i,k}=a_{ik}-\sum_{m=1}^{k-1}l_{j,m}u_{m,k} , i=k,\cdots,n \]

若$\tilde u_{j,k}$是所有$\tilde u_{i,k}$中绝对值最大的，则交换矩阵的$j$行与$k$行。

例 2. 将下面矩阵分解为单位下三角阵与上三角阵的乘积

\[\begin{pmatrix} 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ 4 & 3.5 & 5.5 \end{pmatrix} \]

解. 查看第1列，交换第3行与第1行，得

\[\begin{pmatrix} 4 & 3.5 & 5.5 \\ 3 & 2 & -1.5 \\ 2 & 1 & -3 \\ \end{pmatrix} \]
\[\begin{pmatrix} 4 & 3.5 & 5.5 \\ 3 & 2 & -1.5 \\ 2 & 1 & -3 \\ \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0.75 & 1 & 0 \\ 0.5 & ? & 1 \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & ? & ? \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} 4 & 3.5 & 5.5 \\ 3 & 2 & -1.5 \\ 2 & 1 & -3 \\ \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0.75 & 1 & 0 \\ 0.5 & ? & 1 \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & ? & ? \\ 0 & 0 & ? \end{pmatrix} \]

矩阵中存储为

\[\begin{pmatrix} 4 & 3.5 & 5.5 \\ 0.75 & 2 & -1.5 \\ 0.5 & 1 & -3 \\ \end{pmatrix} \]

注意到

\[\tilde u_{22}=2-0.75*3.5=-0.625 , \tilde u_{32}=1-0.5*3.5=-0.75 \]

$|\tilde u_{22}|<|\tilde u_{32}|$，交换第2行与第3行

\[\begin{pmatrix} 4 & 3.5 & 5.5 \\ 0.5 & 1 & -3 \\ 0.75 & 2 & -1.5 \\ \end{pmatrix} \]

下面计算

\[\begin{pmatrix} * & * & * \\ * & 1 & -3 \\ * & 2 & -1.5 \\ \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.75 & ? & 1 \\ \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & ? & ? \\ 0 & 0 & ? \end{pmatrix} \]

得

\[\begin{pmatrix} * & * & * \\ * & 1 & -3 \\ * & 2 & -1.5 \\ \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.75 & 0.833333 & 1 \\ \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & -0.75 & -5.75 \\ 0 & 0 & ? \end{pmatrix} \]
\[\begin{pmatrix} * & * & * \\ * & 1 & -3 \\ * & 2 & -1.5 \\ \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.75 & 0.833333 & 1 \\ \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & -0.75 & -5.75 \\ 0 & 0 & -0.83333 \end{pmatrix} \]

将得到的下三角阵与上三角阵相乘后，可以得到

\[\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.75 & 0.833333 & 1 \\ \end{pmatrix} \begin{pmatrix} 4 & 3.5 & 5.5 \\ 0 & -0.75 & -5.75 \\ 0 & 0 & -0.83333 \end{pmatrix} =\begin{pmatrix} 4 & 3.5 & 5.5 \\ 2 & 1 & -3 \\ 3 & 2 & -1.5 \\ \end{pmatrix} \]

与交换过程是相匹配的。

Courant 分解

将矩阵$A$分解为下三角阵$L$与单位上三角阵$U$的乘积的过程，称为Courant 分解。即

\[A=\begin{pmatrix} l_{11} & \\ l_{21} & l_{22} & \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{pmatrix} \begin{pmatrix} 1 & u_{12} & \cdots & u_{1n} \\ & 1 & \cdots & u_{2n} \\ & & \ddots & \vdots & \\ & & & 1 \end{pmatrix} \]
  • 将$L$与$U$相乘后，与$A$的元素作比较，即可以顺序解得$L$与$U$的每个元素。
  • 与Dollitte分解类似，可以先比较第1列和第1行，依次到最后一个元素。

三对角阵的追赶法

对于三对角阵，

\[A=\begin{pmatrix} a_1 & b_1 \\ c_2 & a_2 & b_2 \\ & c_3 & \ddots & \ddots \\ & & \ddots & \ddots & b_{n-1} \\ & & & c_{n} & a_n \end{pmatrix} \]

可以分解为

\[A=\begin{pmatrix} \alpha_1 & \\ \gamma_2 & \alpha_2 & \\ & \gamma_3 & \ddots & \\ & & \ddots & \ddots & \\ & & & \gamma_{n} & \alpha_n \end{pmatrix} \begin{pmatrix} 1 & \beta_1 \\ & 1 & \beta_2 \\ & & \ddots & \ddots \\ & & & \ddots & \beta_{n-1} \\ & & & & 1 \end{pmatrix} \]

即有

\[\begin{cases} \gamma_i=c_i, i=2,3,\cdots,n \\ \gamma_i\beta_{i-1}+\alpha_i=a_i, i=1,2,\cdots,n \\ \alpha_i\beta_i=b_i, i=1,2,\cdots,n-1 \end{cases} \]

分解过程为

\[\begin{cases} \gamma_i=c_i, i=2,3,\cdots,n \\ \alpha_i=a_i-c_i\beta_{i-1} , i=1,2,\cdots,n \\ \beta_i=b_i/\alpha_i, i=1,2,\cdots,n-1 \end{cases} \]

解方程$Ax=f$的2步反代过程为

\[\begin{aligned} y_{i}=(f_i-c_iy_{i-1})/\alpha_i , i=1,2,\cdots,n \\ x_{i}=y_{i}-\beta_{i}x_{i+1} , i=n,n-1,\cdots,1 \end{aligned} \]

注意到，反代过程中，求$y_i$的过程，可以放到分解过程中完成。

整合后的算法为

\[\begin{cases} \alpha_i=a_i-c_i\beta_{i-1} \\ \beta_i=b_i/\alpha_i \\ y_{i}=(f_i-c_iy_{i-1})/\alpha_i \\ \end{cases}, i=1,2,\cdots,n \]
\[x_{i}=y_{i}-\beta_{i}x_{i+1} , i=n,n-1,\cdots,1 \]

称为追赶法

对称正定阵的$LDL^{T}$分解

当矩阵$A$对称正定时，则有$A=LL^{T}$，其中$L$是下三角阵。

即有

\[\begin{pmatrix} l_{11} & 0 & & 0 \\ l_{21} & l_{22} & & 0 \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & l_{nn}\\ \end{pmatrix} \begin{pmatrix} l_{11} & l_{21} & \cdots & l_{n1} \\ 0 & l_{22} & \cdots & l_{n2} \\ & & \ddots & \vdots \\ 0 & 0 & & l_{nn} \end{pmatrix} \]

类似与LU分解的过程，可以按一定顺序得到$L$的所有元素。

这样，得到了分解过程为

\[\begin{aligned} l_{j,j} & = \sqrt{ a_{j,j} - \sum_{k=1}^{j-1} l_{j,k}^2 } \\ l_{i,j} & = \frac{1}{l_{j,j}} \left( a_{i,j} - \sum_{k=1}^{j-1} l_{i,k} l_{j,k} \right) \quad \forall i>j. \end{aligned} \]

称为Cholesky分解，或平方根方法

例 3. 分解如下矩阵

\[\left( \begin{array}{*{3}{r}} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \\ \end{array} \right) = \left( \begin{array}{*{3}{r}} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \\ \end{array} \right) \left( \begin{array}{*{3}{r}} 2 & 6 & -8 \\ 0 & 1 & 5 \\ 0 & 0 & 3 \\ \end{array} \right) \]

将下三角阵$L$写成$L=\tilde L\tilde D$，其中$\tilde L$是单位下三角阵，$\tilde D$是对角阵。 则

\[A=\tilde L\tilde D\tilde D^T\tilde L^T \]

注意到，$\tilde D\tilde D^T$仍然是一个对角阵，则有

定理 2.
对称正定阵$A$可以写成$A=LDL^T$，其中$L$是单位下三角阵，$D$是对角阵。

\[A=\begin{pmatrix} 1 & 0 & & 0 \\ l_{21} & 1 & & 0 \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1\\ \end{pmatrix} \begin{pmatrix} d_1 & \\ & d_2 & \\ & & \ddots & \\ & & & d_n \end{pmatrix} \begin{pmatrix} 1 & l_{21} & \cdots & l_{n1} \\ 0 & 1 & \cdots & l_{n2} \\ & & \ddots & \vdots \\ 0 & 0 & & 1 \end{pmatrix} \]

类似与LU分解过程，$LDL^T$的分解过程为

\[\begin{aligned} d_j &= a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2 d_k \\ l_{ij} & = \frac{1}{d_j} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk} d_k \right) \quad \forall i>j \end{aligned} \]
  • 当$A$对称正定时，可以保证$d_j>0, \forall j$。
  • 分解可行的条件是$d_j\neq 0, \forall j$。这样，对一些非正定的对称阵，分解方法也是可行的。

例 4. 分解如下矩阵

\[\left( \begin{array}{*{3}{r}} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \\ \end{array} \right) \]

解. 可以分解为

\[\left( \begin{array}{*{3}{r}} 1 & 0 & 0 \\ 3 & 1 & 0 \\ -4 & 5 & 1 \\ \end{array} \right) \left( \begin{array}{*{3}{r}} 4 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 9 \\ \end{array} \right) \left( \begin{array}{*{3}{r}} 1 & 3 & -4 \\ 0 & 1 & 5 \\ 0 & 0 & 1 \\ \end{array} \right) \]

谢谢

vertical slide 2

目录

本节读完

例 5.

[#ex9-1-0].