最小二乘问题详解6：梯度下降法

1. 引言

在之前的两篇文章《最小二乘问题详解4：非线性最小二乘》、《最小二乘问题详解5：非线性最小二乘求解实例》中，笔者介绍了非线性最小二乘问题，并使用Gauss-Newton方法来进行求解。不过，求解非线性最小二乘问题还有另外一种方法——梯度下降法。
2. 背景

梯度下降法在人工智能的机器学习中使用的非常多，因为机器学习的训练过程通常被形式化为经验风险最小化问题（Empirical Risk Minimization, ERM）：即在训练数据上最小化损失函数。而最小二乘问题其实也是经验风险最小化问题的一种，甚至机器学习的某些任务（比如回归）本身就是最小二乘问题。经验风险最小化问题是一种通用的函数拟合框架，不过损失函数有所不同，通常使用梯度下降法来进行求解。
那么为什么机器学习中使用梯度下降法来求解，而计算机视觉（SLAM、SfM、相机标定、BA）中使用Gauss-Newton/Levenberg-Marquardt来进行求解呢？这是因为机器学习的问题可以只用关心局部的“好解”，而不用像计算机视觉问题那样需要求解全局的“精确最小值”；另外，机器学习问题规模巨大、结构复杂，使用梯度下降法要简单、健壮、高效的多。
3. 求解

接下来就来介绍一下使用梯度下降法求解非线性最小二乘问题。还是先看非线性最小二乘问题的定义：

\[\min_{\theta} S(\theta) = \ \mathbf{r}(\theta)\ ^2 = \sum_{i=1}^m r_i(\theta)^2 \]
其中：
\(\theta \in \mathbb{R}^n\)：待优化的参数向量（比如曲线的系数）
\(\mathbf{r}(\theta) = \begin{bmatrix} r_1(\theta) \\ \vdots \\ r_m(\theta) \end{bmatrix}\)：残差向量，\(r_i(\theta) = y_i - f(x_i; \theta)\)
\(S(\theta)\)：目标函数（损失函数），是我们要最小化的残差平方和
梯度下降法的核心思想是：在当前点，沿着目标函数下降最快的方向走一步，然后重复。而这个“最快下降方向”就是负梯度方向\(-\nabla S(\theta)\)。因此问题的关键在于计算目标函数\(S(\theta) = \ \mathbf{r}(\theta)\ ^2\)的梯度。根据求导的链式法则：

\[\nabla S(\theta) = \frac{\partial}{\partial \theta} \left( \mathbf{r}(\theta)^T \mathbf{r}(\theta) \right) = 2 \, J(\theta)^T \mathbf{r}(\theta) \]
其中：
\(J(\theta)\)：雅可比矩阵（Jacobian），大小为 \(m \times n\)

\[J(\theta) = \begin{bmatrix} \frac{\partial r_1}{\partial \theta_1} & \cdots & \frac{\partial r_1}{\partial \theta_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial r_m}{\partial \theta_1} & \cdots & \frac{\partial r_m}{\partial \theta_n} \end{bmatrix} = \frac{\partial \mathbf{r}}{\partial \theta^T} \]
即目标函数的梯度是：\(\nabla S(\theta) = 2 J(\theta)^T \mathbf{r}(\theta)\)。
另一方面，在每次梯度下降之后，需要更新参数向量：

\[\boxed{ \theta_{k+1} = \theta_k - \alpha \cdot \nabla S(\theta_k) = \theta_k - 2\alpha \cdot J_k^T \mathbf{r}_k } \]
其中：
\(\theta_k\)：第 \(k\) 次迭代的参数
\(\alpha > 0\)：学习率（step size），控制步长
\(J_k = J(\theta_k)\)，\(\mathbf{r}_k = \mathbf{r}(\theta_k)\)
因此，将梯度下降方法完整的流程总结如下：

• 初始化：选一个初始猜测 θ₀
  
• 设置学习率 α（例如 0.01）
  
• 对 k = 0, 1, 2, ... 直到收敛：
  a. 计算残差：\(r_k = y - f(x; θ_k)\)
  b. 计算雅可比矩阵：\(J_k = J(θ_k)\)
  c. 计算梯度：\(g_k = 2 J_k^T r_k\)
  d. 更新参数：\(θ_{k+1} = θ_k - α g_k\)
  e. 检查是否收敛：\(Δθ = θ_{k+1} - θ_k < ε\)或\(g_k < ε\)或\(S(θ)\)变化很小
  
• 输出最终参数 θ
  

4. 实例

从上述求解过程可以看到，梯度下降法其实比之前文章中介绍的Gauss-Newton方法要简单很多，那么这里还是给出一个只使用Eigen实现梯度下降法求解非线性最小二乘问题的例子。例子中模型函数为\(f(x; \boldsymbol{\theta}) = a e ^{bx}\)：
[code]#include  #include  #include  #include  #include   using namespace std; using namespace Eigen;  // 模型函数: y = a * exp(b * x) double model(double x, const Vector2d& theta) {   double a = theta(0);   double b = theta(1);   return a * exp(b * x); }  // 计算残差: r_i = y_i - f(x_i; a, b) VectorXd computeResiduals(const vector& x_data,                           const vector& y_data, const Vector2d& theta) {   int N = x_data.size();   VectorXd r(N);   for (int i = 0; i < N; ++i) {     r(i) = y_data - model(x_data, theta);   }   return r; }  // 计算 Jacobian 矩阵 (N x 2): ∂r_i/∂a, ∂r_i/∂b MatrixXd computeJacobian(const vector& x_data, const Vector2d& theta) {   int N = x_data.size();   MatrixXd J(N, 2);   double a = theta(0);   double b = theta(1);    for (int i = 0; i < N; ++i) {     double x = x_data;     double exp_bx = exp(b * x);  // exp(b*x)      J(i, 0) = -exp_bx;          // ∂r/∂a = -exp(b*x)     J(i, 1) = -a * exp_bx * x;  // ∂r/∂b = -a * exp(b*x) * x   }   return J; }  int main() {   // ========================   // 1. 真实参数   // ========================   Vector2d true_params;   true_params