2.2 2.2 最小二乘与 SVD 算例：传感器标定工单

2.2 最小二乘与 SVD 算例：传感器标定工单

本节摘要：最小二乘是线性代数在工业现场出场率最高的算法。本节用一条压力传感器标定工单走完完整流程：模型设定、正规方程推导、QR 求解、SVD 诊断病态与秩亏，并展示截断 SVD 如何把一个"解出来但不可信"的拟合救回来。

本节要达成的能力

阅读完本节，你应当能够：

1. 从"让残差平方和最小"出发推导正规方程；
2. 用 SVD 的奇异值谱判断数值秩与病态程度；
3. 实施截断 SVD 与岭回归两种稳定化手段并对比结果。

工单：把电压读数变成真实压力

车间送来一只压力传感器与一套活塞式压力计的标定数据：三十二个压力点对应的输出电压。任务是把"电压到压力"的换算公式拟出来，误差带要写进检定报告。设压力 p 与电压 u 的关系为二次多项式（出厂说明书给了依据）：p 等于 c0 加 c1 乘 u 加 c2 乘 u 平方。

三十二个观测、三个未知系数，这是典型的超定系统：观测数多于未知数，不存在精确解，只能求"最接近"的解。最小二乘的定义就是把残差向量的长度平方压到最小。对目标函数求导置零，得到正规方程：系数矩阵的转置乘系数矩阵，再乘未知向量，等于转置乘观测向量。推导只有一行，但这一行背后有个隐患——转置乘出来的方阵的条件数是原矩阵条件数的平方，病态会被加倍。

import numpy as np

rng = np.random.default_rng(11)
u = np.linspace(0.5, 4.5, 32)                       # 电压读数
p_true = 2.0 + 1.3 * u + 0.08 * u**2
p_obs = p_true + rng.normal(0, 0.05, 32)            # 含检定噪声

A = np.column_stack([np.ones_like(u), u, u**2])     # 设计矩阵

# 路线一：正规方程（教科书解法）
c_ne = np.linalg.solve(A.T @ A, A.T @ p_obs)
# 路线二：QR 分解（数值更稳）
q, r = np.linalg.qr(A)
c_qr = np.linalg.solve(r, q.T @ p_obs)
# 路线三：SVD（自带诊断）
u_s, s, vt = np.linalg.svd(A, full_matrices=False)
c_svd = vt.T @ ((u_s.T @ p_obs) / s)

print("正规方程:", np.round(c_ne, 4))
print("QR      :", np.round(c_qr, 4))
print("SVD     :", np.round(c_svd, 4))
print("设计矩阵条件数:", np.linalg.cond(A))
​

这个良态问题三条路线结果一致到小数点后四位。真正拉开差距的是病态场景。

把病情造出来：病态与秩亏

标定现场常见的坑是采样范围太窄——电压只在一个小区间内取值，常数项与线性项的基函数变得几乎线性相关，设计矩阵的条件数飙升。更极端的是拿多项式阶数往上怼：八阶多项式配三十二个点，高阶基函数彼此纠缠，矩阵接近秩亏。

import numpy as np

rng = np.random.default_rng(5)
u = np.linspace(2.0, 2.6, 32)                        # 窄采样区间！
p_obs = 2.0 + 1.3 * u + rng.normal(0, 0.05, 32)

for deg in [2, 5, 8]:
    A = np.vander(u, deg + 1, increasing=True)
    print(f"阶数 {deg}: 条件数 1e{np.log10(np.linalg.cond(A)):.1f}")
    u_s, s, vt = np.linalg.svd(A, full_matrices=False)
    print("  奇异值谱:", np.array2string(s, precision=1))

    # 完整最小二乘 vs 截断 SVD（丢掉最小奇异值方向）
    c_full = vt.T @ ((u_s.T @ p_obs) / s)
    c_trunc = vt.T @ ((u_s.T @ p_obs) / s[:-1].tolist() + [0])
    u_new = np.linspace(2.0, 2.6, 5)
    A_new = np.vander(u_new, deg + 1, increasing=True)
    spread_full = np.ptp(A_new @ c_full)
    spread_trunc = np.ptp(A_new @ c_trunc)
    print(f"  完整解预测摆幅 {spread_full:.3f}, 截断解 {spread_trunc:.3f}")
​

窄区间加高阶时，奇异值谱会出现断崖式下跌——最小的几个奇异值比最大的小十几个数量级。这些方向上的解分量等于把观测噪声除以一个极小的数再放大，表现为系数巨大、正负交替、预测曲线在采样区间内剧烈摆振。截断 SVD 的处理是直接的：这些方向的分量直接置零，相当于承认"数据里根本没有关于这些方向的信息"。

另一种等价的稳定化是岭回归——正规方程的对角线上加一个小量再求解。两者在数学上是近亲：岭回归是软性地压制小奇异方向，截断 SVD 是硬性切除。工程上我更倾向岭回归加扫参数（第 4 章的交叉验证选岭参数），因为"保留多少"交给数据裁决比人工拍截断阈值更稳。

import numpy as np

def ridge(A, y, lam):
    """岭回归：正规方程加对角正则"""
    n = A.shape[1]
    return np.linalg.solve(A.T @ A + lam * np.eye(n), A.T @ y)

rng = np.random.default_rng(5)
u = np.linspace(2.0, 2.6, 32)
p_obs = 2.0 + 1.3 * u + rng.normal(0, 0.05, 32)
A = np.vander(u, 6, increasing=True)                 # 五阶，病态

for lam in [0, 1e-10, 1e-8, 1e-6]:
    c = ridge(A, p_obs, lam)
    coef_norm = np.max(np.abs(c))
    print(f"lam={lam:.0e}  最大系数模长 {coef_norm:.3e}")
​

从输出能看到一个清晰的相变：岭参数从零加到 10 的负 8 次方，系数模长从天文数字骤降几个数量级，而拟合误差几乎不变。"系数塌缩、拟合不变"就是稳定化起效的标志——被压掉的部分本来就不含信息，只含噪声放大器。

⚠️ 常见坑：直接用通用求解器解正规方程还不检查条件数。转置相乘让条件数平方化，原本 10 的 8 次方的设计矩阵会变成 10 的 16 次方的正规矩阵，双精度下等于宣判解不可信。

💡 关键直觉：SVD 告诉你的不只是"解是多少"，还有"哪些方向上根本没有信息"。奇异值谱的断崖是数据在说话：断崖左侧是信号方向，右侧是噪声放大器。读懂这张谱，比背十个公式有用。

标定工单工作流

本节要点回顾

• 正规方程推导自残差平方和最小化，但条件数平方化是它的原罪；
• QR 路线数值更稳，是良态超定问题的默认选择；
• 奇异值谱断崖标识无信息方向，截断 SVD 把噪声放大器置零；
• 岭回归是软性正则，配交叉验证选参数比人工截断更可靠；
• **"系数塌缩、拟合不变"**是稳定化起效的判据，写检定报告时要附上。