Chap 7: Iterative Techniques in Matrix Algebra⚓︎

目标：同第 6 章，还是求解 \(A \bm{x} = \bm{b}\)。

上述思路的优势在于：

• 可以通过迭代次数来控制精度
• 迭代技术被实际运用于求解稀疏的(sparse) 线性方程组

接下来我们需要分析：

• 如何设计一个迭代方案
• 在何种条件下序列将会收敛
• 某个方法的收敛速度有多快
• 误差评估

Norms of Vectors and Matrices⚓︎

Vector Norms⚓︎

在 \(R^n\) 上的向量范数(vector norm) 是一个 \(R^n \rightarrow R\) 的函数 \(\| \cdot \|\)。\(\forall \bm{x}, \bm{y} \in R^n, \alpha \in C\)，它满足以下性质：

• 正定：\(\| \bm{x} \| \ge 0;\ \| \bm{x} \| = 0 \Leftrightarrow \bm{x} = \bm{0}\)

  此“正定”非彼正定

• 齐次(homogeneous)：\(\| \alpha \bm{x} \| = |\alpha| \cdot \| \bm{x} \|\)

• 三角不等式：\(\| \bm{x} + \bm{y}\| \le \| \bm{x} \| + \| \bm{y} \|\)

一些常用的范数：

• \(\| \bm{x} \|_1 = \sum\limits_{i=1}^n |x_i|\)
• \(\| \bm{x} \|_2 = \sqrt{\sum\limits_{i=1}^n |x_i|^2}\)（欧几里得范数）
• \(\| \bm{x} \|_p = \Big( \sum\limits_{i=1}^n |x_i|^p \Big)^{\frac{1}{p}}\)（\(p\) 范数）
• \(\| \bm{x} \|_\infty = \max\limits_{1 \le i \le n} |x_i|\)（无穷范数）

注：\(\lim\limits_{p \rightarrow \infty} \| \bm{x} \|_p = \| \bm{x} \|_{\infty}\)

一些定义和定理

向量的收敛性：

定义 1 定理 1

若 \(\forall \varepsilon > 0\)，\(\exists N(\varepsilon) \in N\)，使得 \(\|\bm{x}^{(k)} - \bm{x}\| < \varepsilon\) 成立，那么在 \(R^n\) 上的向量序列 \(\{ \bm{x}^{(k)} \}_{k=1}^{\infty}\) 关于范数 \(\| \cdot \|\) 收敛到 \(\bm{x}\)。

当且仅当 \(\lim\limits_{k \rightarrow \infty} x_i^{(k)} = x_i\ (i = 1, 2, \dots, n)\) 时，那么在 \(R^n\) 上的向量序列 \(\{\bm{x}\}_{k=1}^{\infty}\) 关于范数 \(\| \cdot \|_{\infty}\) 收敛到 \(\bm{x}\)。

范数的等价性：

定义 2 定理 2

若存在正常数 \(C_1, C_2\)，使得 \(C_1 \|\bm{x}\|_B \le \|\bm{x}\|_A \le C_2 \| \bm{x} \|_B\)，那么 \(\| \cdot \|_A\) 和 \(\| \cdot \|_B\) 是等价的。

所有在 \(R^n\) 上的向量范数都是等价的。

下面证明了 \(\| \cdot \|_2\) 和 \(\| \cdot \|_\infty\) 的等价性：

证明

Matrix Norms⚓︎

对于所有规模为 \(n \times n\) 的矩阵，矩阵范数(matrix norms) 是一个实数值函数 \(\| \cdot \|\)。对于所有规模为 \(n \times n\) 的矩阵 \(A, B\) 以及所有的 \(\alpha \in C\)，满足：

• 正定：\(\| A \| \ge 0;\ \| A \| = 0 \Leftrightarrow A = O\)
• 齐次：\(\| \alpha A \| = |\alpha| \cdot \| A \|\)
• 三角不等式：\(\| A + B \| \le \| A \| + \| B \|\)
• 一致性(consistency)：\(\| AB \| \le \| A \| \cdot \| B \|\)

