数值计算、条件数与可靠实验

方程在理想数系上定义运算；数值程序存储有限近似，执行一系列舍入运算，再返回方程答案的近似。区分输入表示误差、算法算术误差和模型本身的近似。提高精度可以减少舍入，却不改变错误导数、病态逆问题或统计偏差。调试小数前，先确定数学目标、尺度，以及应用真正关心的误差。

正规二进制浮点数包含符号、固定 p 位有效二进制数字及有限范围内的指数。在半开区间 [2^e,2^(e+1)) 内，相邻正规数间距为 2^(e−p+1)。间距随量级增加；同一绝对增量能改变小数，却可能在大数加法中消失。指数控制范围，有效数字控制局部相对精度，两者解决不同存储问题。

常见 IEEE binary16、binary32、binary64 的 p 分别为 11、24、53，包括正规数隐含的首位。NumPy 的 eps 指一上方的间距 2^(1−p)。采用最近舍入且中点取偶时，单位舍入误差 u=eps/2。eps 不是全局最小可分辨变化：间距随量级改变，二的幂上下间距也不同。转换或归约后应检查实际 dtype。

对结果为正规数且没有上溢、下溢的基本运算，标准模型为 fl(a op b)=(a op b)(1+δ)，|δ|≤u。它描述存储输入上精确运算的舍入，不保证存储 a 等于原始实测量，也不保证长算法总相对误差至多 u。精确结果为零时更适合绝对误差分析。接近下溢范围，次正规数允许渐进的小绝对间距，却可能失去相对精度。

上溢是结果超过有限范围，可能产生无穷；后续无穷减无穷或无穷除无穷产生 NaN。下溢可能得到次正规数或零。因此模型赋予正质量的极小概率仍可能存成零。NaN 不能当作普通缺失数参与比较，许多比较会返回假。在约定边界检查有限性，并区分模块 27 中缺失预测支撑导致的数学无穷，与算法失败导致的无穷。

将已舍入 binary32 数组转成 binary64，只保留舍入后的值，不能恢复失去的原始位。先用 binary64 累加、最后才存 binary32，则可能保留低精度累加会丢掉的信息。分别声明输入、累加和输出精度。混合精度需要逐阶段的误差论证，不能只因某处出现大 dtype 名称就保证全流程准确。

参考值 q 与近似 q_hat 的绝对误差为 |q_hat−q|；q 非零时除以 |q| 得相对误差。10^(−8) 相对百万很小，相对 10^(−12) 很大。零目标没有相对误差分母，应按应用尺度选择绝对误差。向量要声明范数；分量检查能发现整体相对范数掩盖的小分量错误。

减法并不必然不准。接近的存储浮点数之差甚至可能精确表示。危险在于操作数已有误差，而所求真差很小；扰动差就可能主导答案。当 |a−b| 远小于 |a|+|b|，即使输入有很多准确位，也未必能确定差的准确位。这种敏感性称消去，不能简单归因于减法指令本身。

考虑小正 x 的 sqrt(1+x)−1：先形成 1+x 可能抹去 x，再减一得到零。乘共轭得到 x/(sqrt(1+x)+1)，在定义域内是同一实函数。分母接近二，避免近数相减，并且不会先把小 x 加入一。实数代数等价式可以有不同数值行为；改写可消除不必要算法误差而不改变目标。

n 个存储输入的和 S，其消去敏感性可用 Σ|x_i|/|S| 度量（S≠0）。该比值大时，各输入微扰可以大幅改变最终相对答案。在通常舍入模型下，顺序求和绝对误差至多 γ_(n−1)Σ|x_i|，其中 γ_k=ku/(1−ku)、ku<1。这是带条件的标准界，不能套用于上溢或任意执行规则。除以 |S| 就解释了消去和为何可能有巨大相对误差。

配对归约将运算依赖深度由 n 降至约 log2 n，在相应条件下改善误差界。补偿方法保存普通加法丢弃的低位。Python math.fsum 可高质量求存储二进制输入的和，却不知道原先十进制意图。线程归约等不同顺序可能改变末位。按量级重排并非通用消去解法，不同归约也不能修复数据问题本身的病态。

