后向误差分析

后向误差分析（backward error analysis）把计算结果看成一个邻近问题的精确解，并衡量为了得到它，需要把输入改动多少。前向误差问结果离原问题答案有多远；后向误差问结果究竟准确解决了哪一个问题。两者通过问题的条件性联系，不能仅凭小残差就断言答案准确。

残差很小但一个分量完全错误

取 $A = (\begin{matrix} 1 & 0 \\ 0 & 1 0^{- 6} \end{matrix}), b = (\begin{matrix} 1 \\ 1 0^{- 6} \end{matrix}) .$ 精确解为 $x = (1, 1)^{T}$ 。假设计算返回 $\hat{x} = (1, 0)^{T}$ ，其残差为 $r = b - A \hat{x} = (0, 1 0^{- 6})^{T}$ ，二范数只有 $1 0^{- 6}$ ；但相对解误差为 $‖ \hat{x} - x ‖_{2} / ‖ x ‖_{2} = 1 / \sqrt{2}$ 。第二分量已经完全错了。

原因是第二行系数很小：将它的右端从 $1 0^{- 6}$ 改成零，只需很小的整体二范数变化，却使对应解从一变为零。矩阵的二范数条件数为 $1 0^{6}$ ，足以把输入小扰动变成大解误差。

小的右端扰动经病态对角矩阵放大，在第二解分量上产生单位误差 — 横轴为第二个右端分量b₂，纵轴为对应解x₂。关系x₂=10⁶b₂使右端从10⁻⁶变为零时，解从一变为零。

只改右端或只改矩阵

由 $A \hat{x} = b - r$ 可知，计算结果是右端改成 $\tilde{b} = b - r$ 后的精确解。只允许改右端时，相对后向误差就是 $‖ r ‖ / ‖ b ‖$ ，前提 $b \neq 0$ 。这是一种误差模型，不能与允许改动矩阵的模型混用。

若只改矩阵，希望 $(A + Δ A) \hat{x} = b$ ，就必须有 $Δ A \hat{x} = r$ 。对于非零 $\hat{x}$ ，二范数给出下界 $‖ Δ A ‖_{2} \geq ‖ r ‖_{2} / ‖ \hat{x} ‖_{2}$ 。取秩一扰动 $Δ A = \frac{r {\hat{x}}^{T}}{‖ \hat{x} ‖_{2}^{2}}$ 便恰好达到下界，因为它将 $\hat{x}$ 映为 $r$ ，且算子二范数为两向量范数的乘积。因此最小相对矩阵后向误差为 $η_{A} = \frac{‖ r ‖_{2}}{‖ A ‖_{2} ‖ \hat{x} ‖_{2}} .$ 主例等于 $1 0^{- 6}$ 。这份扰动一般不保持矩阵的对称性、稀疏模式或物理约束；若要求保结构，需另定义结构化后向误差。

同时允许矩阵与右端扰动

定义最小 $η$ ，使存在 $Δ A, Δ b$ 满足 $(A + Δ A) \hat{x} = b + Δ b, ‖ Δ A ‖_{2} \leq η ‖ A ‖_{2}, ‖ Δ b ‖_{2} \leq η ‖ b ‖_{2} .$ 由 $r = Δ A \hat{x} - Δ b$ 和三角不等式，有 $‖ r ‖_{2} \leq η (‖ A ‖_{2} ‖ \hat{x} ‖_{2} + ‖ b ‖_{2}) .$ 这个下界也能达到。记 $D = ‖ A ‖_{2} ‖ \hat{x} ‖_{2} + ‖ b ‖_{2}$ ， $α = ‖ A ‖_{2} ‖ \hat{x} ‖_{2} / D$ 。取 $Δ A = α \frac{r {\hat{x}}^{T}}{‖ \hat{x} ‖_{2}^{2}}, Δ b = - (1 - α) r .$ 两项相加满足残差关系，其相对范数均为 $‖ r ‖_{2} / D$ ，故 $η = \frac{‖ b - A \hat{x} ‖_{2}}{‖ A ‖_{2} ‖ \hat{x} ‖_{2} + ‖ b ‖_{2}} .$ 这里先假定 $A, b, \hat{x}$ 的所需分母非零；退化情况应按原最小扰动定义单独处理。主例该数约为 $5 \times 1 0^{- 7}$ ，仍远小于前向误差。

条件数怎样连接前后向误差

