Direct Methods for Solving Linear Systems

本章很多内容在数值代数中会更加详细

Linear Systems of Equations

介绍了高斯消元法的过程。高斯消元法的原理就是化上三角 (upper-triangular) 矩阵，或称阶梯形矩阵，随后使用 backward-substitution 过程完成解方程（解上三角矩阵的方法）。

Upper-Triangular Matrices

首先化上三角矩阵。第 \(k\) 步处理右下角的一个 \((n-k+1)\times (n-k+1)\) 的方阵：

需要此时的 \(a_{kk}\neq 0\) （第 \(k\) 步的 \(a_{kk}\) 常写作 \(a_{kk}^{(k)}\)）
计算 \(m_{ik} = a_{ik} / a_{kk}\), \(i=k+1, k+2, \cdots, n\)

每一行要减去时所乘的倍数
更新方阵元素，即有

\[ \begin{aligned} a_{ij}^{(k+1)} &= a_{ij}^{(k)} - m_{ik}a_{kj}^{(k)}\\ b_{ij}^{(k+1)} &= b_{ij}^{(k)} - m_{ik}b_{kj}^{(k)} \end{aligned} \]

\(i, j=k+1, k+2, \cdots, n\)。其中 \(a_{kk}\) 下方第 \(k\) 列的元素已经不需要再考虑了，因为它们已经被消成 0 了。实际进行 LU 分解时，这里用来存储的就是 \(m_{ik}\)，即保存的 \(L\)。

化上三角形矩阵的过程中，总共需要进行 \(n-1\) 步。

Backward-Substitution

Backward-Substitution 的过程就是解系数矩阵为上三角矩阵的线性方程组的过程，即

\[ x_n = \frac{b_n}{a_{nn}},\quad x_i = \frac{b_i - \sum\limits_{j=i+1}^n a_{ij}x_j}{a_{ii}},\quad i=n-1, n-2, \cdots, 1 \]

Amount of Computation

追溯其乘法 / 除法计算量，先考虑比较复杂的化上三角过程。

在第 \(k\) 步，计算 \(m_{ik}\) 需要 \((n-k)\) 次除法，更新 \(a_{ij}\) 和 \(b_i\) 需要考虑 \(i\) 从 \(k+1\) 到 \(n\) 共 \((n-k)\) 行，每行计算 \(n-k\) 个 \(a_{ij}\) 和 \(1\) 个 \(b_i\)（每个都需要 \(1\) 次乘法）。因此，第 \(k\) 步需要的乘除法数量为

\[ (n-k) + (n-k)(n-k+1) = (n-k)(n-k+2) \]

因此，化上三角的过程需要的乘除法数量为

\[ \sum_{k=1}^{n-1}(n-k)(n-k+2) = \frac{n^3}{3} + \frac{n^2}{2} + \frac{5}{6}n \]

而 backward-substitution 过程的乘除法数量就很好算了：

\[ 1 + \sum_{i=1}^{n-1}(n-i+1) = \frac{n^2}{2} + \frac{n}{2} \]

总的乘除法次数为

\[ \frac{n^3}{3} + n^2 + \frac{4}{3}n \]

Pivoting Strategies

Trouble Maker

使用原始的高斯消元法尝试解决

\[ \begin{cases} 0.003000x_1 + 59.14x_2 = 59.17 \\ 5.291x_1 - 6.130x_2 = 46.78 \end{cases} \]

每一步计算都使用 4 位有效数字的四舍五入（使用 matlab 的 round 函数实现），得到的结果是

\[ x_1=-10.0000, \quad x_2=1.0010 \]

精确解本应是

\[ x_1=10, \quad x_2=1 \]

经过分析，在使用 \(4\) 位有效数字的四舍五入的情况下

\[ L = \begin{bmatrix} 1 & 0\\ 1764 & 1 \end{bmatrix} \quad U = \begin{bmatrix} 0.0030 & 59.1400\\ 0 & -104300 \end{bmatrix} \]

这是因为在第一步消元时，\(a_{21}\) 的绝对值比 \(a_{11}\) 的绝对值大了 3 个数量级，在进行除法运算时，造成了较大的误差

为了解决这个问题，需要引入选主元策略 (Pivoting Strategy)

根据高斯消元法的过程，主元 \(a_{kk}^{(k-1)}\) 会作为除法的分母
如果主元太小，舍入误差会很大
因此，为了减少舍入误差，我们需要使 \(a_{kk}^{(k-1)}\) 的值更大。

Complete Pivoting

首先介绍原始的全主元 (Complete Pivoting) 策略，其本质是通过一系列行交换矩阵的复合 \(P\) 和列交换矩阵 \(Q\)，将问题转化为

\[ PAQ(Q^{-1}x)=Pb \]

这样，计算 \(PAQ\) 的 LU 分解 \(PAQ=LU\)，然后再依次解

\[ Ly=Pb,\quad Uz=y \]

最后就可以得到 \(x=Qz\)。在全主元高斯消元的第 \(k\) 轮，选择