数值比较要声明参考和尺度。常用规则 |a−b|≤atol+rtol|b| 将 b 作为参考；NumPy allclose 使用这种不对称形式。接近零时 atol 主导。不同物理单位应分开定容差，不能用一个神奇无量纲常数。小 rtol 配大 atol 可能接受完全错误的微小概率。整数计数、任务 ID、形状和布尔约定若数学要求精确，就应精确检查。

边界与分支尤其要小心。略负特征值既可能来自 PSD 矩阵舍入，也可能反映真正不定输入。结合尺度容差比较，公开容差，并区分数值证据与数学证明。静默裁剪不能使域外输入有效；在取 log 前裁剪概率会把无穷或很大损失改成另一个量，不能称为原目标的稳定求值。

条件性属于带明确输入输出度量的数学问题。可微标量 f 在 x、f(x) 都非零时，局部相对条件数为 |x f’(x)/f(x)|，一阶描述相对输入变化到相对输出变化的放大。x² 的条件数为二，x−1 接近一时为 |x/(x−1)|，可极大。零或不可微点需其他绝对或局部定义。

度量本身很重要。混合米、秒和无量纲系数的向量，没有天然合理的未加权欧氏误差。改变单位可能改变矩阵条件数，却不改变物理实验。先选无量纲坐标或有应用理由的权重。如果科学目标是预测，应单独考察预测敏感性：近共线特征允许系数间不稳定交换，而在观测设计上的合成预测可能相似。新输入离开原设计后又可能暴露不稳定方向。因此应报告目标量，不能以一个条件数裁定整个模型。

算法稳定性描述特定计算对算术误差的响应。前向误差比较输出与预定存储输入的精确答案；后向误差问输出是否精确解了一个邻近问题，以及输入改变多少。病态问题上，后向稳定算法仍可有很大前向误差：它忠实解出的邻近输入可能有很不同的数学解。准确算术与不敏感数学是独立性质。

对非奇异 A 的 Ax=b，二范数条件数 κ2(A)=||A||2||A^(−1)||2=σ_max/σ_min。将 A、b 整体乘同一非零常数不改变该数。非均匀行列缩放可以改变坐标度量和条件性，必须说明怎样恢复原模型。奇异矩阵无普通逆，条件数无穷。数值秩依赖声明的阈值和精度。

固定 A，扰动 b 为 δb，则 δx=A^(−1)δb，故 ||δx||≤||A^(−1)||||δb||。再用 ||b||≤||A||||x||，得到相对解误差至多 κ(A) 倍相对右端扰动。这是最坏方向上界，不是每次扰动等号。沿最小奇异模态的变化放大最多；大条件数说明可能敏感，不证明每次计算必然不准确。

若 A 也改变，写 (A+ΔA)(x+δx)=b+δb，整理 δx=A^(−1)(δb−ΔA x−ΔA δx)。取范数并移项。在 κ ε_A<1 时，ε_A=||ΔA||/||A||、ε_b=||δb||/||b||，相对误差界为 κ(ε_A+ε_b)/(1−κ ε_A)。分母说明何时扰动论证有用；近奇异处，极小变化也可能使这个看似放心的界失效。

残差 r=b−A x_hat 无需真 x 即可算。固定 A 时，x_hat 精确解 b−r，所以 b 非零时 ||r||/||b|| 是右端后向误差度量。前向相对误差≤κ(A)||r||/||b||。因此小相对残差要结合条件数才说明解准。应以充分精度算残差；若沿用低精度计算，可能隐藏差异。

方形非奇异系统应采用适合的分解求解，而非显式形成逆再乘 b。逆需要额外运算和存储，又多一层舍入。一般稠密求解通常用带主元消去；对称正定系统可用 Cholesky。结构要建立，不能由变量名推断。分解失败可能来自无效数学假设、舍入后近奇异输入，或应分析的阈值判断。

