矩阵的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消元法的基本思想是将矩阵$ $通过初等行变换转化为阶梯型矩阵，我们假定该过程中不涉及“对换变换”，变换过程如下：

记

\[ 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\)步，那么可以将矩阵$ $化为上三角矩阵：

\[ \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 \]

那么问题是，如何判断$ \(的前\)n-1\(个主元\)$？

回溯上述过程可知：

\[ \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\)即为$ \(的\)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\) 的逆，但是实际算法求解过程中，三角型的矩阵的乘法和逆是容易计算的，并且甚至不需要额外的空间存储开销 ：算法会将 \(\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分解 ，其伪代码如算法 \(\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 分解（对称三角分解）。