\[ |A(p, q)| = \max_{k\leqslant i\leqslant n,\; k\leqslant j\leqslant n}|A(i, j)| \]

使它到达对角元位置 \(A(k, k)\)，这样就能保证除法的分母在有意义的范围（\(A\) 中仍需要参与 \(LU\) 分解计算的右下方子矩阵）内最大，从而尽可能地减小了舍入误差。

(Scaled) Partial Pivoting

全主元策略选主元的开销是比较大的，作为一种妥协引入了列主元 (Partial Pivoting, Maximal Column Pivoting) 策略，即只进行一系列行交换

\[ PAx=Pb \]

这样，计算 \(PA\) 的 LU 分解 \(PA=LU\)，然后再依次解

\[ Ly=Pb,\quad Ux=y \]

即可解得 \(x\)。而 Scaled Partial Pivoting 策略相比列主元，只是在选择列主元的过程中考虑元素相对行最大值的相对大小(\(\leqslant 1\)) 而不是绝对大小而已

关键的选择上，对于列主元有

\[ |A(p, k)| = \max_{k\leqslant i\leqslant n}|A(i, k)| \]

而对于 Scaled 列主元有

\[ \begin{gathered} s_i = \max_{1\leqslant j\leqslant n}|a_{ij}|\\ \frac{|A(p, k)|}{s_p} = \max_{k\leqslant i\leqslant n}\frac{|A(i, k)|}{s_i} \end{gathered} \]

Scaled Partial Pivoting 相比 Partial Pivoting，预先完成所有缩放因子 \(s_i\) 的计算，随后根据相对缩放因子的大小选择列主元，使 \(\frac{a_{kk}^{(k-1)}}{s_k}\) 为最大值

Remark

在 Partial Pivoting(PP) 策略中，每一列中选取的主元元素是主元以下该列中的最大元素，这种方法可能导致误差放大的程度比 Scaled Partial Pivoting(SPP) 更大
另一方面，SPP 通过将每一列中的元素按照其绝对值的最大值进行缩放来选择主元元素，有助于减小舍入误差的影响，提高数值稳定性
在具有特定的结构的数据下，SPP 数值稳定性将超过 PP，例如

\[ A = \begin{bmatrix} 1 & 10000\\ 1 & 0.0001 \end{bmatrix} ,\quad b = \begin{bmatrix} 10000\\ 1 \end{bmatrix} \]

在 4 位机器上运算，PP 将不会进行行交换，解得数值解 \((0, 1)^{\top}\)，而 SPP 会进行行交换，得到更精确的解 \((1, 1)^{\top}\)。

Amount of Computation

Partial Pivoting: 需要 \(O(n^2)\) 的额外比较
Scaled Partial Pivoting: 需要 \(O(n^2)\) 的额外比较和 \(O(n^2)\) 的额外除法
Complete Pivoting: 需要 \(O(n^3/3)\) 的额外比较

然而，如果 Scaled Partial Pivoting 中的 scaled factors 不是预先计算好的，而是会实时更新，那么还需要 \(O(n^3/3)\) 的额外除法。

Matrix Factorization

实际上在前面的分析中已经默认采用了这里的 LU 分解的视角

要解线性方程组

\[ Ax=b \]

常用的解法是计算 \(A\) 的 LU 分解，即 \(A=LU\)，其中 \(L\) 为下三角矩阵，\(U\) 为上三角矩阵。这样，原方程组 \(Ax=b\) 就可以转化为 \(L(Ux)=b\)，于是可以先解

\[ Ly=b \]

得到 \(y\)，再解

\[ Ux=y \]

就可以得到原方程组的解了。LU 分解的计算可以通过高斯消元，Cholesky 分解等方法进行。

在高斯消元法手算过程中，高斯消元是通过行交换、行倍加、行倍乘三种操作，最终得到阶梯形矩阵

从 LU 分解的角度来看，\(U\) 就是阶梯形矩阵，\(L\) 则是三种基本行变换矩阵的复合，由此可知 \(L\) 对角线元全为 1

因此在算法实现上，为了节约空间开销，可以直接使用原来的 \(A\) 存储 \(LU\) 分解的结果

\[ \left[ \begin{array}{ccccc} u_{11} & u_{12} & u_{13} & \cdots & u_{1 n} \\ l_{21} & u_{22} & u_{23} & \cdots & u_{2 n} \\ l_{31} & l_{32} & u_{33} & \cdots & u_{3 n} \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{n, n-1} & u_{n n} \end{array} \right] \]

值得一提的是，如果 \(U\) 的对角元全为 1（称此时 \(U\) 是 unitary 的），那么这种 LU 分解就称为 Crout's factorization。Crout's factorization 可以通过对 \(A^{\top}\) 进行 LU 分解得到。

Special Types of Matrices

Mathematical Preliminaries

Strictly Diagonally Dominant

对于 \(A=(a_{ij})_{n\times n}\)，如果满足以下条件，就称它是严格对角占优 (strictly diagonally dominant) 的

