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

1. 引言

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

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

2. 问题

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

[min_{theta} |Atheta - b|^2 ]

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

[theta^* = (A^T A)^{-1} A^T b ]

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

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

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

2. 过拟合（Overfitting）：

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

3. 秩亏（Rank Deficiency）：

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

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

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

3. 定义

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

[boxed{ min_{theta} left( |Atheta - b|^2 + lambda |theta|^2 right) } ]

其中：

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

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

4. 求解

将目标函数展开：

[f(theta) = (Atheta - b)^T (Atheta - b) + lambda theta^T theta = theta^T A^T A theta - 2 b^T A theta + b^T b + lambda theta^T theta ]

对(theta)求梯度并令其为零：

[nabla_theta f(theta) = 2 A^T A theta - 2 A^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}} = (A^T A + lambda I)^{-1} A^T b } ]

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

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

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

5. 分解

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

5.1 QR分解

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

[min_{theta} left| begin{bmatrix} A \ sqrt{lambda} I end{bmatrix} theta - begin{bmatrix} b \ 0 end{bmatrix} right|^2= |Atheta - b|^2 + lambda |theta|^2 ]

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

令：

[tilde{A} = begin{bmatrix} A \ sqrt{lambda} I end{bmatrix},quad tilde{b} = begin{bmatrix} b \ 0 end{bmatrix} ]

对(tilde{A})做QR分解：(tilde{A} = tilde{Q} tilde{R})

则解为：

[theta^*_{text{ridge}} = 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}} = (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}} = 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}} = V begin{bmatrix} frac{sigma_1}{sigma_1^2 + lambda} & & \ & ddots & \ & & frac{sigma_n}{sigma_n^2 + lambda} end{bmatrix} U^T b ]

即：

[boxed{ theta^*_{text{ridge}} = sum_{i=1}^r frac{sigma_i}{sigma_i^2 + lambda} (u_i^T b) v_i } ]

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

[theta^* = 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) 的舍入误差），线性方程组 (Atheta = 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 geq n)），只要 (A) 列满秩。条件数 (kappa(A)) 的含义如下：

• (kappa(A) approx 1) 良态，解非常稳定
• (kappa(A) < 10^2) 良态
• (10^2 leq kappa(A) < 10^6) 中等病态
• (kappa(A) geq 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 geq sigma_2 geq cdots geq sigma_r geq 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 A x > 0 ]

则称 (A) 是正定矩阵。

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

[x^T A x ge 0 ]

则称 (A) 是半正定矩阵。

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

[A = Q Lambda Q^T ]

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

那么：

• (A) 正定 (Longleftrightarrow) 所有特征值 (lambda_i > 0)
• (A) 半正定 (Longleftrightarrow) 所有特征值 (lambda_i ge 0)

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

回到岭估计中的 (A^T A + lambda I)：

1. (A^T A) 是半正定的：

   • 因为对任意 (x)，有 (x^T A^T A x = |Ax|^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 + x2      A(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 b   VectorXd theta_ols = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);    // 岭回归解：theta = (A^T A + λI)^{-1} A^T b   MatrixXd 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 ≈ 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 方法求解普通最小二乘，参数估计仍严重偏离真实值。这是因为问题本身的结构特征高度相关，导致参数不可识别。

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