矩阵论：矩阵的LU分解

矩阵的LU分解

Gauss消元法

在介绍LU分解之前，我们先通过Gauss消元法引入。

对于求解以下线性方程组的问题：

$$
\left{\begin{matrix}
a_{11}x_{11} + a_{12}x_{12} + \cdots + a_{1n}x_{1n} = b_1\
\vdots \
a_{n1}x_{n1} + a_{n2}x_{n2} + \cdots + a_{nn}x_{nn} = b_n \
\end{matrix}\right.
$$

其矩阵形式为：

$$
\begin{bmatrix}
a_{11} & a_{12} & \cdots & a_{1n} & b_1\
\vdots & \vdots & \vdots & \vdots & \vdots\
a_{n1} & a_{n2} & \cdots & a_{nn} & b_n
\end{bmatrix} \triangleq (A \vdots b)
$$

Guass消元法的基本思想是将矩阵$ \boldsymbol{A}$通过初等行变换转化为阶梯型矩阵，我们假定该过程中不涉及“对换变换”，变换过程如下：

记

$$
A^{(1)} = A = \begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)}\
\vdots & \vdots & \vdots & \vdots\
a_{n1}^{(1)} & a_{n2}^{(1)} & \cdots & a_{nn}^{(1)}
\end{bmatrix}
$$

若$a_{11}^{(1)} \ne 0$，则：

$$
A^{(1)} \xrightarrow{r_i - \frac{a_{21}^{(1)}}{a_{11}^{(1)}}r_1,i = 2,3,\cdots,n}
\begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)}\
0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)}\
\vdots & \vdots & \vdots & \vdots\
0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)}
\end{bmatrix} \triangleq A^{(2)}
$$

若$a_{22}^{(2)} \ne 0$，则：

$$
A^{(2)} \xrightarrow{r_i - \frac{a_{32}^{(2)}}{a_{22}^{(2)}}r_1,i = 3,4,\cdots,n}
\begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)}& \cdots & a_{1n}^{(1)}\
& a_{22}^{(2)} & a_{23}^{(2)}& \cdots & a_{2n}^{(2)}\
& & a_{33}^{(3)} & \cdots & a_{3n}^{(3)}\
& & \vdots & & \vdots\
& & a_{n3}^{(3)} & \cdots & a_{nn}^{(3)}
\end{bmatrix} \triangleq A^{(3)}
$$

如此继续，如果能进行$n-1$步，那么可以将矩阵$ \boldsymbol{A}$化为上三角矩阵：

$$
\begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)}\
& a_{22}^{(2)} & \cdots & a_{2n}^{(2)}\
& & \ddots & \vdots\
& & & a_{nn}^{(n)}
\end{bmatrix} \triangleq A^{(n)}
$$

上述过程就是Guass消元法的矩阵描述。

在上述过程进行的过程中，易得，如果Guass消元法能进行到底：

$$
\Leftrightarrow a_{11}^{(1)}a_{22}^{(2)}\cdots a_{n-1n-1}^{(n-1)} \ne 0
$$

那么问题是，如何判断$ \boldsymbol{A}$的前$n-1$个主元$\ne 0$？

回溯上述过程可知：

$$
\begin{aligned}
& a_{11}^{(1)} = a_{11}\
& a_{11}^{(1)}a_{22}^{(2)} =
\begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)}\
0 & a_{22}^{(22)}
\end{bmatrix} =
\begin{bmatrix}
a_{11} & a_{12}\
a_{21} & a_{22}
\end{bmatrix}\
& \cdots \
& a_{11}^{(1)}a_{22}^{(2)}\cdots a_{kk}^{(k)} =
\begin{bmatrix}
a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1k}^{(1)}\
& a_{22}^{(2)} & \cdots & a_{2k}^{(2)}\
& & \ddots & \vdots\
& & & a_{kk}^{(k)}
\end{bmatrix} =
\begin{bmatrix}
a_{11} & \cdots & a_{1k}\
\vdots & & \vdots\
a_{k1} & \cdots & a_{kk}
\end{bmatrix}
\end{aligned}
$$

记：

$$
\begin{aligned}
& \Delta_1 = a_{11} \
& \Delta_2 = \begin{bmatrix}
a_{11} & a_{12}\
a_{21} & a_{22}
\end{bmatrix} \
& \cdots \
& \Delta_{n-1} = \begin{bmatrix}
a_{11} & \cdots & a_{1n-1}\
\vdots & & \vdots\
a_{n-11} & \cdots & a_{n-1n-1}
\end{bmatrix}
\end{aligned}
$$

可以发现，$\Delta_k$即为$ \boldsymbol{A}$的$k$阶顺序主子式。

