高斯消元法与 LU 分解

上一节中，我们从n元线性方程组开始引入了矩阵的概念。 并介绍了 Gauss-Jordan 消元法，该方法不仅可以求解线性方程组，同时还可以对输入的矩阵求逆。 一旦获得了矩阵\(\boldsymbol{A}\)的逆\(\boldsymbol{A}^{-1}\)，理论上我们就可以直接写出方程组 \(\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}\) 的解 \(\boldsymbol{x} = \boldsymbol{A}^{-1}\boldsymbol{b}\)。但是实际的工程应用中，很少这样做，主要是因为这种方法效率不高而且不稳定。 本节我们从高斯消元法开始，介绍效率更高用途更广的 LU 分解。

1. 高斯消元法

大多数线性代数的教科书中都是用高斯消元法来求解线性方程组的，这是一种介于 Gauss-Jordan 消元法和 LU 分解之间的方法。 相比于 Gauss-Jordan 消元法，需要把矩阵 \(\boldsymbol{A}\) 消减成一个单位矩阵，高斯消元法只消减出一个上三角的矩阵，所以其计算量相对较少一些。 Gauss-Jordan 消元法需要执行 \(O(n^3)\) 的乘法运算，而高斯消元法则最多需要 \(\frac{1}{3}O(n^3)\)。比如下面左侧一个 \(4 \times 4\) 的矩阵 \(\boldsymbol{A}\) 和向量\(\boldsymbol{b}\)构成的方程组，经过一系列的初等行变换就可以得到右侧的形式:

$$ \underbrace{\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \cdot \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}}_{\boldsymbol{A}\cdot \boldsymbol{x} = \boldsymbol{b}} \Longrightarrow \underbrace{\begin{bmatrix} a_{11}' & a_{12}' & a_{13}' & a_{14}' \\ 0 & a_{22}' & a_{23}' & a_{24}' \\ 0 & 0 & a_{33}' & a_{34}' \\ 0 & 0 & 0 & a_{44}' \end{bmatrix} \cdot \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} b_1' \\ b_2' \\ b_3' \\ b_4' \end{bmatrix}}_{\boldsymbol{A}'\cdot \boldsymbol{x} = \boldsymbol{b}'} $$

其中 \(x_1, x_2, x_3, x_4\) 是要求的四个未知数，我们可以直接写出它们的解。

$$ \begin{cases} x_1 = \frac{1}{a_{11}'}(b_1' - a_{14}'x_4 - a_{13}'x_3 - a_{12}'x_2) \\ x_2 = \frac{1}{a_{22}'}(b_2' - a_{24}'x_4 - a_{23}'x_3) \\ x_3 = \frac{1}{a_{33}'}(b_3' - a_{34}'x_4) \\ x_4 = b_4' / a_{44}' \end{cases} \Longrightarrow x_i = \frac{1}{a_{ii}'} \cdot \left[b_i' - \sum_{j=i+1}^4 a_{ij}'x_j \right] $$

利用上式，虽然我们只需要\(O(n^2)\)的乘法运算就可以直接写出方程的解，但是仍然需要 \(\frac{1}{3}O(n^3)\) 的乘法运算来进行初等行变换。对于同一个矩阵 \(\boldsymbol{A}\)， 不同的向量 \(\boldsymbol{b}\)，Gauss-Jordan 和高斯消元法都需要先经历一遍初等行变换的过程，不存在什么效率的。实际上初等行变换的过程相当于左乘一个矩阵。 比如对于上面的 \(4 \times 4\) 的矩阵 \(\boldsymbol{A}\)，我们可以左乘如下的矩阵，把第一列的 \(a_{21}, a_{31}, a_{41}\) 干掉。

$$ \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \\ -a_{21} & a_{11} & 0 & 0 \\ -a_{31} & 0 & a_{11} & 0 \\ -a_{41} & 0 & 0 & a_{11} \\ \end{bmatrix}}_{L_1} \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & b_1 \\ a_{21} & a_{22} & a_{23} & a_{24} & b_2 \\ a_{31} & a_{32} & a_{33} & a_{34} & b_3 \\ a_{41} & a_{42} & a_{43} & a_{44} & b_4 \end{bmatrix} \Longrightarrow \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & b_1 \\ 0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)} & b_2^{(1)} \\ 0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)} & b_3^{(1)} \\ 0 & a_{42}^{(1)} & a_{43}^{(1)} & a_{44}^{(1)} & b_4^{(1)} \end{bmatrix} $$

