Chap 6: Direct Methods for Solving Linear Systems⚓︎

目标：求解 \(A \bm{x} = \bm{b}\)

Linear Systems for Equations⚓︎

Gaussian Elimination⚓︎

高斯消元法(Gaussian elimination) 的基本思路：

• 首先将矩阵 \(A\) 归约成一个上三角(upper-triangular) 矩阵
• 然后通过回代(backward-substitution) 来求解未知量

先来看消元(elimination) 的实现：先令 \(A^{(1)} = A = (a_{ij}^{(1)})_{n \times n}, \bm{b}^{(1)} = \bm{b} = \begin{bmatrix}b_1^{(1)} \\ \vdots \\ b_n^{(1)}\end{bmatrix}\)

• 第 1 步：

  • 如果 \(a_{11}^{(1)} \ne 0\)，计算 \(m_{i1} = \dfrac{a_{i1}^{(1)}}{a_{11}^{(1)}}, (i = 2, \dots, n)\)

  • 那么增广矩阵 (augmented matrix) 的第 \(i\) 行 \(\text{row}_i\) 为：\(m_{i1} \times \text{row}_1\)，得到

    \[ \left[ \begin{array}{cccc|c} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)} \\ 0 & & A^{(2)} & & \bm{b}^{(2)} \\ \end{array} \right] \]

    其中 \(\begin{cases}a_{ij}^{(2)} = a_{ij}^{(1)} - m_{i1} a_{1j}^{1} \\ b_i^{(2)} = b_i^{(1)} - m_{i1} b_1^{(1)}\end{cases}, (i, j = 2, \dots, n)\)

• 第 k 步：

  • 如果 \(a_{kk}^{(k)} \ne 0\)，计算 \(m_{ik} = \dfrac{a_{ik}^{(k)}}{a_{kk}^{(k)}}, (i = k+1, \dots, n)\)
  • \(\begin{cases}a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik} a_{kj}^{k} \\ b_i^{(k+1)} = b_i^{(k)} - m_{ik} b_k^{(k)}\end{cases}, (i, j = k+1, \dots, n)\)

• n-1 步后：

  \[ \begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \dots & a_{1n}^{(1)} \\ & a_{22}^{(2)} & \dots & a_{2n}^{(2)} \\ & & \dots & \vdots \\ & & & a_{nn}^{(n)}\end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_n\end{bmatrix} = \begin{bmatrix}b_1^{(1)} \\ b_2^{(2)} \\ \vdots \\ b_n^{(n)}\end{bmatrix} \]

接下来看回代：

• \(x_n = \dfrac{b_n^{(n)}}{a_{nn}^{(n)}}\)
• \(x_i = \dfrac{b_i^{(i)} - \sum\limits_{j=i+1}^n a_{ij}^{(i)} x_j}{a_{ii}^{(i)}}, (i = n - 1, \dots, 1)\)
• 我们必须找到最小的整数 \(k \ge i\) 且 \(a_{ki}^{(i)} \ne 0\)，然后交换第 \(k\) 行和第 \(i\) 行

Step 1  for i = 1, ..., n - 1 do Steps 2-4:
    Step 2  Let p be the smallest integer with i <= p <= n and a[p][i] != 0;
            if no integer p can be found
                then Output('no unique solution exists');
                STOP;
            // row exchange
    Step 3  if p != i then perform (E[p]) <-> (E[i]);
            // elimination
    Step 4  for j = i + 1, ..., n do Step 5 and 6:
        Step 5  Set m[j][i] = a[j][i] / a[i][i]
        Step 6  Perform (E[j] - m[j][i] * E[i]) -> (E[i])
    Step 7  if a[n][n] = 0 then Output('no unique solution exists');
                STOP;
            // backward substitution
    Step 8  Set x_n = a[n][n+1] / a[n][n]
    Step 9  for i = n - 1, ..., 1 set x[i] = [a[i][n+1] - sum(j=i+1, n, a[i][j] * x[j])] / a[i][i];
    Step 10 Output(x[1], ..., x[n]);
            STOP;    // success

Amount of Computation⚓︎

现在我们来统计一下计算量（仅考虑乘法 / 除法）。

• 消元：\(\sum\limits_{k=1}^{n-1} (n - k)(n - k + 2) = \dfrac{n^3}{3} + \dfrac{n^2}{2} - \dfrac{5}{6}n\)
• 回代：\(1 + \sum\limits_{i=1}^{n-1}(n - i + 1) = \dfrac{n^2}{2} + \dfrac{n}{2}\)