下面整理了一些常用的范数：

• 弗罗贝尼乌斯范数(Frobenius norm)：\(\| A \|_F = \sqrt{\sum\limits_{i=1}^n \sum\limits_{j=1}^n |a_{ij}|^2}\)

• 自然范数(natural norm)

  • 算子范数(operator norm)（和向量范数 \(\| \cdot \|\) 关联，所以也可称为 \(p\) 范数）

    \[ \| A \|_p = \max\limits_{\bm{x} \ne \bm{0}} \dfrac{\| A \bm{x} \|_p}{\| \bm{x} \|_p} = \max\limits_{\| \bm{x} \|_p = 1} \| A\bm{x} \|_p \]

  • 无穷范数：\(\| A \|_{\infty} = \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}|\)

    证明 \(\|A \|_{\infty} = \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}|\)

    1. 证明 \(\| A \|_{\infty} = \max\limits_{\| \bm{x} \|_{\infty} = 1} \| A \bm{x} \|_{\infty} \le \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}|\)

       \[ \| A \bm{x} \|_{\infty} = \max\limits_{1 \le i \le n} |(A \bm{x})_i| = \max\limits_{1 \le i \le n} |\sum\limits_{j=1}^n a_{ij} x_j| \le \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}| \cdot \max\limits_{1 \le i \le n} |x_j| \]

    2. 证明 \(\| A \|_{\infty} = \max\limits_{\| \bm{x} \|_{\infty} = 1} \| A\bm{x} \|_{\infty} \ge \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}|\)

       • 令第 \(p\) 行为最大行，即满足 \(\sum\limits_{j=1}^n |a_{pj}| = \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}|\)
       • 取一个特殊的单位向量 \(\bm{x}\) 使得 \(x_j = \begin{cases} 1, & \text{if } a_{pj} \ge 0 \\ -1, & \text{if } a_{pj} < 0 \end{cases}\)
       \[ \| A \bm{x} \|_{\infty} = \max\limits_{1 \le i \le n} \Big| \sum\limits_{j=1}^n a_{ij} x_j\Big| \ge \Big| \sum\limits_{j=1}^n a_{pj} x_j \Big| = \Big| \sum\limits_{j=1}^n |a_{pj}| \Big| = \max\limits_{1 \le i \le n} \sum\limits_{j=1}^n |a_{ij}| \]

  • 1- 范数：\(\| A \|_1 = \max\limits_{1 \le j \le n} \sum\limits_{i=1}^n |a_{ij}|\)

  • 2- 范数（又称谱范数(spectral norm)）：\(\| A \|_2 = \sqrt{\lambda_{\max} (A^T A)}\)

Eigenvalues and Eigenvectors⚓︎

矩阵 \(A\) 的谱半径(spectral radius) \(\rho(A) = \max{| \lambda |}\)，其中 \(\lambda\) 是 \(A\) 的特征值（复数范围，所以 \(|\lambda|\) 表示特征值的模长）。

若 \(\forall i, j = 1, 2, \dots, n\)，有 \(\lim\limits_{k \rightarrow \infty} (A^k)_{ij} = 0\)，那么称规模为 \(n \times n\) 的矩阵 \(A\) 是收敛的。

Iterative Techniques for Solving Linear Systems⚓︎

当线性方程组维度不高的时候，很少会用到迭代方法，因为迭代方法达到足够精度所需的时间超过用高斯消元法计算的时间。不过在大型的线性方程组，迭代方法在存储和计算上的效率优势会变得明显。

下面介绍一些常见的求解线性方程组的迭代方法。

Jacobi Iterative Method⚓︎

对于线性方程组 \(\begin{cases}a_{11}x_1 + a_{12}x_2 + \dots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \dots + a_{2n}x_n = b_2 \\ \dots \\ a_{n1}x_1 + a_{n2}x_2 + \dots + a_{nn}x_n = b_n\end{cases}\)，当 \(a_{ii} \ne 0\) 时，可以解得：

