最小二乘法

最小二乘法是最古老的优化方法之一（Gauss 1801 年用它确定谷神星轨道），也是线性回归、曲线拟合、传感器融合的数学基础。它是数值优化基础中最经典的特例：当目标函数是残差平方和时，很多情况可以一步求出精确解。

核心直觉

最小二乘的本质用一个画面就能理解：画一条线，让所有数据点到它的”垂直距离平方和”最小。

“平方”的意义：平方让正负残差不互相抵消（否则离线 5 和离线 -5 就抵消成 0 了），同时自然地惩罚大误差——离线越远，惩罚越重（平方增长）。
正规方程 = 一步到位：当模型是参数的线性函数时（如 $y = ax + b$ ），残差平方和是参数的二次函数——而二次函数有唯一最小值，可以一步求出解析解（正规方程 $\hat{\beta} = (X^T X)^{-1} X^T y$ ）。
梯度下降 = 迭代逼近：当模型对参数非线性（如 $y = a \cdot e^{bx}$ ），没法一步求解，就用梯度下降迭代逼近——这就是非线性最小二乘。
为什么是”最小二乘”而不是”最小绝对值”：平方对应高斯噪声下的最大似然估计——如果测量误差服从正态分布，最小二乘解恰好就是统计最优解。

原理详解：从残差到解析解

以下推导记 $X$ 为设计矩阵（ $n \times m$ ）， $y$ 为观测值向量（ $n \times 1$ ）， $\beta$ 为待求参数向量（ $m \times 1$ ）， $r$ 为残差向量。

1. 问题定义

给定 $n$ 个观测点和 $m$ 个参数，线性模型预测：

\hat{y} = X\beta \quad \text{（预测值，} n \times 1 \text{ 向量）} \\ r = y - X\beta \quad \text{（残差 = 观测值 − 预测值）}

目标：找到 $\beta$ 使残差平方和（Residual Sum of Squares, RSS）最小：

S(\beta) = \sum_i r_i^2 = \|y - X\beta\|^2 = (y - X\beta)^T (y - X\beta)

2. 正规方程推导

展开 $S(\beta)$ ：

S(\beta) = y^T y - 2\beta^T X^T y + \beta^T X^T X \beta

对 $\beta$ 求偏导，令其等于零（极值条件）：

\frac{\partial S}{\partial \beta} = -2X^T y + 2X^T X \beta = 0

化简得到正规方程（Normal Equation）：

X^T X \beta = X^T y \quad \text{（} m \times m \text{ 的线性方程组）}

当 $X^T X$ 可逆（满列秩）时，直接解出：

\hat{\beta} = (X^T X)^{-1} X^T y \quad \text{（线性回归的闭式解）}

为什么叫”正规方程”？ “Normal”在这里不是”正常”，而是”正交”——残差向量 $r$ 垂直（正交）于 $X$ 的列空间，即 $X^T r = 0$ 。展开后正好得到 $X^T X \beta = X^T y$ 。这是线性代数中投影的几何本质。

3. 几何解释：投影

最小二乘解的几何含义是：把 $y$ 投影到 $X$ 的列空间（所有可能的 $X\beta$ 构成的子空间）。

        y (观测点，在列空间之外)
       /|
      / |
     /  | r (残差 = y 到投影的垂直距离)
    /   |
   /----+
  X β_hat (投影点 = 最佳预测值，在列空间内)

投影公式： $\hat{\beta} = (X^T X)^{-1} X^T y$ 正是线性代数中的正交投影矩阵 $X(X^T X)^{-1} X^T$ 作用于 $y$ 的结果。残差 $r = y - X\hat{\beta}$ 垂直于列空间——这就是 " $X^T r = 0$ " 的几何含义。

4. 最大似然估计：为什么”二乘”是统计最优

假设观测值 $y = X\beta_{\text{true}} + \varepsilon$ ，其中 $\varepsilon$ 服从均值为 $0$ 、方差为 $\sigma^2$ 的正态分布（即测量噪声是高斯的）。

最大似然估计（MLE）的目标：找到 $\beta$ 使观测到 $y$ 的概率（似然函数）最大。

p(y \mid \beta) = \prod_i \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(y_i - x_i \beta)^2}{2\sigma^2}\right)

取对数似然：

\log p(y \mid \beta) = -\frac{1}{2\sigma^2} \sum_i (y_i - x_i \beta)^2 + \text{const}