我们需要再左乘如下两个矩阵，就可以得到高斯消元法期望得到的上三角矩阵:

$$ \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -a_{43}^{(2)} & a_{33}^{(2)} \\ \end{bmatrix}}_{L_3} \underbrace{ \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -a_{32}^{(1)} & a_{22}^{(1)} & 0 \\ 0 & -a_{42}^{(1)} & 0 & a_{22}^{(1)} \\ \end{bmatrix}}_{L_2} \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & b_1 \\ 0 & a_{22}^{(1)} & a_{23}^{(1)} & a_{24}^{(1)} & b_2^{(1)} \\ 0 & a_{32}^{(1)} & a_{33}^{(1)} & a_{34}^{(1)} & b_3^{(1)} \\ 0 & a_{42}^{(1)} & a_{43}^{(1)} & a_{44}^{(1)} & b_4^{(1)} \end{bmatrix} \Longrightarrow \begin{bmatrix} a_{11}' & a_{12}' & a_{13}' & a_{14}' & b_1'\\ 0 & a_{22}' & a_{23}' & a_{24}' & b_2'\\ 0 & 0 & a_{33}' & a_{34}' & b_3'\\ 0 & 0 & 0 & a_{44}' & b_4' \end{bmatrix} $$

初等行变换的选择并不是唯一的，只要能达到消元的目的就行，在上面的例子中，我们刻意选择了下三角形式的三个变换矩阵 \(L_3,L_2,L_1\) 它们的乘积也是一个下三角矩阵。 原矩阵 \(A\) 就可以写成 \(LU\) 的形式，其中 \(L\) 和 \(U\) 分别是一个下三角和上三角矩阵。这就是一个 LU 的分解。

2. LU 分解

假设我们把矩阵\(\boldsymbol{A}\)写成了两个矩阵相乘的形式，其中 \(\boldsymbol{L}\) 是一个下三角矩阵， \(\boldsymbol{U}\)则是上三角矩阵。对于一个 \(4 \times 4\) 的矩阵\(\boldsymbol{A}\)，可以写成如下形式:

$$ \begin{equation}\label{equa_lua} \underbrace{\begin{bmatrix} \alpha_{11} & 0 & 0 & 0 \\ \alpha_{21} & \alpha_{22} & 0 & 0 \\ \alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\ \alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \end{bmatrix}}_{\boldsymbol{L}} \underbrace{\begin{bmatrix} \beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\ 0 & \beta_{22} & \beta_{23} & \beta_{24} \\ 0 & 0 & \beta_{33} & \beta_{34} \\ 0 & 0 & 0 & \beta_{44} \end{bmatrix}}_{\boldsymbol{U}} = \underbrace{\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}}_{\boldsymbol{A}} \end{equation} $$

那么线性方程组 \(\boldsymbol{A}\cdot \boldsymbol{x} = \boldsymbol{b}\) 就可以写成如下的形式:

$$ \begin{equation}\label{equa_axb} \boldsymbol{A}\cdot \boldsymbol{x} = \left(\boldsymbol{L}\boldsymbol{U}\right) \cdot \boldsymbol{x} = \boldsymbol{L}\left(\boldsymbol{U}\boldsymbol{x}\right) = \boldsymbol{b} \end{equation} $$

我们可以分两步来求解 \(\boldsymbol{x}\)，先求解式(\(\ref{equa_lyb}\))得到 \(\boldsymbol{y}\)，再求解式(\(\ref{equa_uxy}\))得到最终的解\(\boldsymbol{x}\)

$$ \begin{equation}\label{equa_lyb} \boldsymbol{L}\cdot \boldsymbol{y} = \boldsymbol{b} \end{equation} $$ $$ \begin{equation}\label{equa_uxy} \boldsymbol{U}\cdot \boldsymbol{x} = \boldsymbol{y} \end{equation} $$

由于\(\boldsymbol{L}\)是一个下三角矩阵，我们可以直接写出 \(\boldsymbol{y}\) 的解:

$$ \begin{equation}\label{solution_y} y_1 = \frac{b_1}{\alpha_{11}}, \qquad y_i = \frac{1}{\alpha_{ii}}\left[b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_j \right] \end{equation} $$

利用上三角矩阵\(\boldsymbol{U}\)，可以写出\(\boldsymbol{x}\)的解:

$$ \begin{equation}\label{solution_x} x_N = \frac{y_N}{\beta_{NN}}, \qquad x_i = \frac{1}{\beta_{ii}}\left[y_i - \sum_{j=1+1}^{N} \beta_{ij} y_j \right] \end{equation} $$

