Eigen 中右预条件的显式矩阵与矩阵向量乘法为何迭代差异很大？

2026-11-06数据小许

结论

在精确算术下，如果两种实现使用同一个固定算子，那么显式构造 C = A M^(-1) 和每次计算 C y = A (M^(-1)y)，求解的应是同一个右预条件系统。有限精度会造成差异。如果迭代次数从约 20 次增至 1000 次以上，尤其是 GMRES 下差距明显缩小，建议先检查自定义矩阵向量乘法是否与显式矩阵实现一致，以及 SparseLU::solve() 的输入、输出和误差。

迭代次数不能只靠条件数预测。矩阵的非正规性、特征值分布、残差变化和舍入误差都可能影响 BiCGStab 的收敛。条件数较低，也不保证迭代次数少。

先确认两个算子确实相同

令 y = Mx，右预条件系统为：

C y = b,    C = A M^(-1)

矩阵自由实现对任意输入向量 y 都应计算：

z = M^(-1) y
C_y = A z

也就是先解 Mz = y，再计算 Az。输入必须是当前算子收到的向量 y，不能误用 x、旧的迭代向量或其他缓存值。求解结束后，根据 Mx = y 恢复原变量，即 x = M^(-1)y。

可以选取几个随机向量 v，比较两种实现的结果：

u_explicit = (A M^(-1)) v
z = SparseLU(M).solve(v)
u_matrix_free = A z

计算相对差异，例如 ||u_explicit - u_matrix_free|| / ||u_explicit||。如果差异远大于预期的浮点舍入误差，应先排查算子实现、矩阵维度、向量类型或预条件器求解，再比较迭代性能。

检查 SparseLU 求解

对于固定的 M，通常会先调用一次 compute()，再对不同右端向量多次调用 solve()，复用同一次分解。每次 solve(y) 都会针对当前右端向量求解；如果 M 发生变化，则需要重新计算分解。

同时检查分解和求解状态，并验证解的残差：

Eigen::SparseLU<SparseMatrix> lu;
lu.compute(M);
if (lu.info() != Eigen::Success) {
    // 处理分解失败
}

Vector z = lu.solve(y);
if (lu.info() != Eigen::Success) {
    // 处理求解失败
}

Vector r = M * z - y;

这里的类型名仅作示意，应替换为程序中的稀疏矩阵和向量类型。把 solve() 的结果立即存入具体向量，也有助于避免自定义算子中的延迟求值表达式引用临时对象或已变化的输入。还要确认输入、输出维度和标量类型正确，并确保算子每次调用都执行相同的线性映射。

核对求解容差和残差

比较两种实现时，要确认它们使用相同的初始猜测、停止准则、最大迭代数和精度设置。还要区分求解器检查的预条件系统残差与原系统残差。若迭代得到 y，直接检查 ||b - A y?|| 并不合适，因为变量关系是 x = M^(-1)y。原问题的残差应计算为：

r_original = b - A x
x = M^(-1)y

预条件系统的残差则是：

r_preconditioned = b - A(M^(-1)y)

两者的相对范数不一定相同。不能只根据求解器报告的收敛状态判断原系统是否达到应用所需的精度。

条件数不能单独解释迭代次数

条件数反映了扰动与解敏感性的一部分信息，但不足以可靠预测 BiCGStab 或 GMRES 的收敛速度。矩阵的非正规性、特征值聚集情况，以及有限精度下残差的实际变化，都可能造成不同的收敛表现。因此，A M^(-1) 的条件数在 100–900，并不能保证 BiCGStab 只需少量迭代。

BiCGStab 和 GMRES 的算法特性不同。GMRES 在这个问题上的差距较小，可以作为排查线索，但不能据此认定原因是条件数或 Eigen。应先用同一批向量验证两种算子的输出是否一致，再观察迭代过程中预条件系统残差和原系统残差的变化。

排查顺序

对相同输入向量逐个比较显式算子与矩阵自由算子的结果。
验证 SparseLU 的分解和求解状态，并检查 ||Mz-y|| 的大小。
检查自定义算子的维度、向量类型、输入是否为当前 y，以及是否存在缓存或临时对象生命周期问题。
统一求解器设置，同时记录预条件系统残差和原系统残差。
如果算子一致且 M 的求解残差足够小，再考虑舍入误差、非正规性和算法对问题的敏感性，并用 GMRES 等方法交叉比较。

缺少自定义算子代码和残差记录时，无法断定具体原因。按目前给出的量级，先逐个向量验证矩阵自由算子与显式算子的结果是否一致，最有助于区分问题所在。