最小二乘最小化 ||Xw−y||²，X 有 n 行 d 列。满列秩时正规方程为 XᵀXw=Xᵀy。SVD X=UΣVᵀ 给 Gram 的正特征值 σ_i²，故实数满秩矩阵精确满足 κ2(XᵀX)=κ2(X)²。形成 Gram 时可能已经丢掉小模态，随后求解再准也无法恢复。舍入改变矩阵后，计算条件数未必等于原条件数精确平方。

约化 QR 写 X=QR，Q 列正交归一。将 y 分成 QQᵀy 及正交剩余，勾股定理使目标等价于 ||Rw−Qᵀy||² 加一个常数。R 非奇异时用三角回代解 Rw=Qᵀy，避免形成 Gram 后显式平方奇异值。后向稳定性属于合适实现与条件，不属于所有自称 QR 的手写例程。

SVD 直接展示秩和敏感方向。满秩解 VΣ^(−1)Uᵀy 中，小 σ_i 将投影噪声乘大倒数。秩亏时伪逆在指定阈值下返回最小范数解。截断小奇异值改变目标，可以作为正则化，却必须说明。固定 NumPy lstsq 显式用 rcond=None，返回秩和奇异值；残差数组可能为空，应直接算残差以统一诊断。

大型系统可用迭代避免稠密分解。SPD Hessian H 的二次梯度迭代误差 e_(t+1)=(I−ηH)e_t。每个模态乘 1−ηh_i；固定步长严格收缩全部模态当且仅当 0<η<2/h_max。h_min 很小时，安全步长仍可能使该模态慢。共轭梯度利用 SPD 结构，精确算术下至多维数步终止；有限算术、停止容差和预条件会限定此结论。

预条件改变系统，使迭代进展更均衡，同时保留恢复原解的方式。SPD 问题适当对称变换保留求解结构。学习中的标准化也改变坐标并改善曲率不平衡，但若不对应转换惩罚，新坐标惩罚代表不同几何。声明坐标后测条件与收敛；岭这样改变估计量的方法，不应与纯代数求解修复混淆。

成本也依赖结构。n≥d 稠密 QR 通常为 nd² 量级工作、nd 设计存储；三角求解为 d²。稠密方形分解为 d³，复用分解解其他右端则便宜。稀疏矩阵向量乘可按非零数计成本，但稀疏分解可能产生填充。比较理论操作数与时间前先声明模型和内存形状；硬件、库影响常数和并行。

σ(z)=1/(1+exp(−z)) 在大负 z 时可内部上溢。非负分支直接算，负分支改为 exp(z)/(1+exp(z))，指数参数都不大于零。有限算术仍可能把微小概率舍入为零，大概率舍入为一。返回合理端点不表示其 log 仍忠实；目标需要 log 概率时，应从 logits 直接求。

有限 logits 取 m=max z_j。公共 exp(m) 抵消，softmax=exp(z_j−m)/Σexp(z_k−m)，各指数至多一，至少一个等于一。log-sum-exp=m+logΣexp(z_k−m)。类别 c 的 NLL 可更准确写成 logΣexp(z_k−m)−(z_c−m)，避免加回大偏移后再次相减。推导要求非空有限 logits；全 −无穷或 NaN 需要另外的输入约定。

Bernoulli y∈{0,1} 的 NLL 为 softplus(z)−yz，稳定形式 max(z,0)−yz+log1p(exp(−|z|))。log1p 在 t 很小时保留 log(1+t) 的信息，不先在 1+t 抹去 t。导数为 σ(z)−y。n×d 的 X 与 d 维 w 的均值损失梯度为 Xᵀ(σ(Xw)−y)/n，再加明确惩罚梯度。稳定算术必须保留推导中同一目标和均值归一化。

梯度检查把解析导数与同一平滑标量目标的独立数值近似比较。坐标 j 的中心差分为 [F(w+he_j)−F(w−he_j)]/(2h)。有界三阶导的 Taylor 展开给 h² 阶截断误差；函数求值舍入被 h 除，通常为相关函数尺度乘 u/h。无限减小 h 反而变差；良好缩放坐标中粗略平衡建议 u^(1/3) 量级，不是通用固定步长。