用矩阵形式表示上述线形方程组，并转化为以下形式：

那么：

因此该迭代法的递推公式为：\(\bm{x}^{(k)} = T_j \bm{x}^{(k-1)} + \bm{c_j}\)，其中 \(T_j\) 被称为雅可比迭代矩阵(Jacobi iterative matrix)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

Step 1  Set k = 1;
Step 2  while (k <= N_max) do step 3-6
        Step 3  for i = 1, ..., n
                    Set X_i = (b_i - sum(j=1, j!=i, j<=n, a_ij * X0_j)) / a_ii;  // compute x^k
        Step 4  if norm(X - X0)_infty = max(1<=i<=n, X_i - X0_i) < TOL then Output(X[]);
                STOP;    // successful
        Step 5  for i = 1, ..., n  Set X0[] = X[];  // update X0
        Step 6  Set k++;
Step 7  Output (Maximum number of iterations exceeded);
        STOP.    // unsuccessful

Gauss-Seidel Iterative Method⚓︎

观察线性方程组的解：

用矩阵形式表述为：

因此该迭代法的递推公式为：\(\bm{x}^{(k)} = T_g \bm{x}^{(k-1)} + c_g \bm{b}\)，其中 \(T_g\) 为高斯 - 塞德尔迭代矩阵(Gauss-Seidel iterative matrix)。

Step 1  Set k = 1;
Step 2  while (k <= N_max) do step 3-6
        Step 3  for i = 1, ..., n
                    Set X_i = (-sum(j=1, i-1, a_ij * x_j) - sum(j=i+1, n, a_ij * X0_j) + b_i) / a_ii;  // compute x^k
        Step 4  if norm(X - X0)_infty = max(1<=i<=n, X_i - X0_i) < TOL then Output(X[]);
                STOP;    // successful
        Step 5  for i = 1, ..., n  Set X0[] = X[];  // update X0
        Step 6  Set k++;
Step 7  Output (Maximum number of iterations exceeded);
        STOP.    // unsuccessful

Convergence of Iterative Method⚓︎

现在我们来考察迭代法 \(\bm{x}^{(k)} = T\bm{x}^{(k-1)} + \bm{c}\) 的收敛性。

Relaxation Methods⚓︎

接下来我们从残差向量的角度审视高斯 - 塞德尔方法：

令 \(x_i^{(k)} = x_i^{(k-1)} + \textcolor{red}{\omega} \dfrac{r_i^{(k)}}{a_{ii}}\)。对于正数 \(\omega\) 的某种选择，我们能够减少残差向量 \(\bm{r}\) 的范数，并且得到更快的收敛。这样的方法称为松弛法(relaxation methods)。根据 \(\omega\) 的大小，我们有以下几类松弛法：

• \(0 < \omega < 1\)：欠松弛法(under-relaxation methods)
• \(\omega = 1\)：高斯 - 塞德尔方法
• \(\omega > 1\)：逐次超松弛法(successive over-relaxation methods, SOR)

用矩阵形式可以表述为：

也就是说，SOR 迭代法的递推公式为：\(\bm{x}^{(k)} = T_\omega \bm{x}^{(k-1)} + \bm{c}_\omega \bm{b}\)。

显然，矩阵 \(T_\omega\) 的表示过于复杂，因而很难准确得到它的谱半径。但好在我们还可以走以下“捷径”，同时也有助于我们选择合适的 \(\omega\) 值。

讨论

题目解答

给定 \(A = \begin{bmatrix}2 & 1 \\ 1 & 2\end{bmatrix}, \bm{b} = \begin{bmatrix}1 \\ 2\end{bmatrix}\)，以及迭代法 \(\bm{x^{(k)}} = \bm{x^{(k-1)}} + \omega (A\bm{x^{(k-1)}} - \bm{b})\)，那么：

• 当 \(\omega\) 取什么值时，该方法会收敛？
• 当 \(\omega\) 取什么值时，该方法的收敛速度最快？