所以对于很大的 \(n\)，乘法和除法的总数大约为 \(\textcolor{red}{\dfrac{n^3}{3}}\)。也就是说，高斯消元法的时间复杂度为 \(O(n^3)\)。

Pivoting Strategies⚓︎

除非有特殊说明，以下的 \(k\) 指的是第 \(k\) 次高斯消元。

所以，我们需要选取合适的主元以减小误差。下面给出一些选取主元的策略。

Partial Pivoting⚓︎

部分主元法(partial pivoting)（或称为最大列主元法 (maximal column pivoting)）：找到最小的 \(p\)，使得 \(|a_{pk}^{(k)}| = \max\limits_{k \le i \le n} |a_{ik}^{(k)}|\)，然后交换第 \(p\) 行和第 \(k\) 行。

Step 1  for i = 1, ..., n set NROW(i) = i;    // initialize row pointer
Step 2  for i = 1, ..., n - 1 do Steps 3-6:    // elimination process
    Step 3  Let p be the smallest integer with i <= p <= n and
            |a(NROW(p), i)| = max_{i <= j <= n}|a(NROW(j), i)|;
    Step 4  if a(NROW(p), i) = 0 then Output('no unique solution exists');
                STOP;
    Step 5  if NROW(i) != NROW(p) then set NCOPY = NROW(i);    // simulated row interchange
                                          NROW(i) = NROW(p);
                                          NROW(p) = NCOPY;
    Step 6  for j = i + 1, ..., n do Steps 7-8:
        Step 7  Set m(NROW(j), i) = a(NROW(j), i) / a(NROW(i), i);
        Step 8  Perform (E_NROW(j) - m(NROW(j), i) * E_NROW(i)) -> (E_NROW(j));
    Step 9  if a(NROW(n), n) = 0 then Output('no unique solution exists');
                STOP;
    // start backward substitution
    Step 10 Set x[n] = a(NROW(n), n + 1) / a(NROW(n), n);
    Step 11 for i = n - 1, ..., 1
                set x[i] = (a(NROW(i), n + 1) - sum(j=i+1, n, a(NROW(i), j) * x[j])) / a(NROW(i), i);
    Step 12 Output(x[1], ..., x[n]);    // procedure completed successfully
            STOP;

Scaled Partial Pivoting⚓︎

缩放部分主元法(scaled partial pivoting)（或称为缩放列主元法 (scaled-column pivoting)）：将一行中最大的元素放在主元的位置上。

• 第 1 步：定义每一行的缩放因子(scale factor) \(s_i = \max\limits_{1 \le j \le n} |a_{ij}|\)（即每行中绝对值最大的元素）
• 第 2 步：（对于第 \(k\) 次高斯消元，）找到最小的 \(p \ge k\)，使得 \(\dfrac{|a_{pk}^{(k)}|}{s_p} = \max\limits_{k \le i \le n} \dfrac{|a_{ik}^{(k)}|}{s_i}\)，然后交换第 \(p\) 行和第 \(k\) 行

Step 1  for i = 1, ..., n set s_i = max_{1 <= j <= n}(|a_ij|);
                          if s_i = 0 then Output('no unique solution exists');
                                          STOP;
Step 2  for i = 1, ..., n - 1 do Steps 3-6:    // elimination process
    Step 3  Let p be the smallest integer with i <= p <= n and
            |a(NROW(p), i)| / s(NROW(p)) = max_{i <= j <= n}(|a(NROW(j), i)| / s(NROW(j)));

Complete Pivoting⚓︎

完全主元法(complete pivoting)（或称为最大主元法 (maximal pivoting)）：搜索所有的元素 \(a_{ij} (i, j = k, \dots, n)\)，找出其中数值最大的元素。通过互换(interchange) 行和列，使得该元素来到主元的位置上。

Amount of Computation⚓︎

• 部分主元法：需要 \(O(n^2)\) 次比较
• 缩放部分主元法：需要 \(O(n^2)\) 次比较，以及 \(O(n^2)\) 次除法（初始计算）
• 完全主元法：需要 \(O(\dfrac{n^3}{3})\) 次比较

所以，要想追求更好的稳定性，就要以更大的计算量为代价。

Matrix Factorization⚓︎

高斯消元法是一种简单粗暴的方法，但效率不是很高。因此下面介绍一种基于高斯消元法实现的改进方法——矩阵分解(matrix factorization)。它的计算过程如下：