扫多个 h，同时看绝对和尺度相对误差，检查多个坐标或随机方向。方向检查比较 gᵀv 与沿 v 中心差分，高维时便宜，但要声明 v、h 的尺度。先核对形状、目标归约、正则化、截距，再归咎舍入。遗漏 n 因子会在合理 h 范围持续不符，提高精度不解决。同样错误公式的两份实现可能一致，所以检查独立性重要。

坐标步长也有单位。对百万附近和十亿分之一附近的系数加同一 h，探测的相对变化不同，一个扰动甚至不改变存储值。可在已声明无量纲坐标中从 h_j=h_base max(1,|w_j|) 开始，再扫 h_base，逐项检查。这只是尺度启发式，并非全目标准确性定理。报告实际扰动与有用误差平台；单一通过阈值掩盖截断、舍入和持续代数错误的区别。

在 abs(0) 等尖点，对称差分可为零，普通导数却不存在。分段分支应在边界外检查，除非明确采用次梯度约定。平滑检查时固定抽样、dropout 掩码和数据顺序，否则各函数求值对应不同随机目标。受控小例子的检查给局部实现证据，不证明所有分支输入或优化器收敛。

报告应明确方程、存储输入、敏感阶段 dtype、求解或更新规则，以及诊断容差。解释系数时给条件数或奇异值；区分精确恒等式与测得浮点输出。只显示六位小数不能证明十二位误差界。保留机器可读结果或复现脚本，标注平台相关值。先解释失败，再展示修复。

低精度可能降低内存并提高吞吐，但取决于算术强度、带宽、硬件指令与累加精度。数组存储为 n×d×每元素字节，不含中间量和工作区。转换可能同时保留两份数组；d×d Gram 也可能成瓶颈。对同一问题、停止标准和预算比较数值质量与时间。

近似误差有不同来源和修复。离散化或有限迭代误差可能随细化及更多步下降，舍入却可能随额外操作增加。Monte Carlo 误差按模块 24 抽样律变化；确定性浮点误差依赖算术与条件性；统计估计不确定性来自数据。将它们统称一个“精度”掩盖意义。声明误差是测量、界还是诊断，以及支撑条件。

种子要结合生成器、抽样顺序、软件包和输入生成才标识伪随机实验。它不控制所有线程归约，也不使样本独立。库构建可能选择不同算法，改末位；同一种子在不同生成器接口可生成不同数据。浮点量用有理由的容差回归检查，精确约定用精确检查，并加入概率和、秩假设、残差等结构检查。

记录 Python、包版本、形状、拆分、生成器、阈值、计时是否含生成或分解。相关时区分编译、预热和重复计时，报告时间分布，不能以单次运行作无依据速度声称。教学小型 CPU 例子使推导可核查，不是所有设备精度或最快部署方案的基准证据。

实验 A · 舍入、归约及有理化

使用固定 NumPy 环境，预测每格式大间距边界丢失一。比较顺序和、fsum 与实际二进制输入的 Decimal 参考，解释有理化式；显示舍入不等于数学修复。

下载 lab1_rounding_summation_and_cancellation.py

"""CPU rounding experiments; Decimal references the actual binary inputs."""
import math
from decimal import Decimal, localcontext
import numpy as np

print('NumPy', np.__version__, '; rounding to nearest on the tested CPU')
for dtype in (np.float16, np.float32, np.float64):
    info = np.finfo(dtype)
    large = dtype(2 ** (info.nmant + 1))
    terms = [large, dtype(1), -large]
    total = dtype(0)
    for term in terms:
        total = dtype(total + term)
    print(f'{dtype.__name__}: eps={info.eps:.9e}; spacing(large)={np.spacing(large):g}; sequential={total:g}; exact=1')
    assert total == 0
terms = [1e16, 1., -1e16]
with localcontext() as ctx:
    ctx.prec = 80
    exact = sum(map(Decimal.from_float, terms))
print('binary64 sequential:', sum(terms), '; fsum:', math.fsum(terms), '; binary-input Decimal reference:', exact)
assert math.fsum(terms) == float(exact) == 1
x = 1e-16
naive = math.sqrt(1+x)-1
stable = x/(math.sqrt(1+x)+1)
print(f'sqrt(1+x)-1 at x={x:g}: naive={naive:.9e}; rationalised={stable:.9e}')
assert naive == 0 and stable > 0
print('0.1+0.2:', repr(.1+.2), '; difference from float(0.3):', (.1+.2)-.3)
print('Subnormal != normal minimum:', np.finfo(np.float64).smallest_subnormal < np.finfo(np.float64).tiny)