最大化 $\log p$ 等价于最小化 $\sum_i (y_i - x_i \beta)^2 = S(\beta)$ 。

结论：当噪声服从高斯分布时，最小二乘解 = 最大似然估计 = 统计最优解。这就是”最小二乘”而不是”最小绝对值”的根本原因——高斯噪声假设在自然界中极其普遍（中心极限定理）。

如果噪声是拉普拉斯分布（重尾），则最优的是最小绝对值（LAD），对应鲁棒回归。

5. 加权最小二乘（WLS）

当不同观测的噪声方差不同（异方差）时，给方差小的观测更大权重：

S_w(\beta) = \sum_i w_i (y_i - x_i \beta)^2 = (y - X\beta)^T W (y - X\beta)

加权正规方程：

\hat{\beta} = (X^T W X)^{-1} X^T W y \quad \text{（} W \text{ 是对角权重矩阵）}

6. 岭回归：解决多重共线性

当特征高度相关时， $X^T X$ 接近奇异（行列式趋零），求逆数值不稳定。岭回归加上 L2 正则化：

S_{\text{ridge}}(\beta) = \|y - X\beta\|^2 + \lambda \|\beta\|^2 \quad \text{（} \lambda \text{ 是正则化强度）}

求导令零得到岭回归正规方程：

\beta_{\text{ridge}} = (X^T X + \lambda I)^{-1} X^T y

加上 $\lambda I$ 后，矩阵恒可逆（只要 $\lambda > 0$ ），数值稳定。 $\lambda$ 越大， $\beta$ 的范数越小（防过拟合），但偏差也越大。

原理图解

线性最小二乘的几何直觉

数据点散布在平面上，最小二乘法找到一条直线，让所有点到线的残差（垂直距离）平方和最小：

非线性最小二乘方法谱系

从解析解到迭代法，方法逐步复杂，能处理的问题也更一般：

Gauss-Newton 用一阶 Taylor 展开把非线性残差”局部线性化”，从而避免计算完整的 Hessian 矩阵；Levenberg-Marquardt（LM） 在 Gauss-Newton 和梯度下降之间自适应切换——离最优远时像梯度下降（稳定），近最优时像 Gauss-Newton（快速收敛），是曲线拟合的工业标配。

动手实践

numpy 手写线性回归（正规方程一步求解）

线性最小二乘有解析解 $\hat{\beta} = (X^T X)^{-1} X^T y$ ：

import numpy as np

# 生成带噪声的线性数据：y = 2x + 1 + noise
np.random.seed(42)
x = np.linspace(0, 10, 50)
X = np.column_stack([x, np.ones_like(x)])   # 设计矩阵 [x, 1]
y = 2 * x + 1 + np.random.randn(50) * 1.5   # 真实斜率 2，截距 1

# 正规方程：β = (XᵀX)⁻¹ Xᵀy，一步求出最优参数
beta = np.linalg.inv(X.T @ X) @ X.T @ y
print(f"斜率={beta[0]:.3f}（真值 2）, 截距={beta[1]:.3f}（真值 1）")
# 也可以直接用 np.linalg.lstsq(X, y) 一行搞定

# 更简洁的写法：numpy 内置最小二乘求解器
beta_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)

三种解法对比：正规方程 vs QR 分解 vs SVD

实际上 np.linalg.lstsq 底层使用的是 SVD 分解，而非直接求逆。这三种方法的数值稳定性和适用场景不同：

import numpy as np

np.random.seed(42)
x = np.linspace(0, 10, 50)
X = np.column_stack([x, np.ones_like(x)])
y = 2 * x + 1 + np.random.randn(50) * 1.5

# ---- 方法 1：正规方程 β = (X^T X)^{-1} X^T y ----
# 优点：直觉清晰。缺点：X^T X 接近奇异时数值不稳定
beta_normal = np.linalg.inv(X.T @ X) @ X.T @ y

# ---- 方法 2：QR 分解 ----
# X = QR，则 X^T X = R^T Q^T Q R = R^T R（因为 Q^T Q = I）
# 正规方程变为 R β = Q^T y，R 是上三角矩阵，回代即可
Q, R = np.linalg.qr(X)
beta_qr = np.linalg.solve(R, Q.T @ y)

# ---- 方法 3：SVD 分解（numpy lstsq 底层使用）----
# X = U Σ V^T，则 β = V Σ^+ U^T y（Σ^+ 是伪逆）
beta_svd, *_ = np.linalg.lstsq(X, y, rcond=None)

