QR分解仿真误差增大问题及算法类型咨询

我接手的代码里有一个用于求解Ax=B的QR分解子程序，通过在子程序末尾添加的误差计算发现，经过数百步仿真后误差会急剧增大。
这个QR分解用来求解三二次插值的系数——我们先生成5×5×5的数值网格，再推导空间函数形式f(x,y,z)=a₀+a₁x+a₂y+a₃z+…，最终转化为Ax=b的形式：b是125×1的网格点数值向量，A是125×20的坐标计算矩阵，x是20×1的系数向量（由a₀,a₁,a₂,…组成），然后用下面的QR分解算法求解x。

查资料得知部分QR分解方法稳定性较差，但Householder反射法稳定性优异，现在有两个问题：

1. 我没法识别这段代码里的QR分解实现类型，它没做变量归一化，也和网上的实现对不上，请问这是不是Householder反射法？
2. 如果这个实现正确且应该具备稳定性，麻烦给出排查误差随仿真步数增大原因的思路；如果它不是Householder反射法或者实现不佳，用网上的Householder反射法示例能不能提升精度？

! QR decomposition solving for vector B in AX = B
subroutine qr_eqsystem(A, B, X, M, N, error)

implicit none

integer, intent(in) :: m, n
real, intent(in) :: A(m, n), B(m)

real, intent(out) :: x(n)
real, intent(out) :: error

!locals
real :: q(m, n), r(n, n), aux
integer :: i, j, k

do i = 1,n
  do k = 1,m
    q(k, i) = A(k, i)
  enddo
enddo

do i = 1,n
  do j = 1,n
    r(j, i) = 0.0
  enddo
enddo

do i = 1,n
  do j = 1,i-1
    r(j, i) = 0.0
    do k = 1,m
      r(j, i) = r(j, i) + q(k, j)*A(k, i)
    enddo
    
    do k = 1,m
      q(k, i) = q(k, i) - r(j, i)*q(k,j)
    enddo
  enddo
  
  r(i, i) = 0.0
  do k = 1,m
    r(i, i) = r(i, i) + q(k, i)*q(k, i)
  enddo
  r(i, i) = sqrt(r(i, i))
  
  if(r(i,i).ne.0.0)then
    do k = 1,m
      q(k, i) = q(k, i) / r(i, i)
    endif
  endif

enddo

do i = 1,n
  x(i) = 0.0
  do j = 1,m
    x(i) = x(i) + q(j, i)*B(j)
  enddo
enddo

if(r(n, n).ne.0.0)then
  x(n) = x(n) / r(n,n)
endif

do i = n-1,1,-1
  do j = i+1,n
    x(i) = x(i) - r(i, j) * x(j)
  enddo
  
  if(r(i, i).ne.0.0)then
     x(i) = x(i) / r(i, i)
  endif
enddo

error = 0

! calculate AX and compare with B as an error, report total raw error
do i = 1,m
  aux = 0.0
  do j = 1,n
    aux = aux + A(i,j)*X(j)
  enddo
  error = error + abs(aux - B(i))
enddo

endsubroutine qr_eqsystem

问题1：该实现是否为Householder反射法？ #

不是，这个实现是经典Gram-Schmidt正交化方法，和Householder反射法完全不同：

• Householder反射法是通过构造Householder变换矩阵，逐列将A的下三角部分消为0，过程中会修改原矩阵或用辅助数组存储反射因子，不会像这段代码这样逐列对Q进行正交化。
• 这段代码的逻辑是：先把A的列复制到Q，然后对第i列，用前i-1列已经正交化后的Q列向量，将当前Q列投影正交化，再归一化得到Q的第i列，同时计算R矩阵的元素——这完全是经典Gram-Schmidt的步骤。
• 经典Gram-Schmidt的数值稳定性很差，尤其是当矩阵A的列向量接近线性相关时，误差会快速累积，这和你观察到的“数百步仿真后误差急剧增大”现象高度吻合。

问题2：误差排查思路与替换方案 #

替换为Householder反射法的效果

肯定能提升精度。Householder反射法是数值稳定的QR分解实现，它的误差增长是可控的，远优于经典Gram-Schmidt。只要网上的Householder示例是正确的（针对超定方程组的最小二乘求解，因为你的问题是125个方程解20个变量，属于超定），替换后能显著降低每一步的求解误差，进而缓解仿真步数增加后的误差累积。

若假设当前实现正确，误差增大的排查思路

如果暂时不换实现，排查方向可以从这几个角度入手：

• 插值模型本身的问题：三二次插值的基函数是否和网格点匹配？比如5×5×5的网格，三二次插值需要的基函数数量是否确实是20？有没有可能基函数选择不当，导致矩阵A接近秩亏，放大求解误差？
• 仿真过程的误差传递：每一步仿真是不是用前一步的插值结果作为下一步的输入？如果是，误差会逐步累积，即使单步求解误差很小，多步后也会爆发。可以试试单步求解的误差大小，再对比多步后的误差增长速率，判断是不是传递导致的。
• 数值精度问题：代码里用的是单精度real，如果仿真中数值范围变化大，单精度的精度不够会导致误差快速累积。可以改成double precision试试，看误差增长是否减缓。
• 矩阵A的条件数：计算矩阵A的条件数，如果条件数很大（比如1e6以上），说明矩阵接近奇异，不管用什么QR方法，求解误差都会很大。这种情况需要优化插值的基函数或者网格设计，降低条件数。

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