NumPy 1.26.4 ; rounding to nearest on the tested CPU
float16: eps=9.765625000e-04; spacing(large)=2; sequential=0; exact=1
float32: eps=1.192092896e-07; spacing(large)=2; sequential=0; exact=1
float64: eps=2.220446049e-16; spacing(large)=2; sequential=0; exact=1
binary64 sequential: 0.0 ; fsum: 1.0 ; binary-input Decimal reference: 1
sqrt(1+x)-1 at x=1e-16: naive=0.000000000e+00; rationalised=5.000000000e-17
0.1+0.2: 0.30000000000000004 ; difference from float(0.3): 5.551115123125783e-17
Subnormal != normal minimum: True

实验 B · 敏感性、残差与求解器

计算小右端扰动及放大的解变化；比较精确生成数据的 QR、lstsq 和正规方程。系数数值误差可随库变化；满秩条件性恒等式和小残差反例是数学重点。

下载 lab2_solves_conditioning_and_residuals.py

"""Small deterministic systems: residual, sensitivity and solver choice."""
import numpy as np

epsilon = 1e-8
A = np.array([[1., 1.], [1., 1.+epsilon]])
truth = np.ones(2)
b = A @ truth
delta = np.array([0., 1e-8])
x = np.linalg.solve(A, b+delta)
relative_input = np.linalg.norm(delta)/np.linalg.norm(b)
relative_error = np.linalg.norm(x-truth)/np.linalg.norm(truth)
print(f'cond2(A)={np.linalg.cond(A):.6e}; relative RHS perturbation={relative_input:.6e}; relative solution change={relative_error:.6e}')
wrong = np.array([0., 2.])
print(f'wrong x relative residual={np.linalg.norm(A@wrong-b)/np.linalg.norm(b):.6e}; relative error={np.linalg.norm(wrong-truth)/np.linalg.norm(truth):.6e}')
assert relative_input < 1e-8 and relative_error > .9
t = np.linspace(-1, 1, 30)
X = np.column_stack([np.ones_like(t), t, t+1e-7*t*t])
beta = np.array([1., 2., -1.])
y = X @ beta
Q, R = np.linalg.qr(X, mode='reduced')
qr = np.linalg.solve(R, Q.T@y)
svd, _, rank, singular = np.linalg.lstsq(X, y, rcond=None)
print(f'X rank={rank}; cond2(X)={singular[0]/singular[-1]:.6e}; computed cond2(Gram)={np.linalg.cond(X.T@X):.6e}')
for label, estimate in [('QR', qr), ('lstsq', svd)]:
    print(f'{label}: coefficient error={np.linalg.norm(estimate-beta):.6e}; residual={np.linalg.norm(X@estimate-y):.6e}')
    assert np.linalg.norm(estimate-beta) < 1e-6
try:
    normal = np.linalg.solve(X.T@X, X.T@y)
    print(f'normal equations: coefficient error={np.linalg.norm(normal-beta):.6e}; residual={np.linalg.norm(X@normal-y):.6e}')
except np.linalg.LinAlgError:
    print('normal equations: rounded Gram reported singular')
print('Exact full-column-rank identity: cond2(Gram)=cond2(X)^2; computed Gram may have lost its smallest mode.')
H = np.diag([1., 100.])
for rate in (.01, .03):
    w = np.ones(2)
    for _ in range(20):
        w -= rate * (H@w)
    print(f'quadratic eta={rate:.2f}; w20={w}; stable={0 < rate < 2/100}')

