最小二乘问题详解7：正则化最小二乘

1. 引言

在之前的文章《最小二乘问题详解4：非线性最小二乘》、《最小二乘问题详解5：非线性最小二乘求解实例》和《最小二乘问题详解6：梯度下降法》中分别介绍了使用Gauss-Newton方法（简称GN方法）和梯度下降法求解最小二乘问题之后，让我们插入另一个基础知识：正则化最小二乘（Regularized Least Squares），也就是大家常说的岭估计（Ridge Estimator），因为接下来要介绍的 Levenberg-Marquardt方法会用到这个思想。

本文讨论的岭估计在机器学习中也常称为岭回归（Ridge Regression）

2. 问题

复习《最小二乘问题详解2：线性最小二乘求解》中讨论的标准线性最小二乘问题：

$$\underset{\theta }{\min }\Vert A\theta -b{\Vert }^2$$

其解为正规方程 $A^{T}A\theta =A^{T}b$的解。当$A$列满秩时，解唯一且为：

$${\theta }^{\ast }=(A^{T}A{)}^{-1}{A}^{T}b$$

然而，在实际应用中，直接使用这个解可能会遇到一些问题。

1. 矩阵病态（Ill-conditioning）：

   • 矩阵 $A^{T}A$ 的条件数可能非常大。
   • 这意味着即使数据 $b$ 中有微小的扰动，也会导致解 ${\theta }^{\ast }$ 发生剧烈变化，数值不稳定。
   • 条件数大的一个常见原因是特征之间高度相关（多重共线性），此时 $A^{T}A$ 接近奇异，逆矩阵难以精确计算。

2. 过拟合（Overfitting）：

   • 当模型参数过多或特征维度很高时，标准最小二乘倾向于拟合训练数据中的噪声，导致泛化能力差。
   • 解 ${\theta }^{\ast }$ 的模长（范数）可能非常大，表示模型对某些特征赋予了不合理的高权重。

3. 秩亏（Rank Deficiency）：

   • 若 $A$ 不是列满秩（即 $\text{rank}(A)<n$），则 $A^{T}A$ 是奇异的，无法求逆，正规方程有无穷多解。

矩阵的病态问题详见6.1节。

为了解决这些问题，我们引入正则化（Regularization）——通过在目标函数中加入一个额外的惩罚项，来“约束”解的行为，使其更稳定、更平滑，并避免过拟合。其实，正则化是一个通用的数学和机器学习思想：在优化问题中，通过向目标函数添加一个额外的“惩罚项”（penalty term），来引导解具有某种期望的性质，比如平滑性、稀疏性、小范数、结构简单等，从而提高模型的稳定性、泛化能力或可解释性。

3. 定义

最常用的正则化方法之一是Tikhonov 正则化，其核心思想是在原始的残差平方和基础上，加上一个关于参数$\theta$的L2范数平方的惩罚项：

$$\boxed{\underset{\theta }{\min }\left(\Vert A\theta -b{\Vert }^2+\lambda \Vert \theta {\Vert }^2\right)}$$

其中：

• $\Vert A\theta -b{\Vert }^2$是数据拟合项（data fidelity term），衡量模型对数据的拟合程度。
• $\Vert \theta {\Vert }^2={\theta }^T\theta$是正则化项（regularization term），也叫 Tikhonov 项，它惩罚较大的参数值。
• $\lambda >0$是正则化参数（regularization parameter），控制正则化的强度：
  • $\lambda \to 0$：正则化作用消失，退化为标准最小二乘。
  • $\lambda \to \infty$：正则化主导，迫使 $\theta \to 0$，模型趋于常数零。

直观理解就是：通过加入正则项，不仅要保证预测误差小，还要保证模型参数不要太大。这相当于在“拟合数据”和“保持模型简单”之间做权衡。

4. 求解

将目标函数展开：

$$f(\theta )=(A\theta -b{)}^T(A\theta -b)+\lambda {\theta }^T\theta ={\theta }^T{A}^{T}A\theta -2b^{T}A\theta +b^{T}b+\lambda {\theta }^T\theta$$

对$\theta$求梯度并令其为零：