设 $A$ 非奇异且 $b \neq 0$ ，因此真解非零。只改右端时， $A (\hat{x} - x) = - r$ ，所以 $‖ \hat{x} - x ‖ \leq ‖ A^{- 1} ‖ ‖ r ‖$ 。再用 $‖ b ‖ = ‖ A x ‖ \leq ‖ A ‖ ‖ x ‖$ ，得到 $\frac{‖ \hat{x} - x ‖}{‖ x ‖} \leq κ (A) \frac{‖ r ‖}{‖ b ‖} .$ 因此小后向误差与适度条件数结合，才能保证小前向误差。主例的条件数巨大，前向保证弱；如果变式取 $A = I$ 、 $b = (1, 1)$ 、 $\hat{x} = (1, 1 + 1 0^{- 6})$ ，残差与解误差范数恰好相同，右端相对后向误差和前向相对误差都为 $1 0^{- 6} / \sqrt{2}$ 。

一个算法若在舍入误差尺度内总能解释为轻微扰动输入的精确求解，称为后向稳定。它仍可能在病态问题上产生明显前向误差。反过来，某个特定算例碰巧准确，也不能证明整个算法后向稳定。LAPACK 的线性方程误差章节区分后向误差、条件数和解误差界，其中的范数后向误差采用无穷范数，另给逐分量的误差模型；上面采用二范数的版本则已通过下界与扰动构造证明，不能在不同范数之间直接搬用数值。^[1]

可运行的误差诊断

代码使用二范数，同时返回残差、固定右端模型和允许两种扰动的后向误差，以及给定参考解时的前向误差。它不修改输入；计算谱范数使用稠密奇异值算法，约为 $O (n^{3})$ ，残差乘法为 $O (n^{2})$ 。大型稀疏程序通常用范数或条件数估计器，不能忽略诊断本身的成本。

Python 3

import numpy as np


def linear_error_diagnostics(A, b, approximate, reference=None):
    A = np.array(A, dtype=float, copy=True)
    b, x = np.array(b, dtype=float, copy=True), np.array(approximate, dtype=float, copy=True)
    if b.ndim != 1 or not b.size or x.shape != b.shape or A.shape != (b.size, b.size):
        raise ValueError("matching square system and approximate vector required")
    if not np.isfinite(A).all() or not np.isfinite(b).all() or not np.isfinite(x).all():
        raise ValueError("finite data required")
    def stable_norm(v):
        scale = float(np.max(np.abs(v)))
        if not np.isfinite(scale):
            raise ArithmeticError("vector outside floating-point range")
        if scale == 0:
            return 0.0
        norm = scale*float(np.sqrt(np.sum((v/scale)**2)))
        if not np.isfinite(norm):
            raise ArithmeticError("norm outside floating-point range")
        return norm
    r = b-A@x
    nr, nb, nx = stable_norm(r), stable_norm(b), stable_norm(x)
    na = float(np.linalg.norm(A, 2)); denominator = na*nx+nb
    if not np.isfinite(na) or not np.isfinite(denominator):
        raise ArithmeticError("diagnostic scale outside floating-point range")
    result = {"residual_norm": float(nr),
              "rhs_backward": float(nr/nb) if nb > 0 else (0.0 if nr == 0 else float("inf")),
              "normwise_backward": float(nr/denominator) if denominator > 0 else (0.0 if nr == 0 else float("inf"))}
    if reference is not None:
        exact = np.array(reference, dtype=float, copy=True)
        if exact.shape != x.shape or not np.isfinite(exact).all() or np.max(np.abs(exact)) == 0:
            raise ValueError("finite nonzero matching reference required")
        result["forward_relative"] = stable_norm(x-exact)/stable_norm(exact)
        if not np.isfinite(result["forward_relative"]):
            raise ArithmeticError("forward ratio outside floating-point range")
    if (nb > 0 and not np.isfinite(result["rhs_backward"])) or (denominator > 0 and not np.isfinite(result["normwise_backward"])):
        raise ArithmeticError("backward ratio outside floating-point range")
    return result


print(linear_error_diagnostics([[1, 0], [0, 1e-6]], [1, 1e-6], [1, 0], [1, 1]))
# residual 1e-6; normwise backward about 5e-7; relative forward about 0.7071

分量相对扰动、结构化扰动和范数扰动回答不同问题。主例第二行的右端相对改变百分之百，却在整体范数中很小，正说明整体指标可能掩盖小尺度分量。建模时还应确认数据本身的测量精度，数值求解器无法恢复未被输入可靠记录的信息。

参考资料

↑ LAPACK authors. Error Bounds for Linear Equation Solving 及 Further Details: Error Bounds for Linear Equation Solving，无穷范数后向误差、逐分量后向误差、RCOND、BERR 与 FERR。

[lapack-1] LAPACK authors. Error Bounds for Linear Equation Solving 及 Further Details: Error Bounds for Linear Equation Solving，无穷范数后向误差、逐分量后向误差、RCOND、BERR 与 FERR。

[1]