QR分解成功但误差过高的原因排查求助

问题背景 #

我在Fortran中实现两种QR分解方案求解超定方程组Ax=B：

• 自定义QR分解子程序qr_eqsystem
• 调用Intel Fortran内置LAPACK的dgeqrf、dormqr和dtrtrs函数

其中A为125×21矩阵（由空间坐标构建，模拟过程中相对恒定），B为125×1向量（元素幅值及极差随模拟时间增大），x为21×1待求解向量。

模拟初期，两种方案的残差（sum(abs(matmul(A,x)-B))）均低于1e-5，但随模拟推进残差持续升高且无饱和趋势。已验证分解得到的Q、R矩阵满足Q^TQ近似为单位阵、Q*R可还原A，但无法解释残差升高的原因。此外，MATLAB的x=A\B方案残差更大，尽管三者分解得到的Q、R矩阵几乎一致。

可能的原因分析 #

1. 残差指标未做归一化处理 #

当B的元素幅值和极差随模拟持续增大时，绝对残差sum(abs(Ax-B))会随B的量级同步升高——这并非代表求解精度恶化，只是指标未消除B本身的量级影响。建议改用相对残差评估求解质量：

relative_residual = sum(abs(matmul(A,x)-B)) / sum(abs(B))

或基于2-范数的相对残差norm(Ax-B)/norm(B)，该指标能更准确反映求解的实际精度。

2. 浮点舍入误差随B量级放大 #

A是125×21的超定矩阵，QR分解求解的是最小二乘解x = argmin ||Ax-B||_2：

• 即使Q、R的正交性和A还原精度在可接受范围内，计算Q^TB时的舍入误差会随B的幅值增大而累积；后续回代求解Rx=Q^TB时，这些误差会进一步传递到x中，最终放大Ax-B的残差。
• 若A存在列相关性（秩接近21但不满秩），最小二乘解的条件数会显著增大，微小的B扰动或计算误差都会导致x的显著偏差，进而推高残差。

3. 模拟迭代中的误差累积循环 #

如果B是模拟过程中逐步迭代生成的（比如前一步求解的x参与后续B的计算），那么每一步的x求解误差会传递到下一轮的B中，形成误差累积循环：

• 单步求解的残差可能相对较小，但多次迭代后误差会被不断放大，最终表现为残差随时间持续升高且无饱和趋势。

4. MATLAB A\B的实现差异 #

MATLAB的A\B对于超定方程组默认使用QR分解，但可能采用了不同的数值处理策略：

• 比如对秩亏矩阵的截断规则、中间计算的精度优化（如使用更高精度的临时变量）。即使Q、R矩阵视觉上“几乎一致”，微小的数值差异在B量级较大时会被放大，导致残差表现更明显。

内容的提问来源于stack exchange，提问作者Jesse Feng