基于Krylov子空间的GMRES算法解非对称线性方程组
题目描述
考虑求解大型稀疏非对称线性方程组 \(A\mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{n \times n}\) 为非奇异矩阵（但不一定对称正定），\(\mathbf{b} \in \mathbb{R}^n\)。由于矩阵规模大且稀疏，直接法（如LU分解）可能内存消耗过高或计算量过大。GMRES（广义最小残差法）是一种基于Krylov子空间的迭代算法，它通过最小化残差范数来逼近解，尤其适用于非对称矩阵。本题目详细讲解GMRES算法的原理、构造过程和计算步骤。
解题过程
-
算法核心思想
GMRES的目标是在第 \(m\) 步迭代时，从仿射空间 \(\mathbf{x}_0 + \mathcal{K}_m(A, \mathbf{r}_0)\) 中选择近似解 \(\mathbf{x}_m\)，使得残差范数 \(\|\mathbf{r}_m\|_2 = \|\mathbf{b} - A\mathbf{x}_m\|_2\) 最小。其中：- \(\mathbf{x}_0\) 为初始猜测解（通常取零向量或估计值）。
- \(\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0, A\mathbf{r}_0, \dots, A^{m-1}\mathbf{r}_0\}\) 是由初始残差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\) 生成的Krylov子空间。
-
Arnoldi过程构造正交基
由于Krylov子空间的基向量 \(\{ \mathbf{r}_0, A\mathbf{r}_0, \dots \}\) 可能逐渐线性相关，需通过Arnoldi过程生成一组标准正交基 \(\{\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m\}\) 来稳定计算：- 初始化：归一化初始残差向量 \(\mathbf{v}_1 = \mathbf{r}_0 / \|\mathbf{r}_0\|_2\)。
- 迭代计算（对于 \(j = 1, 2, \dots, m\)）：
- 计算 \(\mathbf{w}_j = A\mathbf{v}_j\)。
- 对 \(\mathbf{w}_j\) 与所有已生成的基向量 \(\mathbf{v}_1, \dots, \mathbf{v}_j\) 进行正交化（使用Modified Gram-Schmidt过程）：
\[ h_{ij} = \mathbf{v}_i^T \mathbf{w}_j, \quad \mathbf{w}_j \leftarrow \mathbf{w}_j - h_{ij}\mathbf{v}_i \quad (i = 1, \dots, j). \]
3. 令 $ h_{j+1,j} = \|\mathbf{w}_j\|_2 $，若 $ h_{j+1,j} = 0 $ 则提前终止（此时精确解已找到）。
4. 生成新的基向量 $ \mathbf{v}_{j+1} = \mathbf{w}_j / h_{j+1,j} $。
该过程可写为矩阵形式：
\[ A V_m = V_{m+1} \bar{H}_m, \]
其中 \(V_m = [\mathbf{v}_1, \dots, \mathbf{v}_m]\) 的列构成Krylov子空间的标准正交基，\(\bar{H}_m \in \mathbb{R}^{(m+1) \times m}\) 为上Hessenberg矩阵（元素 \(h_{ij}\) 存储正交化系数）。
- 最小化残差问题转化
设近似解 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)，其中 \(\mathbf{y}_m \in \mathbb{R}^m\) 为系数向量。则残差可表示为：
\[ \mathbf{r}_m = \mathbf{b} - A\mathbf{x}_m = \mathbf{r}_0 - A V_m \mathbf{y}_m = \mathbf{r}_0 - V_{m+1} \bar{H}_m \mathbf{y}_m. \]
由于 \(\mathbf{r}_0 = \|\mathbf{r}_0\|_2 \mathbf{v}_1 = \|\mathbf{r}_0\|_2 V_{m+1} \mathbf{e}_1\)（其中 \(\mathbf{e}_1 = [1, 0, \dots, 0]^T \in \mathbb{R}^{m+1}\)），残差范数可简化为：
\[ \|\mathbf{r}_m\|_2 = \| \|\mathbf{r}_0\|_2 \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m \|_2. \]
因此，原问题转化为求解最小二乘问题：
\[ \text{最小化 } \| \beta \mathbf{e}_1 - \bar{H}_m \mathbf{y}_m \|_2, \quad \beta = \|\mathbf{r}_0\|_2. \]
- 利用QR分解求解最小二乘问题
对上Hessenberg矩阵 \(\bar{H}_m\) 进行QR分解（通过Givens旋转高效更新）：- 初始化 \(Q_0 = I_{m+1}\)，\(R_0 = \bar{H}_m\)（实际上逐步构造）。
- 对每列 \(j = 1, \dots, m\)，设计Givens旋转矩阵 \(G_j\) 消去 \(\bar{H}_m\) 的第 \(j+1\) 行元素，累积得到正交矩阵 \(Q_m = G_m \cdots G_1\) 和上三角矩阵 \(R_m = Q_m \bar{H}_m \in \mathbb{R}^{(m+1) \times m}\)（最后一行全零）。
- 记 \(\mathbf{g}_m = Q_m (\beta \mathbf{e}_1) = [\gamma_1, \dots, \gamma_{m+1}]^T\)，则最小二乘问题变为：
\[ \text{最小化 } \| \mathbf{g}_m - R_m \mathbf{y}_m \|_2. \]
- 忽略最后一行（零残差），求解上三角系统 \(R_m^{(1:m, 1:m)} \mathbf{y}_m = \mathbf{g}_m^{(1:m)}\) 得到 \(\mathbf{y}_m\)。
-
算法流程与重新启动策略
- 输入：矩阵 \(A\)，向量 \(\mathbf{b}\)，初始解 \(\mathbf{x}_0\)，最大子空间维数 \(m\)，容差 \(\epsilon\)。
- 步骤：
- 计算初始残差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)，\(\beta = \|\mathbf{r}_0\|_2\)，若 \(\beta < \epsilon\) 则输出 \(\mathbf{x}_0\)。
- 令 \(\mathbf{v}_1 = \mathbf{r}_0 / \beta\)。
- 进行 \(m\) 步Arnoldi过程，构建 \(V_m\) 和 \(\bar{H}_m\)。
- 对 \(\bar{H}_m\) 应用Givens旋转，更新QR分解并计算残差范数 \(\|\mathbf{r}_m\|_2 = |\gamma_{m+1}|\)。
- 若 \(\|\mathbf{r}_m\|_2 < \epsilon\) 或子空间维度达到 \(m\)，则求解最小二乘问题得到 \(\mathbf{y}_m\)，计算近似解 \(\mathbf{x}_m = \mathbf{x}_0 + V_m \mathbf{y}_m\)。
- 若不收敛且未达到最大迭代次数，设 \(\mathbf{x}_0 \leftarrow \mathbf{x}_m\) 并重新启动（从步骤1开始）。
-
关键细节与注意事项
- 存储成本：需存储 \(m+1\) 个长度为 \(n\) 的基向量，当 \(m\) 较大时内存可能受限。因此常采用重新启动策略（如GMRES(\(m\))），每 \(m\) 步重启一次以控制内存。
- 预处理技术：为提高收敛速度，通常使用预处理矩阵 \(M\) 将原系统转化为 \(M^{-1} A \mathbf{x} = M^{-1} \mathbf{b}\)，使矩阵更接近单位矩阵（例如不完全LU分解预处理）。
- 收敛性：对于正定矩阵（即 \(A + A^T\) 正定），GMRES保证收敛；对于一般矩阵，收敛理论较复杂，可能停滞（如遇到不变子空间）。
总结
GMRES算法通过Arnoldi过程构建Krylov子空间的正交基，将最小化残差问题转化为小规模最小二乘问题，利用QR分解高效求解。其优势在于适用于任意非奇异矩阵，但需注意内存管理和预处理以加速收敛。该算法是求解大型稀疏非对称线性方程组的重要工具。