基于Krylov子空间的QMR算法解非对称线性方程组
算法描述
QMR（Quasi-Minimal Residual，拟极小残差）算法是一种用于求解大规模、稀疏、非对称线性方程组 \(Ax = b\) 的迭代方法。它属于Krylov子空间投影法家族，旨在克服GMRES算法存储需求随迭代次数线性增长，以及BiCG算法可能出现的迭代过程不稳定（如除零）的问题。QMR算法的核心思想是在非对称Lanczos双正交化过程生成的Krylov子空间中，寻找一个使得残量向量的“拟极小”的近似解，其残量范数曲线通常比BiCG更平滑。该算法特别适用于大规模稀疏非对称矩阵，且通常与适当的预处理技术结合以加速收敛。
解题过程详解
1. 问题形式与Krylov子空间
我们要求解线性方程组：
\[Ax = b \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是大规模稀疏非对称矩阵，\(b \in \mathbb{R}^{n}\) 是已知右端向量，\(x \in \mathbb{R}^{n}\) 是未知解向量。
给定初始猜测解 \(x_0\)，初始残差 \(r_0 = b - A x_0\)。目标是构建近似解序列 \(\{x_k\}\)，使得 \(x_k \in x_0 + \mathcal{K}_k(A, r_0)\)，其中 Krylov 子空间为：
\[\mathcal{K}_k(A, r_0) = \text{span}\{r_0, A r_0, A^2 r_0, \dots, A^{k-1} r_0\}. \]
2. 非对称Lanczos双正交化过程
QMR算法的基础是非对称Lanczos过程，它生成两组双正交的向量序列 \(\{v_1, v_2, \dots, v_k\}\) 和 \(\{w_1, w_2, \dots, w_k\}\)，满足：
\[W_k^T V_k = I_k, \quad V_k = [v_1, \dots, v_k], \quad W_k = [w_1, \dots, w_k]. \]
并且它们满足三项递推关系：
\[\begin{aligned} A V_k &= V_k T_k + \beta_{k+1} v_{k+1} e_k^T, \\ A^T W_k &= W_k T_k^T + \gamma_{k+1} w_{k+1} e_k^T, \end{aligned} \]
其中 \(T_k\) 是一个三对角矩阵（由于非对称性，实际上是“上Hessenberg”形式，但Lanczos过程强制其为三对角）：
\[T_k = \begin{bmatrix} \alpha_1 & \gamma_2 \\ \beta_2 & \alpha_2 & \gamma_3 \\ & \beta_3 & \alpha_3 & \ddots \\ & & \ddots & \ddots & \gamma_k \\ & & & \beta_k & \alpha_k \end{bmatrix}. \]
在理想情况下，\(\gamma_j = \beta_j\)，但在实际浮点运算中，数值误差可能导致双正交性丢失，因此需要适当处理。
具体步骤：
- 初始化：选择初始向量 \(v_1 = r_0 / \|r_0\|_2\)，选择 \(w_1\) 使得 \(w_1^T v_1 = 1\)（通常取 \(w_1 = v_1\) 或随机向量再正交化）。
- 对 \(j = 1, 2, \dots, k\) 进行迭代：
- 计算 \(\alpha_j = w_j^T A v_j\)。
- 更新向量 \(\hat{v}_{j+1} = A v_j - \alpha_j v_j - \beta_j v_{j-1}\)（其中 \(v_0 = 0, \beta_1 v_0 = 0\)）。
- 更新向量 \(\hat{w}_{j+1} = A^T w_j - \alpha_j w_j - \gamma_j w_{j-1}\)（其中 \(w_0 = 0, \gamma_1 w_0 = 0\)）。
- 计算 \(\beta_{j+1} = \|\hat{v}_{j+1}\|_2\)，\(\gamma_{j+1} = w_j^T A v_{j+1} / \alpha_j\) 或通过双正交化条件确定。
- 令 \(v_{j+1} = \hat{v}_{j+1} / \beta_{j+1}\)，\(w_{j+1} = \hat{w}_{j+1} / \gamma_{j+1}\)（需保证 \(w_{j+1}^T v_{j+1} = 1\)，可能需要一次额外的重新正交化）。
3. QMR近似解的构造
在每一步 \(k\)，我们寻找近似解 \(x_k = x_0 + V_k y_k\)，其中 \(y_k \in \mathbb{R}^k\) 是待定系数向量。
对应的残差向量为：
\[r_k = b - A x_k = r_0 - A V_k y_k. \]
由 Lanczos 关系 \(A V_k = V_{k+1} \underline{T}_k\)，其中 \(\underline{T}_k \in \mathbb{R}^{(k+1) \times k}\) 是将 \(T_k\) 下方添加一行 \([0, \dots, 0, \beta_{k+1}]\) 得到的矩阵。
又因为 \(r_0 = \|r_0\|_2 v_1 = \|r_0\|_2 V_{k+1} e_1\)（其中 \(e_1 = [1, 0, \dots, 0]^T \in \mathbb{R}^{k+1}\)），所以：
\[r_k = V_{k+1} \left( \|r_0\|_2 e_1 - \underline{T}_k y_k \right). \]
由于 \(V_{k+1}\) 的列是正交的（在精确算术意义下），我们有：
\[\|r_k\|_2 = \left\| \|r_0\|_2 e_1 - \underline{T}_k y_k \right\|_2. \]
拟极小残差 的思想是：不直接最小化 \(\|r_k\|_2\)（那是GMRES的做法，但这里 \(V_{k+1}\) 的列不一定正交，且 \(T_k\) 不是上三角阵，无法简单用最小二乘），而是通过求解一个“拟极小”问题。QMR算法实际上求解：
\[y_k = \arg\min_{y \in \mathbb{R}^k} \left\| \|r_0\|_2 e_1 - \underline{T}_k y \right\|_2, \]
即忽略 \(V_{k+1}\) 可能非正交的影响，只最小化系数向量的2-范数。这样，原问题转化为一个规模为 \((k+1) \times k\) 的最小二乘问题。
4. 通过LQ分解高效求解
为了高效求解该最小二乘问题，对 \(\underline{T}_k\) 进行 LQ分解（即QR分解的转置形式）：存在正交矩阵 \(Q_k \in \mathbb{R}^{(k+1) \times (k+1)}\) 和上三角矩阵 \(L_k \in \mathbb{R}^{(k+1) \times k}\)（实际上是下三角矩阵的转置，但这里习惯称为LQ），使得：
\[\underline{T}_k = L_k Q_k. \]
但更常见的实现方式是：在迭代过程中，利用Givens旋转逐次将 \(\underline{T}_k\) 化为上三角矩阵 \(R_k\)（即QR分解），从而更新最小二乘解。具体地：
- 在第 \(k\) 步，新增一列到 \(\underline{T}_k\)，得到 \(\underline{T}_{k+1}\) 的前 \(k+1\) 行已三角化，新增的最后一行需要被旋转消元。
- 应用Givens旋转 \(J_1, J_2, \dots, J_k\) 将新行化为上三角形式，同时更新右端项 \(\|r_0\|_2 e_1\)。
- 最终得到上三角系统 \(R_k y_k = g_k\)（其中 \(g_k\) 是旋转后的右端项前 \(k\) 个分量），通过回代即可得到 \(y_k\)。
5. 算法迭代步骤
-
初始化：
给定 \(x_0\)，计算 \(r_0 = b - A x_0\)。
令 \(v_1 = r_0 / \|r_0\|_2\)，选择 \(w_1 = v_1\)（或任意满足 \(w_1^T v_1 = 1\) 的向量）。
设 \(\xi = \|r_0\|_2\)，\(g = \xi e_1\)（用于记录最小二乘右端项）。 -
对 \(k = 1, 2, \dots\) 直到收敛，执行：
- Lanczos双正交化步骤：
- 计算 \(\alpha_k = w_k^T A v_k\)。
- 更新 \(\hat{v} = A v_k - \alpha_k v_k - \beta_k v_{k-1}\)（若 \(k=1\)，则 \(\beta_1 v_0 = 0\)）。
- 更新 \(\hat{w} = A^T w_k - \alpha_k w_k - \gamma_k w_{k-1}\)（若 \(k=1\)，则 \(\gamma_1 w_0 = 0\)）。
- 计算 \(\beta_{k+1} = \|\hat{v}\|_2\)，\(\gamma_{k+1} = \hat{w}^T \hat{v} / \beta_{k+1}\)（确保双正交性）。
- 令 \(v_{k+1} = \hat{v} / \beta_{k+1}\)，\(w_{k+1} = \hat{w} / \gamma_{k+1}\)。
- 扩展矩阵 \(\underline{T}_k\)：
将新的一列 \([\gamma_k, \alpha_k, \beta_{k+1}]^T\) 添加到 \(\underline{T}_k\) 的右侧（注意下标：第 \(k\) 列是 \([\gamma_k, \alpha_k, \beta_{k+1}]^T\)，但通常存储为三对角形式，这里为最小二乘需构造稠密矩阵）。 - 更新QR分解：
对上一步得到的增广矩阵（大小为 \((k+1) \times k\)）应用Givens旋转，使其保持上三角形式 \(R_k\)，并相应更新右端向量 \(g\)。 - 求解最小二乘：
解 \(R_k y_k = g(1:k)\) 得 \(y_k\)。 - 更新近似解：
\(x_k = x_0 + V_k y_k\)（实际计算中，不必存储全部 \(V_k\)，可利用递推关系更新解向量）。 - 检查收敛：
计算残差范数估计 \(\|r_k\|_2 \approx |g(k+1)|\)（即旋转后右端项的第 \(k+1\) 个分量）。若小于给定容差，则终止迭代。
- Lanczos双正交化步骤：
6. 算法特点与注意事项
- 存储需求：QMR需要存储两组 Lanczos 向量 \(v_j, w_j\) 以及递推系数，存储量固定（不随 \(k\) 增长），适合大规模问题。
- 收敛性：通常比 BiCG 更平滑，残差范数不会出现剧烈振荡，但收敛速度类似。
- 数值稳定性：双正交性可能因舍入误差丢失，可引入“回溯正交化”或“选择性正交化”改善。
- 预处理：为加速收敛，常使用左预处理 \(M^{-1} A x = M^{-1} b\) 或左右预处理，此时需相应修改 Lanczos 过程。
总结
QMR算法通过非对称Lanczos双正交化构建Krylov子空间，并在该子空间中求解一个拟极小残差问题，避免了GMRES的存储增长和BiCG的不稳定性。其核心步骤包括：Lanczos双正交化、递推更新三对角矩阵、利用Givens旋转进行QR分解以高效求解最小二乘子问题，从而迭代更新近似解。该算法是大规模稀疏非对称线性方程组求解的重要工具。