$${\nabla }_{\theta }f(\theta )=2A^{T}A\theta -2A^{T}b+2\lambda \theta =0$$

整理得：

$$(A^{T}A+\lambda I)\theta =A^{T}b$$

这就是岭估计的正规方程。

由于$\lambda >0$，矩阵$A^{T}A$是半正定的，而$\lambda I$是正定的，因此$A^{T}A+\lambda I$是严格正定的，从而总是可逆的，无论$A$是否列满秩！于是，岭估计的闭式解为：

$$\boxed{{\theta }_{\text{ridge}}^{\ast }=(A^{T}A+\lambda I{)}^{-1}{A}^{T}b}$$

可见，岭估计只是在$A^{T}A$的对角线上加了一个小量$\lambda$，这相当于“抬升”了所有特征值，使得原本接近零的特征值变得远离零，从而显著改善了矩阵的条件数，提高了数值稳定性。

关于正定矩阵的问题详见6.2节。

从几何角度看，标准最小二乘是寻找使 $A\theta$ 投影最接近 $b$ 的点；而岭估计则在此基础上“收缩”参数空间，使 $\theta$ 向零靠近。这可以理解为给 $\theta$ 的空间加上一个“弹簧”，防止参数跑得太远。

5. 分解

虽然有了闭式解，但在实际数值计算中，仍然不推荐直接计算$(A^{T}A+\lambda I{)}^{-1}$，因为$A^{T}A$可能病态，显式构造$A^{T}A$会损失精度。

5.1 QR分解

将正则化最小二乘问题转化为一个更大的最小二乘问题：

$$\underset{\theta }{\min }{\Vert \left[\begin{array}{c}A\\ \sqrt{\lambda }I\end{array}\right]\theta -\left[\begin{array}{c}b\\ 0\end{array}\right]\Vert }^2=\Vert A\theta -b{\Vert }^2+\lambda \Vert \theta {\Vert }^2$$

这是一个标准的最小二乘问题，可以用 QR 分解高效求解。

令：

$$\tilde{A}=\left[\begin{array}{c}A\\ \sqrt{\lambda }I\end{array}\right],\tilde{b}=\left[\begin{array}{c}b\\ 0\end{array}\right]$$

对$\tilde{A}$做QR分解：$\tilde{A}=\tilde{Q}\tilde{R}$

则解为：

$${\theta }_{\text{ridge}}^{\ast }={\tilde{R}}^{-1}{\tilde{Q}}_1^{T}b$$

其中${\tilde{Q}}_1$ 是 $\tilde{Q}$的前 $m$ 行（对应 $A$ 部分）。

5.2 SVD分解

设 $A=U\Sigma V^T$ 是 $A$ 的 SVD。

代入岭估计的闭式解：

$${\theta }_{\text{ridge}}^{\ast }=(V{\Sigma }^T{U}^{T}U\Sigma V^T+\lambda I{)}^{-1}V{\Sigma }^T{U}^{T}b=(V({\Sigma }^T\Sigma +\lambda I)V^T{)}^{-1}V{\Sigma }^T{U}^{T}b$$

利用 $V^{T}V=I$，得：

$${\theta }_{\text{ridge}}^{\ast }=V({\Sigma }^T\Sigma +\lambda I{)}^{-1}{V}^{T}V{\Sigma }^T{U}^{T}b=V({\Sigma }^T\Sigma +\lambda I{)}^{-1}{\Sigma }^T{U}^{T}b$$

记 ${\Sigma }^T\Sigma =\text{diag}({\sigma }_1^2,\dots ,{\sigma }_n^2)$，则：

$${\theta }_{\text{ridge}}^{\ast }=V\left[\begin{array}{ccc}\frac{{\sigma }_1}{{\sigma }_1^2+\lambda }& & \\ & \ddots & \\ & & \frac{{\sigma }_n}{{\sigma }_n^2+\lambda }\end{array}\right]U^{T}b$$

即：

$$\boxed{{\theta }_{\text{ridge}}^{\ast }=\sum _{i=1}^r\frac{{\sigma }_i}{{\sigma }_i^2+\lambda }(u_i^{T}b)v_i}$$

对比标准最小二乘解（SVD 形式）：