\[ |a_{ii}| > \sum\limits_{j\neq i}|a_{ij}| \]

对角线元素的绝对值比这一行所有其他元素的绝对值之和还要大。

特别地，\(>\) 改成 \(\geqslant\) 就是对角占优 (diagonally dominant) 的

Properties of Strictly Diagonally Dominant Matrices

严格对角占优的矩阵是非奇异的
高斯消元过程中不需要进行行 / 列交换
高斯消元过程的计算误差相对于舍入误差而言是稳定的

Proof

等价于证明 \(Ax=0\) 只有零解。

假设有非零解 \(x\)，取能令 \(x\) 取到无穷范数的下标 \(k\)，即

\[ x_k=\max_{1\leqslant i\leqslant n}|x_i| \]

然后考虑 \(Ax=0\) 的第 \(k\) 个方程，有

\[ \sum_{j=1}^n a_{kj}x_j = 0 \Rightarrow a_{kk} = - \sum_{j\neq k}a_{kj}\frac{x_j}{x_k} \]

从而有

\[ |a_{kk}| = \left|\sum_{j\neq k}a_{kj}\frac{x_j}{x_k}\right|\leqslant \sum_{j\neq k}|a_{kj}| \frac{|x_j|}{|x_k|} \leqslant \sum_{j\neq k}|a_{kj}| \]

与严格对角占优矛盾，得证。

关键在于证明 \(a_{kk}^{(k)}\neq 0\)
留做习题证明略

Properties of Positive Definite Matrices

设 \(A\) 是 \(n\times n\) 正定矩阵，则

\(A\) 是非奇异的（\(A^{-1}\) 存在）
\(\forall\: i=1, \cdots, n\), \(a_{ii}>0\)
\(\max\limits_{1\leqslant k, j\leqslant n}a_{kj}\leqslant \max\limits_{1\leqslant i\leqslant n}a_{ii}\)
\(\forall\: i\neq j\), \((a_{ij})^2< a_{ii}a_{jj}\)

留做习题证明略

Factorization of Symmetric Positive Definite Matrices

可以证明，若 \(A\) 是对称正定矩阵，则存在以下两种实际上等价的分解：

\(LDL^{\top}\) 分解：\(A=LDL^{\top}\)，其中 \(L\) 是对角线全为 1 的下三角矩阵 (unitary lower-triangular matrix)，\(D\) 是对角矩阵
Cholesky 分解：\(A=LL^{\top}\)，其中 \(L\) 是下三角矩阵

\(LDL^{\top}\) 分解的乘除法计算量为

\[ \frac{1}{6}n^3 + n^2 - \frac{7}{6}n \]

而 Cholesky 分解的乘除法计算量为

\[ \frac{1}{3}n^3 + \frac{1}{2}n^2 - \frac{2}{3}n \]

看似 Cholesky 分解占据优势，实际上 Cholesky 分解还需要计算 \(n\) 次平方根。不过这个计算是 \(O(n)\) 的，只是常数相对乘除法比较大而已。

Factorization of Tridiagonal Matrices

Crout's factorization 可以用于三对角矩阵的分解，即 \(A\) 是三对角矩阵，\(A=LU\)，其中 \(L\) 是下三角矩阵，\(U\) 是对角线全为 1 的上三角矩阵。Crout's factorization 的计算量为 \(O(n)\)。

\[ L=\begin{bmatrix} l_{11} \\ l_{21} & l_{22} \\ & \ddots & \ddots \\ & & l_{n-1, n-1} & l_{n-1, n} \\ & & & l_{nn} \end{bmatrix},\quad U=\begin{bmatrix} 1 & u_{12} \\ & 1 & u_{23} \\ & & \ddots & \ddots \\ & & & 1 & u_{n-1, n} \\ & & & & 1 \end{bmatrix} \]

先解 \(Ly=b\)，再解 \(Ux=y\)。

\[ \begin{aligned} y_1 &= \frac{b_1}{l_{11}}\\ y_i &= \frac{b_i - l_{i, i-1}y_{i-1}}{l_{ii}},\quad i=2, 3, \cdots, n\\ x_n &= y_n\\ x_i &= y_i - u_{i, i+1}x_{i+1},\quad i=n-1, n-2, \cdots, 1 \end{aligned} \]

可以注意到 Crout's factorization 要求 \(l_{ii}\neq 0\)。以下有三个条件可以保证这一点：

\(A\) 是对称正定的
\(A\) 是严格对角占优的
\(a_{i, i-1}\neq 0\), \(a_{i, i+1}\neq 0\)，\(|a_{11}|>|a_{12}|\), \(|a_{nn}|>|a_{n, n-1}|\), \(|a_{ii}|\geqslant |a_{i, i-1}| + |a_{i, i+1}|\)

可以注意到，第三个条件实际上就是比第二个条件多了次对角线全部非零，以及中间的行从严格对角占优退化为对角占优。