print(f"正规方程: 斜率={beta_normal[0]:.4f}")
print(f"QR 分解:  斜率={beta_qr[0]:.4f}")
print(f"SVD:     斜率={beta_svd[0]:.4f}")
# 三者结果几乎相同（2.00 左右），但 SVD 在病态矩阵上最鲁棒

选择建议：

特征数 < 1000 且无共线性：正规方程最快
中等规模通用场景：QR 分解（数值稳定，速度好）
病态矩阵或特征数 > 样本数：SVD（最鲁棒，能处理秩亏）

梯度下降求解线性回归

当特征数极大（如百万维）时，矩阵求逆不可行，改用梯度下降迭代逼近：

import numpy as np

np.random.seed(42)
x = np.linspace(0, 10, 50)
X = np.column_stack([x, np.ones_like(x)])
y = 2 * x + 1 + np.random.randn(50) * 1.5

# 损失函数 S(β) = ||y - X β||² 的梯度
# ∂S/∂β = -2 X^T (y - X β) = 2 X^T (X β - y)
def gradient(X, y, beta):
    return 2 * X.T @ (X @ beta - y)

# 批量梯度下降
beta = np.zeros(2)              # 初始化参数
lr = 0.001                     # 学习率
for step in range(5000):
    grad = gradient(X, y, beta)
    beta = beta - lr * grad     # 参数更新
    if step % 1000 == 0:
        loss = np.sum((y - X @ beta) ** 2)
        print(f"step {step}: loss={loss:.2f}, β=[{beta[0]:.3f}, {beta[1]:.3f}]")
print(f"最终: 斜率={beta[0]:.3f}, 截距={beta[1]:.3f}")
# 与正规方程结果一致，但不需要矩阵求逆

scipy 非线性曲线拟合

当模型对参数非线性时（如指数衰减 $y = a \cdot e^{-bx}$ ），用 scipy.optimize.curve_fit：

from scipy.optimize import curve_fit

def model(x, a, b): return a * np.exp(-b * x)   # 非线性模型
y_nonlin = 5 * np.exp(-0.3 * x) + np.random.randn(50) * 0.2
popt, _ = curve_fit(model, x, y_nonlin, p0=[1, 0.1])  # p0 是初始猜测
print(f"a={popt[0]:.3f}（真值 5）, b={popt[1]:.3f}（真值 0.3）")

numpy 手写 Gauss-Newton 迭代

理解非线性最小二乘的最快方式是从零实现 Gauss-Newton。以指数模型 y = a·e^(-bx) 为例：

import numpy as np

np.random.seed(42)
x = np.linspace(0, 10, 50)
y_true = 5 * np.exp(-0.3 * x)
y_noisy = y_true + np.random.randn(50) * 0.2

def model(params, x):
    """非线性模型: y = a * exp(-b * x)"""
    a, b = params
    return a * np.exp(-b * x)

def jacobian(params, x):
    """雅可比矩阵 J[i][j] = ∂r_i/∂β_j
    r_i = y_i - model(β, x_i)
    ∂r/∂a = -exp(-b*x),  ∂r/∂b = a*x*exp(-b*x)
    """
    a, b = params
    J = np.column_stack([
        -np.exp(-b * x),          # ∂r/∂a
        a * x * np.exp(-b * x),   # ∂r/∂b
    ])
    return J

# Gauss-Newton 迭代
beta = np.array([1.0, 0.1])       # 初始猜测（重要！）
for iteration in range(20):
    r = y_noisy - model(beta, x)  # 残差向量
    J = jacobian(beta, x)         # 雅可比矩阵
    # Gauss-Newton 更新: Δβ = (J^T J)^{-1} J^T r
    delta = np.linalg.solve(J.T @ J, J.T @ r)
    beta = beta + delta
    if np.max(np.abs(delta)) < 1e-8:
        break

print(f"Gauss-Newton: a={beta[0]:.3f}（真值 5）, b={beta[1]:.3f}（真值 0.3）, 迭代 {iteration+1} 次")

Gauss-Newton 的核心思想：把非线性残差 $r(\beta)$ 在当前点做一阶 Taylor 展开， $r(\beta + \Delta\beta) \approx r(\beta) + J\Delta\beta$ ，然后求最小二乘解得到增量 $\Delta\beta$ ——这把非线性问题转化为每一步求解一个线性最小二乘子问题。Levenberg-Marquardt 在此基础上加入阻尼项 $(J^T J + \lambda I)\Delta\beta = J^T r$ ， $\lambda$ 大时退化为梯度下降（稳定但慢）， $\lambda$ 小时接近 Gauss-Newton（快但可能发散），自适应调节 $\lambda$ 就是 LM 算法的精髓。