3. Crout's 算法

根据式(\(\ref{equa_lua}\))的关系，我们可以写出如下关于 \(\alpha_{ij}, \beta_{jk}\) 的 \(N \times N\) 个方程，但是\(\alpha_{ij}, \beta_{jk}\) 一共有 \(N^2 + N\) 个未知数， 所以为了保证有解，我们通常认定 \(\alpha_{ii} = 1\)，即矩阵 \(\boldsymbol{L}\) 的对角线全为 1。

$$ \begin{equation}\label{equa_crout} \left[ \begin{array}{l} \alpha_{11}\beta_{11} & \alpha_{11}\beta_{12} & \alpha_{11}\beta_{13} & \alpha_{11}\beta_{14} \\ \alpha_{21}\beta_{11} & \alpha_{21}\beta_{12} + \alpha_{22}\beta_{22} & \alpha_{21}\beta_{13} + \alpha_{22}\beta_{23} & \alpha_{21}\beta_{14} + \alpha_{22}\beta_{24} \\ \alpha_{31}\beta_{11} & \alpha_{31}\beta_{12} + \alpha_{32}\beta_{22} & \alpha_{31}\beta_{13} + \alpha_{32}\beta_{23} + \alpha_{33}\beta_{33} & \alpha_{31}\beta_{14} + \alpha_{32}\beta_{24} + \alpha_{33}\beta_{34} \\ \alpha_{41}\beta_{11} & \alpha_{41}\beta_{12} + \alpha_{42}\beta_{22} & \alpha_{41}\beta_{13} + \alpha_{42}\beta_{23} + \alpha_{43}\beta_{33} & \alpha_{41}\beta_{14} + \alpha_{42}\beta_{24} + \alpha_{43}\beta_{34} + \alpha_{44}\beta_{44} \end{array} \right] = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \end{equation} $$

如此一来，我们就可以逐列写出 \(\alpha_{ij}, \beta_{jk}\) 了。比如说第一列，有 \(\alpha_{11} = 1 \Rightarrow \beta_{11} = a_{11}\)，因此 \(\alpha_{i1} = a_{i1} / \beta_{11}\)。 再解出 \(\alpha_{i1}\) 的基础上，我们可以写出第二列中，\(\beta_{12} = a_{12}\)，\(\beta_{22} = a_{22} - \alpha_{21}\beta_{12}\)； \(\alpha_{22} = 1\)，\(\alpha_{i2} = (a_{i2} - \alpha_{i1}\beta_{12}) / \beta_{22}\)，这正是 Crout's 算法，它可以形式化描述如下左侧。 另外，由于 \(\alpha_{ii} = 1\)，在计算机里，我们可以把\(\boldsymbol{L},\boldsymbol{U}\)保存到一个 \(N\times N\) 的矩阵中，即下面右侧的形式。

• 对于\(L\)对角线 \(i \in [1, \cdots, N]\)，令 \(\alpha_{ii} = 1\)
• 对于每一列 \(j \in [1, \cdots, N]\)，先后完成如下两个计算:
  1. 对于\(\boldsymbol{U}\)中第 \(j\) 列中每一行 \(i \in [1, \cdots, j]\)，有 \(\beta_{ij} = a_{ij} - \sum_{k = 1}^{i} \alpha_{ik}\beta_{kj}\)
  2. 对于\(\boldsymbol{L}\)中第 \(j\) 列中每一行 \(i \in [j + 1, \cdots, N]\)，有 \(\alpha_{ij} = (a_{ij} - \sum_{k=1}^{j} \alpha_{ik}\beta_{kj})/\beta_{jj}\)

以上的所有讨论中，都忽略了一个问题，如果矩阵 \(\boldsymbol{A}\) 的对角线上存在元素 0，那么得到的 \(\boldsymbol{LU}\) 中将有一些非数字的元素。 这并不是我们期望看到的结果。为了规避这一问题，我们就需要逐列找一个最大值作为主元(pivot)，然后将主元所在行与当前行交换。 但是这样一来就改变了方程\(\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}\)中\(\boldsymbol{b}\)中的元素顺序，所以还需要额外记录下发生交换的行，在后续的求解过程中恢复出来。 这种带主元的 LU 分解也被称为 PLU。

4. 完

我们在例程中定义了一个模板类， 在它的构造函数中完成矩阵 \(\boldsymbol{A}\) 的带主元 LU分解，并保存到成员变量 mLU 中。同时提供了函数 Solve 和 Inverse 分别用于求解方程组和原矩阵的逆。