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

1. 引言

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

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

2. 问题

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

\[\min_{\theta} \|A\theta - b\|^2\]
其解为正规方程 \(A^T A \theta = A^T b\)的解。当\(A\)列满秩时，解唯一且为：

\[\theta^* = (A^T A)^{-1} A^T b\]
然而，在实际应用中，直接使用这个解可能会遇到一些问题。

• 矩阵病态（Ill-conditioning）：
  
  
  
  • 矩阵 \(A^T A\) 的条件数可能非常大。
    
  • 这意味着即使数据 \(b\) 中有微小的扰动，也会导致解 \(\theta^*\) 发生剧烈变化，数值不稳定。
    
  • 条件数大的一个常见原因是特征之间高度相关（多重共线性），此时 \(A^T A\) 接近奇异，逆矩阵难以精确计算。
    
  
  
• 过拟合（Overfitting）：
  
  
  
  • 当模型参数过多或特征维度很高时，标准最小二乘倾向于拟合训练数据中的噪声，导致泛化能力差。
    
  • 解 \(\theta^*\) 的模长（范数）可能非常大，表示模型对某些特征赋予了不合理的高权重。
    
  
  
• 秩亏（Rank Deficiency）：
  
  
  
  • 若 \(A\) 不是列满秩（即 \(\text{rank}(A) < n\)），则 \(A^T A\) 是奇异的，无法求逆，正规方程有无穷多解。
    
  
  

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

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

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

\[\boxed{\min_{\theta} \left( \|A\theta - b\|^2 + \lambda \|\theta\|^2 \right)}\]
其中：

• \(\|A\theta - 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) = (A\theta - b)^T (A\theta - 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节。

从几何角度看，标准最小二乘是寻找使 \(A\theta\) 投影最接近 \(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=\|A\theta - 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\) 的舍入误差），线性方程组 \(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 \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\)：

• \(A^T A\) 是半正定的：
  
  
  
  • 因为对任意 \(x\)，有 \(x^T A^T A x = \|Ax\|^2 \ge 0\)
    
  • 所以它的所有特征值 \(\ge 0\)
    
  
  
• 加上 \(\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\)接近线性相关，那么普通方法求得的解不稳定，可以使用岭估计来给出稳定解。代码实现如下：
[code]#include #include #include #include using namespace Eigen;using namespace std;int main() {  // 设置随机数生成器  default_random_engine gen;  normal_distribution 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 svd(A, ComputeThinU | ComputeThinV);  cout