实践要点

正规方程 vs 梯度下降：特征数少（< 10000）时正规方程一步到位更快；特征数极大时 $(X^T X)^{-1}$ 求逆代价 $O(n^3)$ 太贵，改用梯度下降或 QR 分解。
多重共线性问题：当特征之间高度相关时， $X^T X$ 接近奇异（不可逆），正规方程数值不稳定——这时用岭回归（L2 正则化）或改用 SVD/QR 分解求解。
鲁棒回归：数据中有离群点时，平方损失会被异常值主导。改用 Huber 损失（小误差用平方、大误差用线性）或 RANSAC（随机抽样一致）能获得更鲁棒的拟合。
非线性拟合要给好的初值：curve_fit 的 p0 参数很重要，初值离最优太远会陷入局部最优或不收敛——先画图估算量级是良好习惯。

典型应用

线性回归：最小二乘法的直接应用，统计学和机器学习中最基础的回归方法。
曲线拟合 / 趋势线：Excel 里的”添加趋势线”、科学实验中的经验公式拟合，底层都是最小二乘。
GPS 定位：GPS 接收器同时收到 4+ 颗卫星的距离信号，用最小二乘从超定方程组解出三维坐标 + 时钟偏差——这是你手机定位的数学核心。
相机标定：从棋盘格图像中估计相机内参（焦距、畸变系数），本质是非线性最小二乘优化。
SLAM 后端优化：机器人同时定位与建图（SLAM）的后端，用最小二乘（图优化 / Bundle Adjustment）对所有位姿和路标点做全局一致性优化。
传感器融合 / 卡尔曼滤波：卡尔曼滤波的更新步骤，数学本质是最小二乘——融合多个带噪声的传感器读数得到最优估计。

典型类库与工具

类库	语言	说明
numpy.linalg.lstsq	Python	最基础的线性最小二乘求解器，底层用 SVD，数值稳定
scipy.optimize.curve_fit	Python	非线性曲线拟合，基于 Levenberg-Marquardt，接口简洁
scipy.optimize.least_squares	Python	更通用的非线性最小二乘，支持边界约束和鲁棒损失函数
statsmodels	Python	提供 OLS / WLS / 鲁棒回归，附带完整的统计检验（R²、p 值、置信区间）
torch.linalg.lstsq	Python	PyTorch 的可微最小二乘求解器，支持梯度穿过 lstsq 层反向传播
Ceres Solver	C++	Google 开源的非线性最小二乘库，SLAM / 计算机视觉的标准后端
GTSAM	C++	Georgia Tech 因子图优化库，SLAM 核心工具，支持增量求解

术语表

术语	英文	解释
残差	Residual	观测值与模型预测值之差，最小二乘最小化残差平方和
正规方程	Normal Equation	线性最小二乘的解析解公式 β̂ = (XᵀX)⁻¹ Xᵀy
设计矩阵	Design Matrix	由特征构成的矩阵 X，每行一个样本、每列一个特征
Gauss-Newton	Gauss-Newton Method	非线性最小二乘的经典迭代法，用一阶展开近似避免完整 Hessian
Levenberg-Marquardt	LM Method	GN 与梯度下降的自适应混合，曲线拟合的工业标准算法
鲁棒回归	Robust Regression	用对异常值不敏感的损失（如 Huber）代替平方损失，抵抗离群点
RANSAC	RANSAC	随机抽样一致算法，从含大量离群点的数据中鲁棒地拟合模型
加权最小二乘	WLS	给不同样本赋予不同权重（如按测量精度），处理异方差数据
岭回归	Ridge Regression	加 L2 正则化的最小二乘，解决多重共线性导致的不稳定
雅可比矩阵	Jacobian Matrix	非线性函数的一阶导数矩阵，Gauss-Newton 迭代中用于线性化残差
条件数	Condition Number	衡量矩阵接近奇异程度，越大则求逆越不稳定
可微最小二乘	Differentiable Least Squares	嵌入深度学习框架的最小二乘，梯度可穿过求解层反向传播
随机化 SVD	Randomized SVD	用随机投影加速的大规模矩阵低秩近似分解方法