• 第 1 步：

  • \(m_{i1} = \dfrac{a_{i1}}{a_{11}} (a_{11} \ne 0)\)
  • 令 \(L_1 = \begin{bmatrix}1 & & & \\ -m_{21} & 1 & & \\ \vdots & & \ddots \\ -m_{n1} & & &1\end{bmatrix}\)（这就是第一高斯变换矩阵(first Gaussian transformation matrix)），那么 \(L_1 [A^{(1)} \quad \bm{b}^{(1)}] = \begin{bmatrix}a_{11}^{(1)} & \dots a_{1n}^{(1)} & b_1^{(1)} \\ 0 & A^{(2)} & \bm{b}^{(2)}\end{bmatrix}\)

• 第 k 步：

  • 第 k 高斯变换矩阵(kth Gaussian transformation matrix)：\(L_k = \begin{bmatrix}1 & & & & \\ & \ddots & & & \\ & & 1 & & \\ & & -m_{k+1, k} & & \\ & & \vdots &\ddots & \\ & & -m_{n, k} & & 1\end{bmatrix}\)（空的地方都是 0）

• 第 n-1 步：

  \[ L_{n-1}L_{n-2} \dots L_1 [A \quad \bm{b}] = \begin{bmatrix}a_{11}^{(1)} & a_{12}^{(1)} & \dots & a_{1n}^{(1)} & b_1^{(1)}\\ & a_{22}^{(2)} & \dots & a_{2n}^{(2)} & b_2^{(2)} \\ & & \dots & \vdots & \vdots\\ & & & a_{nn}^{(n)} & b_n^{(n)}\end{bmatrix} \]

Step 1  Select l_11 and u_11 satisfying l_11 * u_11 = a_11;
        if l_11 * u_11 = 0 then Output('Factorization impossible');
            Stop;
Step 2  for j = 2, ..., n set u_1j = a_1j / l_11;    // first row of U
                              l_j1 = a_j1 / u_11;    // first column of L
Step 3  for i = 2, ..., n - 1 do Steps 4 and 5:
    Step 4  Select l_ii and u_ii safisfying l_ii * u_ii = a_ii - sum(k=1, i-1, l_ik * u_ki);
            if l_ii * u_ii = 0 then Output('Factorization impossible');
                Stop;
    Step 5  for j = i + 1, ..., n:
                set u_ij = 1 / l_ii * (a_ij - sum(k=1, i-1, l_ik * u_kj));    // ith row of U
                set l_ji = 1 / u_ii * (a_ji - sum(k=1, i-1, l_jk * u_ki));    // ith column of L
Step 6  Select l_nn and u_nn satisfying l_nn * u_nn = a_nn - sum(k=1, n-1, l_nk * u_kn);
// if l_nn * u_nn = 0, then A = LU but A is singular
Step 7  Output(l_ij for j = 1, ..., i and i = 1, ..., n);
        Output(u_ij for j = i, ..., n and i = 1, ..., n);
        STOP;

Special Types of Matrices⚓︎

Strictly Diagonally Dominant Matrix⚓︎

严格对角占优矩阵(strictly diagonally dominant matrix) 满足：

即对角线元素的绝对值是其所在行内所有元素的绝对值中的最大者。

Choleski's Method for Positive Definite Matrix⚓︎

定理

正定矩阵 \(A\) 的性质：

a. \(A\) 是非奇异的

b. \(a_{ii} > 0, i = 1, 2, \dots n\)

c. \(\max_{1 \le k, j \le n} |a_{kj}| \le \max_{1 \le i \le n} |a_{ii}|\)

d. \((a_{ij})^2 < a_{ii} a_{jj}, i \ne j\)

PPT 上还有这些性质：

• \(A^{-1}\) 也是正定的
• \(A\) 的每个前导主子矩阵(leading principal submatrices) \(A_k\) 的行列式 (determinant) 都是正的

前导主子矩阵

矩阵 \(A\) 的前导主子矩阵为 \(A_k = \begin{bmatrix}a_{11} & a_{12} & \dots & a_{1k} \\ a_{21} & a_{22} & \dots & a_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ a_{k1} & a_{k2} & \dots & a_{kk} \end{bmatrix}, 1 \le k \le n\)

证明

我们将 \(A = LU\) 中的 \(U\) 进一步拆分成对角矩阵 \(D\) 和单位上三角矩阵 \(\widetilde{U}\)：