$${\theta }^{\ast }=\sum _{i=1}^r\frac{1}{{\sigma }_i}(u_i^{T}b)v_i$$

可见，岭估计通过因子 $\frac{{\sigma }_i}{{\sigma }_i^2+\lambda }$ 对小奇异值方向进行了压制。当 ${\sigma }_i\to 0$ 时，标准解会爆炸（$1/{\sigma }_i\to \infty$），但岭估计中该项趋于 0，从而避免了对噪声方向的过度放大。

6. 补充

有一些《线性代数》、《矩阵论》相关的知识笔者也忘记了，这里就总结一下。如果有的读者很熟悉，可以直接略过。

6.1 病态矩阵

病态矩阵指的是矩阵条件数（Condition Number）很大的矩阵。矩阵病态会导致输入数据的微小扰动（如 $b$ 或 $A$ 的舍入误差），线性方程组 $A\theta =b$ 解就会剧烈变化。换句话说，系统对噪声极度敏感，数值计算中结果不可靠。

对于一个可逆矩阵 $A\in {\mathbb{R}}^{n\times n}$，其谱条件数（Spectral Condition Number）定义为：

$$\kappa (A)=\frac{{\sigma }_{\max }(A)}{{\sigma }_{\min }(A)}$$

其中：
${\sigma }_{\max }(A)$ 是 $A$ 的最大奇异值，
${\sigma }_{\min }(A)$ 是 $A$ 的最小奇异值。

注意这个定义也适用于非方阵 $A\in {\mathbb{R}}^{m\times n}$（$m\ge n$），只要 $A$ 列满秩。条件数 $\kappa (A)$ 的含义如下：

• $\kappa (A)\approx 1$ 良态，解非常稳定
• $\kappa (A)<10^2$ 良态
• $10^2\le \kappa (A)<10^6$ 中等病态
• $\kappa (A)\ge 10^6$ 严重病态，浮点计算可能完全失效
• $\kappa (A)>10^{14}$ 双精度浮点数（约16位有效数字）完全失效

任何实矩阵 $A\in {\mathbb{R}}^{m\times n}$ 都可以进行奇异值分解（SVD）：

$$A=U\Sigma V^T$$

其中：
$U\in {\mathbb{R}}^{m\times m}$ 是左奇异向量矩阵（正交），
$V\in {\mathbb{R}}^{n\times n}$ 是右奇异向量矩阵（正交），
$\Sigma \in {\mathbb{R}}^{m\times n}$ 是对角矩阵，对角线元素为奇异值(Singular Value) ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r\ge 0$，$r=\min (m,n)$。

奇异值表示矩阵 $A$ 在各个正交方向上的“拉伸”程度。如果 ${\sigma }_{\min }\to 0$，说明 $A$ 把某个方向“压扁”了，接近奇异，也就是矩阵病态。

6.2 正定矩阵

设 $A\in {\mathbb{R}}^{n\times n}$ 是一个实对称矩阵（即 $A=A^T$），我们有如下定义：

如果对所有非零向量 $x\in {\mathbb{R}}^n,x\ne 0$，都有：

$$x^{T}Ax>0$$

则称 $A$ 是正定矩阵。

如果对所有非零向量 $x\in {\mathbb{R}}^n,x\ne 0$，都有：

$$x^{T}Ax\ge 0$$

则称 $A$ 是半正定矩阵。

对于实对称矩阵 $A$，我们可以进行谱分解（特征值分解）：

$$A=Q\Lambda Q^T$$

其中 $Q$ 是正交矩阵（列是特征向量），$\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_n)$ 是对角矩阵，对角元素是特征值。

那么：

• $A$ 正定 $⟺$ 所有特征值 ${\lambda }_i>0$
• $A$ 半正定 $⟺$ 所有特征值 ${\lambda }_i\ge 0$

既然正定矩阵的所有特征值都 大于零，那么意味着正定矩阵是满秩，也就是正定矩阵一定可逆。

回到岭估计中的 $A^{T}A+\lambda I$：

1. $A^{T}A$ 是半正定的：

   • 因为对任意 $x$，有 $x^T{A}^{T}Ax=\Vert Ax{\Vert }^2\ge 0$
   • 所以它的所有特征值 $\ge 0$

