数值线性代数1

今天晚上突然被人提到一下，然后就想起来了。原来我还可以翻译姐姐的博客来着。这一篇来自这里，是介绍的LU分解。

同样也是捣鼓MathJax也花费了大量的时间呢。

线性系统

线性方程的一般形式可以写成：

$${\displaystyle {\begin{aligned}a_{11}x_{1}&+a_{12}x_{2}+\cdots +a_{1n}x_{n}=b_{1}\\ a_{21}x_{1}&+a_{22}x_{2}+\cdots +a_{2n}x_{n}=b_{2}\\ \vdots &\\ a_{m1}x_{1}&+a_{m2}x_{2}+\cdots +a_{mn}x_{n}=b_{m},\end{aligned}}}.$$

这等价于矩阵方程 $ Ax = b $， 其中：

$$A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\ a_{21}&a_{22}&\cdots &a_{2n}\\ \vdots &\vdots &\ddots &\vdots \\ a_{m1}&a_{m2}&\cdots &a_{mn}\end{bmatrix}},\quad {\mathbf {x}}={\begin{bmatrix}x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end{bmatrix}},\quad {\mathbf {b}}={\begin{bmatrix}b_{1}\\ b_{2}\\ \vdots \\ b_{m}\end{bmatrix}}$$

LU分解

我们常用把一个矩阵分解为简单形式来解矩阵。例如：上三角矩阵，对角矩阵和酉矩阵。我将要在这里介绍臭名昭著的LU分解。LU分解是高斯消元的另一种说法，它通过在矩阵左侧应用线性变换（初等行变换）把一个完整的矩阵转换为一个上三角矩阵。

假设$A \in \Bbb{C}^{m,m}$是一个复数方阵¹，这表示这个系统有相同数量的方程和未知数。LU分解的想法是通过引入对角线元素以下的0将$A$转换为$m\times m$的上三角矩阵$U$。这可以通过减去每一行和他们相应的倍数做到。这个消元过程就是高斯消元，它等价于在左边给$A$乘上一系列的下三角矩阵 $L_{k}$ 。

$$L_{m-1}\cdots L_{2} L_{1} A = L^{-1} A = U.$$

另$L = L_{1}^{-1}L_{2}^{-1} \cdots L_{m-1}^{-1}$，我们有：

$$A=LU,$$

其中L是下三角矩阵，U是上三角矩阵。我们常常把L做成一个单位下三角矩阵，这样就是一次单独的LU分解。

例子：

代码可以在这里找到。

译注：这里的代码要求在Julia项目的master下才能运行，如果在Julia0.6版本，可以把其中这一行的uninitialized,删除，则可以正常运行。

-------------------         A        -------------------
  3  -7  -2    2
 -3   5   1    0
  6  -4   0   -5
 -9   5  -5  -12
-------------------       L_{1}     -------------------
  1.0   ⋅    ⋅    ⋅
  1.0  1.0   ⋅    ⋅
 -2.0  0.0  1.0   ⋅
  3.0  0.0  0.0  1.0
------------------- L_{1}A = U_{1} -------------------
 3.0   -7.0   -2.0   2.0
 0.0   -2.0   -1.0   2.0
 0.0   10.0    4.0  -9.0
 0.0  -16.0  -11.0  -6.0
-------------------       L_{2}     -------------------
 1.0    ⋅    ⋅    ⋅
 0.0   1.0   ⋅    ⋅
 0.0   5.0  1.0   ⋅
 0.0  -8.0  0.0  1.0
------------------ L_{2}U_{1} = U_{2} ------------------
 3.0  -7.0  -2.0    2.0
 0.0  -2.0  -1.0    2.0
 0.0   0.0  -1.0    1.0
 0.0   0.0  -3.0  -22.0
-------------------       L_{3}     -------------------
 1.0   ⋅     ⋅    ⋅
 0.0  1.0    ⋅    ⋅
 0.0  0.0   1.0   ⋅
 0.0  0.0  -3.0  1.0
-------------------         U       -------------------
 3.0  -7.0  -2.0    2.0
  ⋅   -2.0  -1.0    2.0
  ⋅     ⋅   -1.0    1.0
  ⋅     ⋅     ⋅   -25.0
-------------------      L^{-1}    -------------------
  1.0    ⋅     ⋅    ⋅
  1.0   1.0    ⋅    ⋅
 -2.0   5.0   1.0   ⋅
  3.0  -8.0  -3.0  1.0

现在为了取得完整的LU分解$A = LU$，我们需要计算$L = L_{1}^{-1} L_{2}^{-1}\cdots L_{m-1}^{-1}$。$L_{i}$的逆矩阵就是对角线元素下所有元素取负的它自己。而且$\prod_{k=1}^{k=m-1}L_{k}^{-1}$是具有非零子对角线元素的单位下三角矩阵。

因此我们有：

$$\begin{bmatrix} 1 & & & \\ 1 & 1 & & \\ -2 & 0 & 1 & \\ 3 & 0 & 0 & 1 \end{bmatrix}^{-1} = \begin{bmatrix} 1 & & & \\ -1 & 1 & & \\ 2 & 0 & 1 & \\ -3 & 0 & 0 & 1 \end{bmatrix}$$

最终我们得到了：

$$\begin{bmatrix} 3 & -7 & -2 & 2 \\ -3 & 5 & 1 & 0 \\ 6 & -4 & 0 & -5 \\ -9 & 5 & -5 & -12 \end{bmatrix} = \begin{bmatrix} 1 & & & \\ -1 & 1 & & \\ 2 & -5 & 1 & \\ -3 & 8 & 3 & 1 \end{bmatrix} \begin{bmatrix} 3 & -7 & -2 & 2 \\ & -2 & -1 & 2 \\ & & -1 & 1 \\ & & & -25 \end{bmatrix}$$

证明

$L_{i}$的逆矩阵就是对角线元素下所有元素取负的它自己。而且$\prod_{k=1}^{k=m-1}L_{k}^{-1}$是具有非零子对角线元素的单位下三角矩阵。

证明1:

我们可以将$L_{k}$重写为

$$\begin{align} L_{k} &= I - l_{k}e_{k}^{\dagger} \quad \text{其中 } l_{k} = \begin{bmatrix} 0\\ \vdots\\ 0\\ l_{k+1,k}\\ \vdots\\ l_{m,k} \end{bmatrix} \end{align}$$

很明显$e_{k}^{\dagger}l_{k} = 0$，因此：

$$(I-e_{k}^{\dagger}l_{k})(I+e_{k}^{\dagger}l_{k}) = I - l_{k}e_{k}^{\dagger}l_{k}e_{k}^{\dagger} = I.$$

也就是说$L_{k}^{-1} = I+l_{k}e_{k}^{\dagger}$. $\blacksquare$

证明2

从$l_{k+1}$的稀疏排列中，我们可以知道$e_{k}^{\dagger}l_{k+1} = 0$。因此

$$L_{k}^{-1}L_{k+1}^{-1} = (I+e_{k}^{\dagger}l_{k})(I+e_{k+1}^{\dagger}l_{k+1}) = I + l_{k}e_{k}^{\dagger} + l_{k+1}e_{k+1}^{\dagger}.$$

因此$L_{k}^{-1}L_{k-1}^{-1}$是一个带有$L_{k}^{-1}$和$L_{k-1}^{-1}$的非零元素的单位下三角矩阵。$\blacksquare$