2. 线性最小二乘的几何与代数基础

在反演问题中，我们几乎不可能让模型预测与观测数据完全一致。观测有噪声，模型有离散化误差，参数化本身也不完备。因此“求解 Gm = d”往往是一个超定方程组：方程个数多于未知数个数，精确解不存在。最小二乘给出的回答是：在所有候选模型 m 中，选一个让残差 r = d − Gm 的平方和最小的。要真正理解这个“最小”意味着什么、什么时候有解、解是否唯一，必须回到几何与代数两条线索上。

一、几何视角：把数据看成向量空间中的点

设数据向量 d ∈ R^N，模型参数 m ∈ R^M，正演矩阵 G 为 N×M。矩阵 G 的每一列可以看作一个“数据模式”向量，所有列的线性组合构成 G 的列空间 col(G)，它是 R^N 中的一个子空间。当 m 遍历 R^M 时，Gm 恰好扫过整个列空间。于是问题变成：在列空间这个“可达预测集合”里，找一个离观测点 d 最近的点。

欧氏距离下，最近点就是 d 在 col(G) 上的正交投影，记作 d̂ = Gm̂。投影的几何特征是：残差向量 r = d − Gm̂ 必须与列空间中的每一个向量都正交。由于列空间的基就是 G 的各列，这等价于 r 与每一列正交，即 Gᵀr = 0。展开就是 Gᵀ(d − Gm̂) = 0，也就是著名的正规方程 GᵀG m̂ = Gᵀd。几何上的“垂直”直接翻译成了代数上的“左乘转置为零”，这是整个最小二乘理论最核心的一步。

二、代数推导：从正交条件到正规方程

把目标函数写成 Φ(m) = ‖d − Gm‖² = (d − Gm)ᵀ(d − Gm)。展开得 Φ(m) = dᵀd − 2mᵀGᵀd + mᵀGᵀGm。对 m 求梯度并令其为零：∂Φ/∂m = −2Gᵀd + 2GᵀGm = 0，即得正规方程。这里 GᵀG 是 M×M 的对称半正定矩阵，Gᵀd 是 M 维向量。注意正规方程是“方阵系统”，但它未必可逆，这正是存在性与唯一性讨论的入口。

从二阶条件看，Φ 的 Hessian 是 2GᵀG，半正定，所以任何满足正规方程的 m̂ 都是全局极小点，最小二乘问题不会陷入局部极小。这一点在非线性反演中不成立，但在本章的线性情形下是可靠的结论。

三、存在性与唯一性：秩决定一切

正规方程 GᵀG m̂ = Gᵀd 一定有解，因为 Gᵀd 总落在 GᵀG 的列空间内（这是投影的代数体现）。所以最小二乘解总是存在。但唯一性取决于 G 的秩：

• 列满秩（rank(G) = M，即 M ≤ N 且各列线性无关）：GᵀG 正定可逆，解唯一，m̂ = (GᵀG)⁻¹Gᵀd。此时称问题超定且良态。
• 列秩亏（rank(G) = r < M）：GᵀG 奇异，有无穷多个解，它们相差一个零空间向量。此时最小二乘解不唯一，但投影 Gm̂ 和最小残差平方和是唯一的。
• 欠定（M > N）：必然秩亏，通常有无穷多解，需要附加约束（如最小范数解）才能挑出唯一模型。

工程上判断唯一性不能只看 M 与 N 的大小，而要看奇异值。若最小奇异值与最大奇异值之比（条件数倒数）接近机器精度，矩阵在数值上就是秩亏的，直接求逆会放大误差。此时应使用奇异值分解或带阻尼的求解方式，而不是硬解正规方程。

四、一个可手算的示例

设 G = [[1,0],[0,1],[1,1]]，d = [1,2,2]ᵀ。计算 GᵀG = [[2,1],[1,2]]，Gᵀd = [3,4]ᵀ。解 2m₁ + m₂ = 3、m₁ + 2m₂ = 4，得 m̂ = (2/3, 5/3)。预测 Gm̂ = (2/3, 5/3, 7/3)，残差 r = (1/3, 1/3, −1/3)。验证正交性：Gᵀr = (1/3 − 1/3, 1/3 − 1/3) = (0,0)，确实与两列都正交。残差平方和为 1/3。这个例子中 G 列满秩，解唯一。

五、实操中的数值注意点

理论上 m̂ = (GᵀG)⁻¹Gᵀd 很漂亮，但实际计算中直接构造 GᵀG 会把条件数平方。若 G 的条件数是 10⁴，GᵀG 就是 10⁸，在单精度下几乎不可用。更稳妥的做法是对 G 做 QR 分解或 SVD：QR 下解 Rm̂ = Qᵀd，SVD 下用奇异值阈值判断有效秩并构造伪逆。下面是一段示意性代码，用于检查秩与条件数：

import numpy as np

G = np.array([[1.,0.],[0.,1.],[1.,1.]])
d = np.array([1.,2.,2.])

# 正规方程直接求解（仅作对照）
m_ne = np.linalg.solve(G.T @ G, G.T @ d)

# SVD 求解，稳健判断有效秩
U, s, Vt = np.linalg.svd(G, full_matrices=False)
tol = s[0] * max(G.shape) * np.finfo(float).eps
rank = np.sum(s > tol)
m_svd = Vt.T @ np.diag(1/s) @ U.T @ d

r = d - G @ m_svd
print("奇异值:", s, "有效秩:", rank)
print("解:", m_svd, "残差平方和:", r @ r)
print("正交性检查 G^T r:", G.T @ r)

运行这类检查时，重点看三件事：奇异值谱是否出现明显断层（断层说明秩亏）、Gᵀr 是否在数值精度内接近零（验证投影条件）、残差平方和是否与理论最小值一致。若 Gᵀr 明显不为零，通常说明求解过程有数值问题或矩阵构造有误，而不是最小二乘理论失效。

把这一节的几何图像记住，后面讨论残差诊断时就会非常自然：残差不是“剩下的垃圾”，而是观测中无法被模型解释、且与模型预测空间严格正交的那一部分。它的结构、大小和分布，正是判断模型是否合适、数据中是否还有系统误差的重要线索。