今天晚上突然被人提到一下,然后就想起来了。原来我还可以翻译姐姐的博客来着。这一篇来自这里,是介绍的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}$是一个复数方阵[^1],这表示这个系统有相同数量的方程和未知数。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, 删除,则可以正常运行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
-------------------         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$