所以，得到如下命题：

命题 1
当$\Delta_k \ne 0,k = 1,2,\cdots,n-1\Leftrightarrow a_{kk}^{(k)} \ne 0$

于是，自然得到：

定理 1
Guass消元法能进行到底 $\Leftrightarrow$ $\boldsymbol{A}$的前$n-1$个顺序主子式$\ne 0$；

矩阵LU分解的步骤推导

从高斯消元法的过程，每一步初等行变换的过程相当于一个高斯变换，即左乘一个矩阵：

$$
L_1 = \begin{bmatrix}
1 & & & \
l_1 & 1 & &\
l_2 & & 1 & & \
\vdots & & & \ddots & \
l_n & & & & 1
\end{bmatrix}
$$

其中，$l_i = -\frac{a_{i1}^{(1)}}{a_{11}^{(1)}},i = 2,3,\cdots,n$。

易得，高斯消元法的过程可以表示为：

$$
L_{n-1}L_{n-2}\cdots L_{1}A = U
$$

其中，$\boldsymbol{L}_i$为下三角矩阵，$\boldsymbol{U}$为上三角矩阵。

令 $\boldsymbol{L} = L_1^{-1}\cdots L_{n-2}^{-1}L_{n-1}^{-1}$，则有：

$$
\boldsymbol{A} = \boldsymbol{L} \boldsymbol{U}
$$

这就是矩阵的LU分解。

求得 $\boldsymbol{L}$ 和 $\boldsymbol{U}$ 后，如何计算 $x$？

$$
Ax = LUx = b
$$

取中间变量 $y = Ux$，则 $Ly = b$。

那么，为什么LU分解能加速线性方程组的求解呢？

实际上，虽然从矩阵形式上看，$\boldsymbol{L}$ 需要计算 $n-1$ 个 $L_i$ 的逆，但是实际算法求解过程中，三角型的矩阵的乘法和逆是容易计算的，并且甚至不需要额外的空间存储开销 \cite{Matrix_Computations}：算法会将 $\boldsymbol{L}$ 和 $\boldsymbol{U}$ 的元素直接存储到原矩阵 $\boldsymbol{A}$ 的相应位置上：

• $\boldsymbol{L}$ 的下三角部分（不包括对角线）会存储在 $\boldsymbol{A}$ 的下三角部分。
• $\boldsymbol{U}$ 的上三角部分和对角线会存储在 $\boldsymbol{A}$ 的上三角部分和对角线上。

例如，对于矩阵

$$
\boldsymbol{A} =\begin{bmatrix}
a_{11} & a_{12} & a_{13} \
a_{21} & a_{22} & a_{23} \
a_{31} & a_{32} & a_{33}
\end{bmatrix}
$$

经过 LU 分解后，它的存储会变为：

$$
\boldsymbol{A} =\begin{bmatrix}
u_{11} & u_{12} & u_{13} \
l_{21} & u_{22} & u_{23} \
l_{31} & l_{32} & u_{33}
\end{bmatrix}
$$

其中：

• $l_{ij}$ 为第 $i$ 行，第 $j$ 列的乘子；
• $u_{ij}$ 为第 $i$ 行，第 $j$ 列的上三角元素；

通过这种方式，所有分解信息都存储在原矩阵中，无需额外空间。LU分解既降低了时间复杂度又降低了空间复杂度。

不过，从我们进行高斯消元法的过程就可以发现，要保证高斯消元法能进行到底，需要避免乘子中的除数（$A_{kk}$）出现 0（实际应用中，很小的情况也应该避免），所以一般不直接应用LU分解，而是应用列主元PLU分解，它通过交换行来避免出现相对很小的主元，实现稳定的LU分解 \cite{数值分析5}，其伪代码如算法 \ref{列选主元PLU分解法} 所示。

算法：列选主元PLU分解法

输入： 矩阵 $A \in \mathbb{R}^{n \times n}$
输出： 矩阵 $P, L, U$，使得 $P A = L U$

初始化：
U = A, L = I, P = I

for k = 1 to n-1 do:
    1. 寻找主元：
       pivot = max_{i=k, ..., n} |U[i, k]|
       记录主元行号 row
       
    2. 若 row ≠ k:
       构造置换矩阵 P_k，交换 U 的第 k 行与第 row 行：
       P_k = 交换 k 行与 row 行的单位矩阵
       
       更新 U: U = P_k U
       更新 P: P = P_k P
       
       若 k > 1，则更新 L: L = P_k L
    
    3. 对于 i = k+1 到 n:
       计算消元因子：
       L[i, k] = U[i, k] / U[k, k]
       
       更新 U：
       U[i, :] = U[i, :] - L[i, k] * U[k, :]
       