cond2(A)=4.000000e+08; relative RHS perturbation=3.535534e-09; relative solution change=1.000000e+00
wrong x relative residual=3.535534e-09; relative error=1.000000e+00
X rank=3; cond2(X)=4.444812e+07; computed cond2(Gram)=2.444580e+15
QR: coefficient error=3.437643e-09; residual=1.464485e-15
lstsq: coefficient error=2.904832e-09; residual=6.466036e-15
normal equations: coefficient error=7.856742e-02; residual=9.681657e-09
Exact full-column-rank identity: cond2(Gram)=cond2(X)^2; computed Gram may have lost its smallest mode.
quadratic eta=0.01; w20=[0.81790694 0.        ]; stable=True
quadratic eta=0.03; w20=[5.43794343e-01 1.04857600e+06]; stable=False

实验 C · 稳定损失与导数扫描

修复原始指数上溢，在同一均值目标比较三个差分步长。最小步长因舍入失败。单独诊断遗漏均值因子，解释绝对值不可微例子。

下载 lab3_stable_losses_and_gradient_checks.py

"""Repair overflow and audit a smooth central-difference gradient check."""
import numpy as np

def sigmoid(z):
    z = np.asarray(z, dtype=np.float64)
    out = np.empty_like(z)
    positive = z >= 0
    out[positive] = 1/(1+np.exp(-z[positive]))
    e = np.exp(z[~positive])
    out[~positive] = e/(1+e)
    return out

logits = np.array([1000., 1001., 1002.])
with np.errstate(over='ignore', invalid='ignore'):
    broken = np.exp(logits)/np.exp(logits).sum()
shifted = logits-logits.max()
softmax = np.exp(shifted)/np.exp(shifted).sum()
loss = np.log(np.exp(shifted).sum())-shifted[0]
print('raw softmax finite:', bool(np.isfinite(broken).all()))
print('shifted softmax:', softmax, '; class-0 NLL:', f'{loss:.9f}')
assert np.isfinite(softmax).all() and np.isclose(softmax.sum(), 1)
extreme = np.array([-1000., 1000.])
labels = np.array([1., 0.])
losses = np.maximum(extreme, 0)-labels*extreme+np.log1p(np.exp(-np.abs(extreme)))
print('stable sigmoid:', sigmoid(extreme), '; stable Bernoulli losses:', losses)
assert np.all(losses == 1000)
X = np.array([[1., .2], [1., -1.], [1., 2.]])
y = np.array([1., 0., 1.])
w = np.array([.3, -.4])
def objective(v):
    z = X@v
    return np.mean(np.logaddexp(0, z)-y*z)
analytic = X.T@(sigmoid(X@w)-y)/len(y)
for step in (1e-2, 1e-5, 1e-16):
    numeric = np.array([(objective(w+step*np.eye(2)[j])-objective(w-step*np.eye(2)[j]))/(2*step) for j in range(2)])
    error = np.linalg.norm(numeric-analytic)
    print(f'h={step:.0e}; numerical={numeric}; error={error:.6e}')
    if step == 1e-5:
        assert error < 1e-9
print('missing mean factor discrepancy:', np.linalg.norm(3*analytic-analytic))
print('abs at zero: central difference=0 while derivative does not exist; a symmetric slope is not a differentiability certificate.')

raw softmax finite: False
shifted softmax: [0.09003057 0.24472847 0.66524096] ; class-0 NLL: 2.407605964
stable sigmoid: [0. 1.] ; stable Bernoulli losses: [1000. 1000.]
h=1e-02; numerical=[-0.13316435 -0.66738056]; error=2.981183e-06
h=1e-05; numerical=[-0.13316411 -0.66738353]; error=1.168331e-11
h=1e-16; numerical=[0. 0.]; error=6.805391e-01
missing mean factor discrepancy: 1.3610781824396445
abs at zero: central difference=0 while derivative does not exist; a symmetric slope is not a differentiability certificate.

核查官方固定版 NumPy 1.26 finfo、lstsq、cond 的精度、秩阈值与范数。SciPy 1.14 logsumexp 是选学库比较；实验只需 NumPy。上面的条件性与导数论证在本课推导。

结业任务：给出间距例子、消去修复、残差敏感性反例、求解器理由、稳定概率计算和平滑导数扫描。

可以继续：解释计算建立什么、还剩什么不确定性。泛化下一课连接风险数学与有限数据、核模型。回到课程概览。