基于Krylov子空间的TFQMR算法解非对称线性方程组
1. 问题描述
在科学计算中，我们经常需要求解大规模非对称线性方程组：
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个非对称的大型稀疏矩阵，\(\mathbf{b} \in \mathbb{R}^n\) 是已知向量，\(\mathbf{x} \in \mathbb{R}^n\) 是待求解向量。当 \(A\) 规模很大时，直接解法（如LU分解）往往因内存和计算成本过高而不可行。这时，基于Krylov子空间的迭代法成为首选。
TFQMR（Transpose-Free Quasi-Minimal Residual）算法是此类方法中的一种，它不需要显式计算矩阵 \(A\) 的转置与向量的乘积，适用于 \(A\) 难以显式转置或 \(A^T\) 不易计算的情况。TFQMR算法能高效求解非对称线性方程组，并在迭代过程中保持较平稳的收敛性。
2. 算法背景与动机
TFQMR源于双共轭梯度法（BiCG）和最小残差法（QMR）的思想：
- BiCG算法：利用两个Krylov子空间 \(\mathcal{K}_m(A, \mathbf{r}_0)\) 和 \(\mathcal{K}_m(A^T, \tilde{\mathbf{r}}_0)\) 构建迭代，但需要计算 \(A^T\mathbf{v}\)，这在某些应用中不便实现。
- QMR算法：通过将BiCG的迭代向量组合，使残差范数在某种意义下最小化，但仍需 \(A^T\)。
- TFQMR：通过巧妙的数学变换，避免了显式计算 \(A^T\)，同时保持了类似QMR的“拟最小残差”性质，提高了数值稳定性和收敛效率。
3. 算法推导与步骤
TFQMR的核心思想是利用BiCG的迭代参数，但通过重新组合迭代向量来构造一组新的近似解序列，使残差范数尽可能小。
步骤1：初始化
给定初始猜测解 \(\mathbf{x}_0\)，计算初始残差：
\[\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0 \]
选择任意向量 \(\tilde{\mathbf{r}}_0\) 满足 \(\tilde{\mathbf{r}}_0^T \mathbf{r}_0 \neq 0\)（通常取 \(\tilde{\mathbf{r}}_0 = \mathbf{r}_0\)）。
设定初始辅助变量：
\[\mathbf{w}_0 = \mathbf{u}_0 = \mathbf{r}_0, \quad \mathbf{d}_0 = 0, \quad \tau_0 = \|\mathbf{r}_0\|_2, \quad \theta_0 = 0, \quad \eta_0 = 0 \]
初始化近似解 \(\mathbf{x}_0\)。
步骤2：迭代过程（对于 \(m = 0, 1, 2, \dots\)）
TFQMR的迭代分为两个循环：
- 外层循环（索引 \(m\)）：每两次迭代产生一个近似解。
- 内层循环（索引 \(j = 1, 2\)）：每次计算BiCG类型的向量，并更新解。
第 \(m\) 次外层迭代的具体计算：
- BiCG步骤（计算 \(\alpha_m, \rho_m, \beta_m\)）：
\[ \rho_m = \tilde{\mathbf{r}}_0^T \mathbf{r}_m, \quad \beta_m = (\rho_m / \rho_{m-1}) \cdot (\alpha_{m-1} / \omega_{m-1}) \]
其中 \(\omega_{m-1}\) 是前一次迭代的辅助标量（初始 \(\omega_{-1}=1\)）。
- 如果 \(\rho_m = 0\)，算法失败（出现“中断”）。
- 更新方向向量 \(\mathbf{u}_m\)：
\[ \mathbf{u}_m = \mathbf{r}_m + \beta_m (\mathbf{u}_{m-1} - \omega_{m-1} A\mathbf{u}_{m-1}) \]
- 矩阵向量乘：
\[ \mathbf{v}_m = A\mathbf{u}_m \]
- 计算 \(\alpha_m\)：
\[ \alpha_m = \rho_m / (\tilde{\mathbf{r}}_0^T \mathbf{v}_m) \]
- 计算中间残差 \(\mathbf{r}_{m+1/2}\)：
\[ \mathbf{r}_{m+1/2} = \mathbf{r}_m - \alpha_m \mathbf{v}_m \]
- 矩阵向量乘：
\[ \mathbf{w}_{m+1/2} = A \mathbf{r}_{m+1/2} \]
- 计算 \(\omega_m\)：
\[ \omega_m = (\mathbf{w}_{m+1/2}^T \mathbf{r}_{m+1/2}) / (\mathbf{w}_{m+1/2}^T \mathbf{w}_{m+1/2}) \]
如果分母为0，则取 \(\omega_m = 0\)。
- 更新残差 \(\mathbf{r}_{m+1}\)：
\[ \mathbf{r}_{m+1} = \mathbf{r}_{m+1/2} - \omega_m \mathbf{w}_{m+1/2} \]
- TFQMR特有关键步骤：构造新的解向量序列
定义辅助向量：
\[ \mathbf{d}_{2m} = \mathbf{u}_m + \theta_{2m}^2 (\omega_m / \alpha_m) \mathbf{d}_{2m-1} \]
\[ \mathbf{d}_{2m+1} = \mathbf{r}_{m+1/2} + \theta_{2m+1}^2 (\omega_m / \alpha_m) \mathbf{d}_{2m} \]
其中标量 \(\theta_k\) 和 \(\eta_k\) 通过以下递推更新：
- 对于 \(j = 1, 2\)：
\[ \tau_{2m+j} = \tau_{2m+j-1} \cdot \| \mathbf{v}_{m,j} \|_2 / \sqrt{1 + \theta_{2m+j-1}^2} \]
\[ \theta_{2m+j} = \tau_{2m+j} / (\eta_{2m+j-1} \sqrt{1 + \theta_{2m+j-1}^2}) \]
\[ \eta_{2m+j} = \eta_{2m+j-1} \cdot \theta_{2m+j-1} / \theta_{2m+j} \]
（这里的 $\mathbf{v}_{m,1} = \mathbf{v}_m$, $\mathbf{v}_{m,2} = \mathbf{w}_{m+1/2}$）
- 更新近似解：
每完成一次内层循环（即每计算一个 \(\mathbf{d}_k\)），就更新解：
\[ \mathbf{x}_k = \mathbf{x}_{k-1} + \eta_k \mathbf{d}_k \]
其中 $k = 2m, 2m+1$。
步骤3：收敛检查
检查残差范数 \(\|\mathbf{r}_{m+1}\|_2\) 是否小于预设容差 \(\epsilon\)。若满足则停止，输出当前近似解 \(\mathbf{x}_k\)；否则继续迭代。
4. 算法特点与优势
- 转置无关性：不需要计算 \(A^T\) 与向量的乘积，适用于 \(A^T\) 不可用或计算代价高的情况。
- 拟最小残差性质：通过组合BiCG的迭代向量，使残差范数在某种加权意义下最小，收敛曲线通常比BiCG更平滑。
- 计算开销：每步迭代需要 两次矩阵向量乘法（\(A\mathbf{u}_m\) 和 \(A\mathbf{r}_{m+1/2}\)），以及一些向量内积和标量运算，总体与BiCGSTAB相当。
- 数值稳定性：通过递推更新 \(\theta_k, \eta_k\) 避免了直接计算残差范数可能出现的数值不稳定问题。
5. 实际应用中的注意事项
- 预处理：为了加速收敛，通常使用预处理技术，即求解等价系统 \(M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}\)，其中 \(M\) 是近似于 \(A\) 且易于求逆的矩阵（如不完全LU分解）。
- 收敛性：TFQMR可能在某些问题中出现停滞或振荡，这时可结合重启策略或切换其他算法（如GMRES）。
- 实现细节：在实际代码中，需注意向量存储和更新顺序，避免不必要的重复计算。
6. 简单示例
考虑一个小规模非对称矩阵：
\[A = \begin{pmatrix} 4 & -1 & 0 \\ -2 & 5 & -1 \\ 0 & -1 & 3 \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \]
取初始猜测 \(\mathbf{x}_0 = \mathbf{0}\)，容差 \(\epsilon = 10^{-6}\)。
手动执行TFQMR迭代（略去具体数值），通常会在10步以内收敛到真解 \(\mathbf{x} \approx (0.4, 0.8, 1.2)^T\)。
7. 总结
TFQMR算法是求解大规模非对称线性方程组的一种高效Krylov子空间方法。它通过避免矩阵转置运算、结合拟最小残差策略，在保持较低计算成本的同时，提供了较为平滑的收敛行为。该算法在计算流体力学、电磁场仿真等领域的稀疏线性系统求解中有着广泛应用。