考虑 \(T = I + \omega A\) 的特征值，解得 \(\lambda_1 = 1 + \omega, \lambda_2 = 1 + 3 \omega\)

• 要使该方法收敛，需满足 \(\rho(T) < 1 \quad \Rightarrow -\dfrac{2}{3} < \omega < 0\)
• 现在考虑 \(\rho(T) = \max\{|1 + \omega|, |1 + 3\omega|\}\) 在何时取最小值（即收敛速度最快的时候）——我们可以画图研究：

发现当 \(\omega = -\dfrac{1}{2}\) 时取值最小。

Step 1  Set k = 1;
Step 2  while (k <= N_max) do step 3-6
        Step 3  for i = 1, ..., n
                    Set X_i = (1 - omega) * X0_i + (omega * (-sum(j=1, i-1, a_ij * x_j) - sum(j=i+1, n, a_ij * X0_j) + b_i)) / a_ii  // compute x^k
        Step 4  if norm(X - X0)_infty = max(1<=i<=n, X_i - X0_i) < TOL then Output(X[]);
                STOP;    // successful
        Step 5  for i = 1, ..., n  Set X0[] = X[];  // update X0
        Step 6  Set k++;
Step 7  Output (Maximum number of iterations exceeded);
        STOP.    // unsuccessful

Error Bounds and Iterative Refinement⚓︎

目标：研究 \(A\) 和 \(\bm{b}\) 的误差会如何影响 \(A \bm{x} = \bm{b}\) 的解 \(\bm{x}\)？

假设 \(A\) 是准确的，\(\bm{b}\) 有误差 \(\delta \bm{b}\)，那么带有误差的解可以写成 \(\bm{x} + \delta \bm{x}\)。可以得到：

其中 \(\| A \| \cdot \| A^{-1} \|\) 被称为相对放大因子(relative amplification factor)。

假如 \(\bm{b}\) 是准确的，\(A\) 有误差 \(\delta A\)，那么带有误差的解可以写成 \(\bm{x} + \delta \bm{x}\)。可以得到：

其中 \(\| A \| \cdot \| A^{-1} \|\) 是误差放大的关键因子，被称为条件数(condition number) \(K(A)\)。

• 如果 \(K(A)\) 接近 1，那么矩阵 \(A\) 是良态的(well-conditioned)
• 如果 \(K(A)\) 比 1 大很多，那么矩阵 \(A\) 是病态的(ill-conditioned)

讨论

问题解答

对于给定的 \(A = \begin{bmatrix}1 & 0.99 \\ 0.99 & 0.98\end{bmatrix}, \bm{b} = \begin{bmatrix}1.99 \\ 1.97\end{bmatrix}\)

1. 计算 \(K(A)_2\)
2. 给 \(\bm{b}\) 一个很小的扰动 \(\delta \bm{b} = \begin{bmatrix}-0.97 \times 10^{-4} \\ 0.106 \times 10^{-3}\end{bmatrix}\)，原方程的解会有什么变化？

直接从前辈笔记上 copy 下来的：

Iterative Refinement⚓︎

迭代优化(iterative refinement) 的步骤为：

1. \(A \bm{x} = \bm{b} \Rightarrow\) 近似解 \(\bm{x}_1\)
2. \(\bm{r}_1 = \bm{b} - A \bm{x}_1\)
3. \(A \bm{d}_1 = \bm{r}_1 \Rightarrow \bm{d}_1\)
   • 如果 \(\bm{d}_1\) 是精确的，那么 \(\bm{x}_2 = \bm{x}_1 + A^{-1} (\bm{b} - A\bm{x}_1) = A^{-1} \bm{b}\)，\(\bm{x}_2\) 也是精确的。
4. \(\bm{x}_2 = \bm{x}_1 + \bm{d}_1\)

之后重复 2-4 步。

算法：迭代优化

对应的作业练习📝

对应小测 5 和 6💯

评论区