形式陈述
考虑可逆的线性系统 公理库 线性方程组 System of linear equations 可写为矩阵方程 Ax=b 的有限个一次方程系统。 A x = b x ∗ = A − 1 b u f x ^ 0 M = A + Δ k 残差 公理库 残差、误差估计与停止准则 Residual and error estimation · Stopping criterion 区分可计算残差与未知真误差，并说明把缺陷转成误差界和停止证书所需的条件。 LU 因子 公理库 Gaussian 消元与 LU 分解 Gaussian elimination · LU factorization · PA equals LU 把 Gaussian 消元保存为可复用的带置换 LU 分解，再用三角求解处理一个或多个右端。
r k = b − A x ^ k , M d k = r k , x ^ k + 1 = x ^ k + d k . 这里 M d k = r k A d k = r k M − 1 A , b
固定因子的精确递推
统一取误差符号 e k = x ∗ − x ^ k 矩阵范数 公理库 矩阵范数与诱导算子范数 Matrix norm · Induced matrix norm · Operator norm of a matrix 用诱导范数和常用可计算矩阵范数度量线性映射的放大能力，并区分算子范数、Frobenius 范数与谱半径。 r k = A e k
e k + 1 = e k − M − 1 A e k = ( I − M − 1 A ) e k = M − 1 Δ e k . 记 T = I − M − 1 A = M − 1 Δ q = ‖ T ‖ q < 1
‖ e k + 1 ‖ ≤ q ‖ e k ‖ , ‖ e k ‖ ≤ q k ‖ e 0 ‖ ⟶ 0. 这说明足够好的固定因子可以反复消除误差，而精确因子 M = A T = 0 q < 1 q ≥ 1 T ρ ( T ) < 1 T k → 0
因子误差还可以给出一个明确但保守的充分条件。令 C = A − 1 Δ a = ‖ C ‖ a < 1 M = A ( I + C ) I + C ‖ ( I + C ) − 1 ‖ ≤ 1 / ( 1 − a )
T = ( I + C ) − 1 C , q ≤ a 1 − a . 因此 a < 1 / 2 ε = ‖ Δ ‖ / ‖ A ‖ 线性系统条件数 公理库 线性方程组的条件数与扰动 Conditioning of linear systems · Matrix condition number 把一般问题条件性具体化为可逆线性系统的右端、系数矩阵与联合扰动界。 a ≤ κ ( A ) ε κ ( A ) ε < 1 / 2 κ ( A ) T
残差、求解和更新扰动下的平台定理
实际计算中，把三类误差定义为以下恒等式中的扰动：
r ^ k = b − A x ^ k + η k , d ^ k = M − 1 r ^ k + σ k , x ^ k + 1 = x ^ k + d ^ k + ξ k . η k σ k ξ k
e k + 1 = T e k − M − 1 η k − σ k − ξ k . 假设每一轮都有非负常数 R 0 , R 1 , W 0 , W 1
‖ η k ‖ ≤ R 0 + R 1 ‖ e k ‖ , ‖ σ k + ξ k ‖ ≤ W 0 + W 1 ‖ e k ‖ . 定义
α = q + ‖ M − 1 ‖ R 1 + W 1 , β = ‖ M − 1 ‖ R 0 + W 0 . 若 α < 1
及 ‖ e k ‖ ≤ α k ‖ e 0 ‖ + β ( 1 − α k ) 1 − α 及 lim sup k → ∞ ‖ e k ‖ ≤ β 1 − α . 证明只需先对精确递推取范数，得到 E k + 1 ≤ α E k + β E k = ‖ e k ‖ E 0 E k ≤ α k E 0 + β ∑ j = 0 k − 1 α j k → ∞ ‖ η k ‖ ≤ R ‖ σ k ‖ ≤ D ‖ ξ k ‖ ≤ U R 1 = W 1 = 0 R 0 = R W 0 = D + U ( ‖ M − 1 ‖ R + D + U ) / ( 1 − q )
直觉 低精度直接解先给出大体正确的答案，高精度残差再指出它代回原方程后还差多少。已有 LU 因子把方程空间里的缺口翻译成解空间中的校正。T η k , σ k , ξ k
精度分工因此有三种角色：用 u f u u r A x ^ k b
图片加载失败 迭代精化的混合精度校正回路 例子与边界 同一个系统，两组因子的相反轨迹
取
A = ( 2 1 1 2 ) , b = ( 3 3 ) , x ∗ = ( 1 1 ) , x ^ 0 = 0. A 3 , 1 κ 2 ( A ) = 3 A − 1 = 1 3 ( 2 − 1 − 1 2 ) κ ∞ ( A ) = 3
第一组因子满足
L g = ( 1 0 1 / 3 1 ) , U g = ( 3 1 0 8 / 3 ) , L g U g = ( 3 1 1 3 ) = A + I = M g . 对任意残差 r = ( r 1 , r 2 ) T L g y = r y 1 = r 1 y 2 = r 2 − r 1 / 3 U g d = y
d 2 = 3 8 y 2 = 3 r 2 − r 1 8 , d 1 = r 1 − d 2 3 = 3 r 1 − r 2 8 . 因此 T g = M g − 1 = 1 8 ( 3 − 1 − 1 3 ) q ∞ = 1 / 2 v = ( 1 , 1 ) T T g v = v / 4
e k = 4 − k v , x ^ k = ( 1 − 4 − k ) v , r k = 3 ⋅ 4 − k v . 范数界 2 − k ‖ e 0 ‖ ∞ 4 − k ( 1 , − 1 ) T 1 / 2
第二组特意提供的劣质因子仍然完全可逆：
L b = ( 1 0 1 / 2 1 ) , U b = ( 2 / 3 1 / 3 0 1 / 2 ) , L b U b = ( 2 / 3 1 / 3 1 / 3 2 / 3 ) = A / 3 = M b . 此时前代给出 y 1 = r 1 y 2 = r 2 − r 1 / 2 d 2 = 2 y 2 = 2 r 2 − r 1 d 1 = 3 2 ( r 1 − d 2 / 3 ) = 2 r 1 − r 2 M b − 1 = 3 A − 1 T b = I − M b − 1 A = − 2 I
e k = ( − 2 ) k v , x ^ k = ( 1 − ( − 2 ) k ) v , r k = 3 ( − 2 ) k v . 下面每个近似解栏列出两个相同坐标中的一个；残差由 r k = A e k
k 良好因子的 x ^ k
‖ r k ‖ ∞ 劣质因子的 x ^ k
‖ r k ‖ ∞
0
0 3 0 3
1
3 / 4 3 / 4 3 6
2
15 / 16 3 / 16 − 3 12
3
63 / 64 3 / 64 9 24
4
255 / 256 3 / 256 − 15 48
两次计算的 A , b T
持续更新偏差留下非零极限
仍取良好因子、x ^ 0 = 0 ξ k = − ϵ v ϵ > 0 v
e k + 1 = 1 4 e k + ϵ v , e k = [ 4 − k + 4 ϵ 3 ( 1 − 4 − k ) ] v ⟶ 4 ϵ 3 v . 这给出了增加轮数仍不能消除的误差。统一范数定理取 α = 1 / 2 β = ϵ 2 ϵ 4 ϵ / 3
Hilbert 系统的混合精度实验
考虑五阶 Hilbert 矩阵：
( H 5 ) i j = 1 i + j − 1 , 1 ≤ i , j ≤ 5. 本次实验使用 NumPy 2.3.5 与 SciPy 1.17.0，可运行完整复算脚本 查看库版本与每轮输出。先将矩阵元素存为 float64（binary64），令参考向量 x ∗ = ( 1 , 1 , 1 , 1 , 1 ) T b = H5 @ x_ref 构造右端。调用 scipy.linalg.lu_factor(H5.astype(np.float32)) 做带部分主元的 binary32 LU 分解，再用 lu_solve 对转为 binary32 的右端求初值，并将初值转回 binary64。
每轮以 binary64 计算 r = b - H5 @ x，将残差转为 binary32 后复用同一组因子调用 lu_solve，最后将校正量转回 binary64 并更新解。下表中的残差范数也在 binary64 中计算。
这里的 x ∗
精化轮数 k
‖ b − H 5 x ^ k ‖ ∞ / ‖ b ‖ ∞ ‖ x ^ k − x ∗ ‖ ∞ / ‖ x ∗ ‖ ∞
0
6.87 × 10 − 8 3.60 × 10 − 3
1
1.94 × 10 − 11 1.02 × 10 − 5
2
6.35 × 10 − 14 2.82 × 10 − 8
3
1.94 × 10 − 16 8.26 × 10 − 11
4
9.72 × 10 − 17 1.78 × 10 − 11
κ 2 ( H 5 ) ≈ 4.77 × 10 5
推论与应用 把抽象扰动界用于浮点程序时，必须为各项提供来源。在标准浮点误差模型 公理库 浮点算术标准误差模型 Standard floating-point arithmetic model 以每次基本运算的小相对扰动和 gamma 记号组织多步浮点误差分析。 ( n + 1 ) u r < 1
‖ η k ‖ ∞ ≤ γ n + 1 ( r ) ( ‖ b ‖ ∞ + ‖ A ‖ ∞ ‖ x ^ k ‖ ∞ ) , γ n + 1 ( r ) = ( n + 1 ) u r 1 − ( n + 1 ) u r . 由 ‖ x ^ k ‖ ∞ ≤ ‖ x ∗ ‖ ∞ + ‖ e k ‖ ∞
R 0 = γ n + 1 ( r ) ( ‖ b ‖ ∞ + ‖ A ‖ ∞ ‖ x ∗ ‖ ∞ ) , R 1 = γ n + 1 ( r ) ‖ A ‖ ∞ , 再单独加上精度转换的误差界。三角求解与更新也必须给出相应的 W 0 , W 1 M k ‖ I − M k − 1 A ‖ ‖ M k − 1 ‖ M
混合精度讨论中常见 κ ( A ) u f < 1 u κ ( A ) u r κ ( A ) ε < 1 / 2 u f ε u f
停止时可使用尺度化 normwise 后向误差
berr k = ‖ b − A x ^ k ‖ ∞ ‖ A ‖ ∞ ‖ x ^ k ‖ ∞ + ‖ b ‖ ∞ , 并检查校正量相对于解的大小。分子、分母同时为零时可约定该指标为零。达到容差后，应以显式高精度残差复核；连续数轮不降、校正量增长、因子求解失败、非有限值或达到最大轮数，则应返回相应退出状态。过去几轮校正量之比可以帮助诊断，却不能自动充当约束所有后续轮次的收缩常数。残差与误差估计 公理库 残差、误差估计与停止准则 Residual and error estimation · Stopping criterion 区分可计算残差与未知真误差，并说明把缺陷转成误差界和停止证书所需的条件。 线性系统条件数 公理库 线性方程组的条件数与扰动 Conditioning of linear systems · Matrix condition number 把一般问题条件性具体化为可逆线性系统的右端、系数矩阵与联合扰动界。
稠密问题只做一次 O ( n 3 ) O ( n 2 ) nnz ( A )
参考资料
Erin Carson and Nicholas J. Higham, “Accelerating the Solution of Linear Systems by Iterative Refinement in Three Precisions” , SIAM Journal on Scientific Computing 40(2), 2018, A817–A847；§§2–3 区分求解器假设、收缩和误差平台，§7 给出 LU 专门化分析。
James Demmel et al., “Error Bounds from Extra Precise Iterative Refinement” , LAPACK Working Note 165, March 2004，§2 的残差、求解与更新扰动分析；这是 2006 年 ACM Transactions on Mathematical Software 论文之前的报告版本。
Nicholas J. Higham, “What Is Iterative Refinement?” , 2023。
Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms , 2nd ed., SIAM, 2002, Ch. 12.
LAPACK Users’ Guide, 3rd ed., SIAM, 1999, Further Details: Error Bounds for Linear Equations .