返回 P, L, U

$$
\boldsymbol{A} = \boldsymbol{L} \boldsymbol{U}
$$

其中，若矩阵 $\boldsymbol{A}$ 满秩，即 $\operatorname{rank} A = r$，则可利用 LU 分解或 PLU 分解构造其满秩分解：

• LU 分解构造满秩分解

若 $\boldsymbol{A}$ 可分解为 $\boldsymbol{A} = \boldsymbol{L} \boldsymbol{U}$，则从 $\boldsymbol{L}$ 提取前 $r$ 列，从 $\boldsymbol{U}$ 提取前 $r$ 行：
$$
\boldsymbol{F} = \boldsymbol{L}[:, 0:r],\quad \boldsymbol{G} = \boldsymbol{U}[0:r, :]
$$

其中，$\boldsymbol{F}$ 为列满秩矩阵（$m \times r$），$\boldsymbol{G}$ 为行满秩矩阵（$r \times n$），从而得到矩阵的满秩分解。

• PLU 分解构造满秩分解

设 $\boldsymbol{A}$ 通过列主元 LU 分解得：

$$
\boldsymbol{A} = \boldsymbol{P} \boldsymbol{L} \boldsymbol{U}
$$

同样取：
$$
\boldsymbol{B} = \boldsymbol{L}[:, 0:r],\quad \boldsymbol{C} = \boldsymbol{U}[0:r, :]
$$

代入分解式可得：
$$
\boldsymbol{A} = \boldsymbol{P} \boldsymbol{B} \boldsymbol{C}
$$

其中，$\boldsymbol{P}$ 为置换矩阵，显然满秩。

进一步推广，若方阵 $\boldsymbol{A}$ 可分解为：
$$
\boldsymbol{A} = \boldsymbol{L} \boldsymbol{D} \boldsymbol{U}
$$

其中，$\boldsymbol{D}$ 为对角矩阵，则称其为 LDU 分解。

由高斯消元法的推导，可得如下定理及推论：

定理：矩阵 $\boldsymbol{A} = (a_{i j})_{n \times n}$ 存在唯一分解：
$$
\boldsymbol{A} = \boldsymbol{L} \boldsymbol{D} \boldsymbol{U}
$$

其中：

• $\boldsymbol{L}$ 为单位下三角矩阵，

• $\boldsymbol{U}$ 为单位上三角矩阵，

• $\boldsymbol{D}$ 为对角矩阵，其对角元：$d_{k} = \frac{\Delta_{k}}{\Delta_{k-1}}, \quad k=1,2,\dots,n \quad (\Delta_0 = 1)$

其中 $\Delta_k$ 为 $\boldsymbol{A}$ 的顺序主子式。

推论：若 $\boldsymbol{A}$ 为非奇异矩阵，则其可进行三角分解：
$$
\boldsymbol{A} = \boldsymbol{L} \boldsymbol{U}
$$

当且仅当 $\boldsymbol{A}$ 的顺序主子式：$\Delta_k \neq 0, \quad k=1,2,\dots,n$

进一步推广：对 Doolittle、Crout、Cholesky 分解，定义如下：

• Doolittle 分解：设 $\boldsymbol{A}$ 具有唯一 LDU 分解，令：
$$
\hat{\boldsymbol{U}} = \boldsymbol{D} \boldsymbol{U}
$$

则有：
$$
\boldsymbol{A} = \boldsymbol{L} \hat{\boldsymbol{U}}
$$

• Crout 分解：令：
$$
\hat{\boldsymbol{L}} = \boldsymbol{L} \boldsymbol{D}
$$

则有：
$$
\boldsymbol{A} = \hat{\boldsymbol{L}} \boldsymbol{U}
$$

• Cholesky 分解（平方根分解）：设 $\boldsymbol{A}$ 为实对称正定矩阵，则其可分解为：
$$
\boldsymbol{A} = \boldsymbol{L} \widetilde{\boldsymbol{D}}^2 \boldsymbol{U}
$$

其中：
$$
\widetilde{\boldsymbol{D}} = \operatorname{diag}(\sqrt{d_1},\sqrt{d_2},\dots,\sqrt{d_n})
$$

进一步，可得：
$$
\boldsymbol{A} = \boldsymbol{L} \widetilde{\boldsymbol{D}}^2 \boldsymbol{L}^\top
$$

令：
$$
\boldsymbol{G} = \boldsymbol{L} \widetilde{\boldsymbol{D}}
$$

则有：
$$
\boldsymbol{A} = \boldsymbol{G} \boldsymbol{G}^\top
$$

这被称为 Cholesky 分解（对称三角分解）。