2. 加上 $\lambda I$（其中 $\lambda >0$）后：

   • 新矩阵是 $A^{T}A+\lambda I$
   • 它的特征值变为：原来的每个特征值 $\lambda_i $ 变成 ${\lambda }_i+\lambda$
   • 原来最小可能是 0，现在最小是 $\lambda >0$
   • 所以所有特征值 $>0$ ⇒ 正定 ⇒ 可逆

因此，$(A^{T}A+\lambda I{)}^{-1}$ 总是存在的，无论 $A$ 是否列满秩。这正是岭估计的精髓所在：通过加一个小的正数 $\lambda$ 到对角线上，把原本可能奇异（不可逆）的 $A^{T}A$ “修复”成一个严格正定、可逆的矩阵，从而保证了解的唯一性和数值稳定性。

7. 实例

如果线性最小二乘问题的设计矩阵$A$接近线性相关，那么普通方法求得的解不稳定，可以使用岭估计来给出稳定解。代码实现如下：

#include <Eigen/Dense>
#include <iostream>
#include <random>
#include <vector>using namespace Eigen;
using namespace std;int main() {// 设置随机数生成器default_random_engine gen;normal_distribution<double> noise(0.0, 0.5);  // 噪声 ~ N(0, 0.5)// 数据生成：y = x1 + x2 + 0.1*x3 + noise// 但 x3 ≈ x1 + x2 （高度相关，造成病态）int n_samples = 20;int n_features = 3;MatrixXd A(n_samples, n_features);VectorXd b(n_samples);for (int i = 0; i < n_samples; ++i) {double x1 = (double)rand() / RAND_MAX * 10;double x2 = (double)rand() / RAND_MAX * 10;double x3 = x1 + x2 + (double)rand() / RAND_MAX * 0.1;  // x3 ≈ x1 + x2A(i, 0) = x1;A(i, 1) = x2;A(i, 2) = x3;double true_y = 1.0 * x1 + 1.0 * x2 + 0.1 * x3;b(i) = true_y + noise(gen);  // 加噪声}//使用 SVD 计算条件数BDCSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);cout << "Condition number of A: "<< svd.singularValues()(0) / svd.singularValues()(2)<< endl;  // 假设 3 个奇异值// 普通最小二乘解：theta = (A^T A)^{-1} A^T bVectorXd theta_ols = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);// 岭回归解：theta = (A^T A + λI)^{-1} A^T bMatrixXd AtA = A.transpose() * A;VectorXd Atb = A.transpose() * b;double lambda = 0.01;MatrixXd ridge_matrix =AtA + lambda * MatrixXd::Identity(n_features, n_features);VectorXd theta_ridge = ridge_matrix.ldlt().solve(Atb);// 输出结果cout << "True weights: [1.0, 1.0, 0.1]" << endl;cout << "OLS weights:  " << theta_ols.transpose() << endl;cout << "Ridge weights:" << theta_ridge.transpose() << endl;return 0;
}

这里可以看到，我们使第三个已知参数向量，可以用第一个已知参数向量和第二个已知参数向量接近线性表示（$x3\approx x1+x2$），那么设计矩阵就接近线性相关，这个设计矩阵就是病态的。分别用普通最小二乘求解和岭估计来求解，最后结果如下：

Condition number of A: 715.286
True weights: [1.0, 1.0, 0.1]
OLS weights:    2.00569   2.03977 -0.938346
Ridge weights: 1.02399  1.05935 0.038262

理论上来说，使用QR/SVD方法可以一定程度上解决矩阵的病态问题。但是从这个实例结果可以看到，即使使用最稳健的 SVD 方法求解普通最小二乘，参数估计仍严重偏离真实值。这是因为问题本身的结构特征高度相关，导致参数不可识别。

相比之下，岭估计通过引入正则化项$\lambda I$，显著改善了矩阵的条件数，给出了稳定且接近真实的参数估计。岭估计通过牺牲一定的精度，提升模型的泛化能力，保证在噪声存在下仍能稳定预测。岭估计的另外一个问题是正则化参数$\lambda$不太好进行给定，往往需要一些先验经验的辅助才能确定，有机会再进一步进行讨论了。