可以推导出：\(A\) 是对称矩阵（\(A = A^T \rightarrow LU = LD \widetilde{U} = \widetilde{U^T} DL^T\)）\(\Rightarrow L = \widetilde{U}^T \Rightarrow A = LDL^T\)。这样我们得到了另一种矩阵分解—— \(LDL^T\) 分解：

Step 1  for i = 1, ..., n do Steps 2-4:
    Step 2  for j = 1, ..., i - 1, set v_j = l_ij * d_j;
    Step 3  Set d_i = a_ii - sum(j=1, i-1, l_ij * v_j);
    Step 4  for j = i + 1, ..., n, set l_ji = (a_ji - sum(k=1, i-1, l_jk * v_k)) / d_i;

令 \(D^{\frac{1}{2}} = \begin{bmatrix}\sqrt{u_{11}} & & & \\ & \sqrt{u_{22}} & & \\ & & & & \\ & & & \sqrt{u_{nn}}\end{bmatrix}\)，\(\widetilde{L} = LD^{\frac{1}{2}}\) 仍然是一个上三角矩阵，因此 \(A = \widetilde{L} \widetilde{L}^T\)

因为前导主子矩阵都是正的，所以 \(u_{ii} > 0\)

综上，若 \(A\) 是正定矩阵，那么：

• 当 \(L\) 是一个对角线元素均为 1 的下三角矩阵，并且 \(D\) 是一个对角项均为正数的对角矩阵时，\(A\) 可被分解为 \(LDL^T\)
• 当 \(L\) 是一个对角线上均为非零元素的下三角矩阵时，\(A\) 可被分解为 \(LL^T\)

Step 1  set l_11 = sqrt(a_11);
Step 2  for j = 2, ..., n, set l_j1 = a_j1 / l_11;
Step 3  for i = 2, ..., n - 1 do steps 4 and 5
    Step 4  set l_ii = sqrt(a_ii - sum(pow(l_ik, 2), 1, i - 1))
    // LDL^T is faster, but must be modified to solve Ax = b

    Step 5  for j = i + 1, ..., n, set l_ji = (a_ji - sum(l_jk * l_ik, 1, i - 1)) / l_ii;
Step 6  set l_nn = sqrt(a_nn - sum(pow(l_nk, 2), 1, n - 1))
Step 7  output (l_ij for j = 1, ..., i and i = 1, ..., n);
Stop.

Crout Reduction for Tridiagonal Linear System⚓︎

当 \(p = q = 2\) 时，\(w = 3\)，此时的矩阵称为三对角矩阵(tridiagonal matrix)，其形式如下：

对于上述形式的线性方程组（\(A \bm{x} = \bm{f}\)），我们采用一种特殊的 \(LU\) 分解，称为 Crout 分解。下面给出具体求解步骤：

1. 寻找矩阵 \(A\) 的 Crout 分解，\(L, U\) 分别为：

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

2. 求解 \(L\bm{y} = \bm{f} \Rightarrow y_1 = \dfrac{f_1}{l_{11}}, y_i = \dfrac{(f_i - l_{i, i-1} y_{i-1})}{l_{ii}} (i = 2, \dots, n)\)

3. 求解 \(U\bm{x} = \bm{y} \Rightarrow x_n = y_n, x_i = y_i - u_{i,i+1} x_{i+1}\)

Step 1  Set l_11 = a_11;
            u_12 = a_12 / l_11;
            z_1 = a_{1,n+1} / l_11;
Step 2  for i = 2, ..., n-1 set l_{i, i-1} = a_{i, i-1};    // ith row of L
                                l_ii = a_ii - l_{i, i-1} * u_{i-1, i};
                                u_{i, i+1} = a_{i, i+1} / l_ii;    // (i+1)th column of U
                                z_i = (a_{i, n+1} - l_{i, i-1} * z_{i-1}) / l_ii;
Step 3  Set l_{n, n-1} = a_{n, n-1};    // nth row of L
            l_nn = a_nn - l_{n, n-1} * u_{n-1, n};
            z_n = (a_{n, n+1} - l_{n, n-1} * z_{n-1}) / l_nn;
// Step 4 and 5 solve Ux = z
Step 4  x_n = z_n;
Step 5  for i = n-1, ..., 1 set x_i = z_i - u_{i, i+1} * x_{i+1};
Step 6  Output(x_1, ..., x_n);
        STOP;

对应的作业练习📝

对应小测 3 和 4💯