基于Krylov子空间的QMR算法解非对称线性方程组

我将为你讲解基于Krylov子空间的QMR（拟最小残差）算法，这是一种用于求解非对称线性方程组的迭代方法。QMR算法通过结合Lanczos双正交化过程和最小残差思想，有效处理非对称问题。

问题描述
考虑非对称线性方程组：Ax = b，其中A是n×n非对称矩阵，b是已知向量。当A规模较大且稀疏时，直接法求解成本高昂。QMR算法通过构建Krylov子空间来寻找近似解，特别适合这类问题。

算法原理
QMR的核心思想是：在Krylov子空间K_m(A, r₀) = span{r₀, Ar₀, ..., A^(m-1)r₀}中寻找近似解x_m，其中r₀ = b - Ax₀是初始残差。通过Lanczos双正交化过程生成基底，并构建三对角矩阵，将原问题转化为小型最小二乘问题。

详细步骤

步骤1：初始化

1. 选择初始猜测解x₀，计算初始残差r₀ = b - Ax₀
2. 设置β = ||r₀||₂，v₁ = r₀/β（第一个Krylov向量）
3. 选择任意向量w₁满足(w₁, v₁) = 1（通常取w₁ = v₁/||v₁||₂²）
4. 初始化ξ = β，η₀ = 0，c₀ = -1，s₀ = 0
5. 设置ρ₀ = 1，θ₀ = 0，d₀ = 0

步骤2：Lanczos双正交化过程
对于j = 1, 2, ..., m执行：

1. 计算α_j = (w_j, Av_j)
2. 计算v̂_{j+1} = Av_j - α_jv_j - β_jv_{j-1}（其中β₁v₀ ≡ 0）
3. 计算ŵ_{j+1} = Aᵀw_j - α_jw_j - δ_jw_{j-1}（其中δ₁w₀ ≡ 0）
4. 设置δ_{j+1} = |(v̂_{j+1}, ŵ_{j+1})|^(1/2)
5. 计算β_{j+1} = (v̂_{j+1}, ŵ_{j+1})/δ_{j+1}
6. 标准化：v_{j+1} = v̂_{j+1}/β_{j+1}，w_{j+1} = ŵ_{j+1}/δ_{j+1}

步骤3：构建三对角矩阵
Lanczos过程生成三对角矩阵T_m：

T_m = [ α₁  β₂   0   ...   0   ]
      [ δ₂  α₂  β₃  ...   0   ]
      [ 0   δ₃  α₃  ...   0   ]
      [ ... ... ... ...  β_m  ]
      [ 0   0   0   δ_m  α_m  ]

步骤4：QMR最小化过程

1. 前向旋转：应用Givens旋转消除T_m的下对角元素

   • 对于j = 1, 2, ..., m：
     • 计算旋转参数：c_j = ρ_{j-1}/σ_j，s_j = δ_{j+1}/σ_j
       其中σ_j = √(ρ²_{j-1} + δ²_{j+1})
     • 更新：ρ_j = √(ρ²_{j-1} + δ²_{j+1})
     • 计算θ_j = s_jα_j，ρ_j = c_jα_j
     • 计算η_j = c_jξ_j，ξ_{j+1} = -s_jξ_j

2. 向后替换：解上三角系统T_m y_m = η（其中η = [η₁, η₂, ..., η_m]ᵀ）

   • 从j = m递减到1：
     • y_j = (η_j - β_{j+1}y_{j+1})/ρ_j（其中β_{m+1} ≡ 0）

步骤5：更新近似解

1. 计算解更新：x_m = x₀ + V_m y_m
   其中V_m = [v₁, v₂, ..., v_m]是Krylov基向量矩阵
2. 计算残差范数：||r_m||₂ = |ξ_{m+1}|

步骤6：收敛检查

1. 如果||r_m||₂ < ε（预设容差），则停止迭代，输出x_m作为解
2. 否则，继续下一次迭代或重启算法

算法特点

• 通过拟最小残差方法避免直接最小化残差范数的高计算成本
• 利用Lanczos双正交化确保数值稳定性
• 适合大规模稀疏非对称线性方程组
• 相比GMRES需要存储更少的基向量（但可能收敛稍慢）

这个算法通过巧妙的数学变换，将原问题转化为易于求解的小型三对角系统，在保持数值稳定性的同时实现了高效计算。