当前位置：首页 > news >正文

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

news 2025/10/29 7:59:05

1. 引言

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

2. 背景

梯度下降法在人工智能的机器学习中使用的非常多，因为机器学习的训练过程通常被形式化为经验风险最小化问题（Empirical Risk Minimization, ERM）：即在训练数据上最小化损失函数。而最小二乘问题其实也是经验风险最小化问题的一种，甚至机器学习的某些任务（比如回归）本身就是最小二乘问题。经验风险最小化问题是一种通用的函数拟合框架，不过损失函数有所不同，通常使用梯度下降法来进行求解。

那么为什么机器学习中使用梯度下降法来求解，而计算机视觉（SLAM、SfM、相机标定、BA）中使用Gauss-Newton/Levenberg-Marquardt来进行求解呢？这是因为机器学习的问题可以只用关心局部的“好解”，而不用像计算机视觉问题那样需要求解全局的“精确最小值”；另外，机器学习问题规模巨大、结构复杂，使用梯度下降法要简单、健壮、高效的多。

3. 求解

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

$min⁡θS(θ)=r(θ)2=∑i=1mri(θ)2\min_{\theta} S(\theta) = \ \mathbf{r}(\theta)\ ^2 = \sum_{i=1}^m r_i(\theta)^2$

其中：
$θ∈Rn\theta \in \mathbb{R}^n$ ：待优化的参数向量（比如曲线的系数）
$r(θ)=[r1(θ)⋮rm(θ)]\mathbf{r}(\theta) = \begin{bmatrix} r_1(\theta) \\ \vdots \\ r_m(\theta) \end{bmatrix}$ ：残差向量， $ri(θ)=yi−f(xi;θ)r_i(\theta) = y_i - f(x_i; \theta)$
$S(θ)S(\theta)$ ：目标函数（损失函数），是我们要最小化的残差平方和

梯度下降法的核心思想是：在当前点，沿着目标函数下降最快的方向走一步，然后重复。而这个“最快下降方向”就是负梯度方向 $−∇S(θ)-\nabla S(\theta)$ 。因此问题的关键在于计算目标函数 $S(θ)=r(θ)2S(\theta) = \ \mathbf{r}(\theta)\ ^2$ 的梯度。根据求导的链式法则：

$∇S(θ)=∂∂θ(r(θ)Tr(θ))=2J(θ)Tr(θ)\nabla S(\theta) = \frac{\partial}{\partial \theta} \left( \mathbf{r}(\theta)^T \mathbf{r}(\theta) \right) = 2 \, J(\theta)^T \mathbf{r}(\theta)$

其中：
$J(θ)J(\theta)$ ：雅可比矩阵（Jacobian），大小为 $\times n$
$J(θ)=[∂r1∂θ1⋯∂r1∂θn⋮⋱⋮∂rm∂θ1⋯∂rm∂θn]=∂r∂θTJ(\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}$
即目标函数的梯度是： $∇S(θ)=2J(θ)Tr(θ)\nabla S(\theta) = 2 J(\theta)^T \mathbf{r}(\theta)$ 。

另一方面，在每次梯度下降之后，需要更新参数向量：

$θk+1=θk−α⋅∇S(θk)=θk−2α⋅JkTrk\boxed{ \theta_{k+1} = \theta_k - \alpha \cdot \nabla S(\theta_k) = \theta_k - 2\alpha \cdot J_k^T \mathbf{r}_k }$

其中：
$θk\theta_k$ ：第 $k$ 次迭代的参数
$α>0\alpha > 0$ ：学习率（step size），控制步长
$Jk=J(θk)J_k = J(\theta_k)$ ， $rk=r(θ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实现梯度下降法求解非线性最小二乘问题的例子。例子中模型函数为 $\boldsymbol{\theta}) = a e ^{bx}$ ：

#include <Eigen/Dense>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>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<double>& x_data,const vector<double>& y_data, const Vector2d& theta) {int N = x_data.size();VectorXd r(N);for (int i = 0; i < N; ++i) {r(i) = y_data[i] - model(x_data[i], theta);}return r;
}// 计算 Jacobian 矩阵 (N x 2): ∂r_i/∂a, ∂r_i/∂b
MatrixXd computeJacobian(const vector<double>& 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[i];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 << 2.0, -0.3;  // a=2.0, b=-0.3 → y = 2 * exp(-0.3 * x)cout << "真实参数: a = " << true_params(0) << ", b = " << true_params(1)<< endl;// ========================// 2. 生成带噪声的数据// ========================int N = 20;vector<double> x_data(N), y_data(N);random_device rd;mt19937 gen(rd());normal_distribution<double> noise(0.0, 0.05);  // 小噪声for (int i = 0; i < N; ++i) {x_data[i] = -2.0 + i * 0.4;  // x 从 -2 到 6double y_true = model(x_data[i], true_params);y_data[i] = y_true + noise(gen);}// ========================// 3. 初始化参数// ========================Vector2d theta;theta << 1.0, 0.0;  // 初始猜测: a=1.0, b=0.0cout << "初始猜测: a = " << theta(0) << ", b = " << theta(1) << endl;// ========================// 4. 梯度下降法// ========================int max_iter = 500;double alpha = 5e-3;  // 学习率double tol = 1e-6;cout << "\n开始梯度下降...\n";cout << "迭代\t残差平方和\t\t参数 a\t\t参数 b\n";cout << "----\t----------\t\t------\t\t------\n";for (int iter = 0; iter < max_iter; ++iter) {// 计算残差VectorXd r = computeResiduals(x_data, y_data, theta);double cost = r.squaredNorm();// 计算梯度MatrixXd J = computeJacobian(x_data, theta);Vector2d gradient = 2.0 * J.transpose() * r;// 打印当前状态（每10次）if (iter % 10 == 0) {cout << iter << "\t" << cost << "\t\t" << theta(0) << "\t\t" << theta(1)<< endl;}// 终止条件if (gradient.norm() < tol) {cout << "收敛！梯度范数: " << gradient.norm() << endl;break;}// 更新参数theta -= alpha * gradient;}// ========================// 5. 输出结果// ========================cout << "\n--- 拟合完成 ---" << endl;cout << "估计参数: a = " << theta(0) << ", b = " << theta(1) << endl;cout << "真实参数: a = " << true_params(0) << ", b = " << true_params(1)<< endl;return 0;
}

运行结果如下：

真实参数: a = 2, b = -0.3
初始猜测: a = 1, b = 0开始梯度下降...
迭代    残差平方和              参数 a          参数 b
----    ----------              ------          ------
0       22.7591         1               0
10      1.11435         1.72284         -0.345
20      0.100641                1.93634         -0.301778
30      0.0326195               1.99193         -0.294493
40      0.0286004               2.00545         -0.292882
50      0.0283681               2.0087          -0.292503
60      0.0283548               2.00948         -0.292413
70      0.028354                2.00967         -0.292391
80      0.0283539               2.00971         -0.292386
90      0.0283539               2.00972         -0.292385
100     0.0283539               2.00972         -0.292384
110     0.0283539               2.00973         -0.292384
120     0.0283539               2.00973         -0.292384
收敛！梯度范数: 9.36104e-07--- 拟合完成 ---
估计参数: a = 2.00973, b = -0.292384
真实参数: a = 2, b = -0.3

求解的关键还是在于计算雅可比矩阵，对于问题模型函数 $\boldsymbol{\theta}) = a e ^{bx}$ 来说，雅可比矩阵应该是：

$J(θ)=[∂(y1−aebx1)∂a∂(y1−aebx1)∂b∂(y2−aebx2)∂a∂(y2−aebx2)∂b⋮⋮∂(ym−aebxm)∂a∂(ym−aebxm)∂b]=−[ebx1aebx1x1ebx2aebx2x2⋮⋮ebxmaebxmxm]J(\theta) = \begin{bmatrix} \frac{\partial (y_1-ae^{bx_1})}{\partial a} & \frac{\partial (y_1-ae^{bx_1})}{\partial b} \\ \frac{\partial (y_2-ae^{bx_2})}{\partial a} & \frac{\partial (y_2-ae^{bx_2})}{\partial b} \\ \vdots & \vdots \\ \frac{\partial (y_m-ae^{bx_m})}{\partial a} & \frac{\partial (y_m-ae^{bx_m})}{\partial b} \\ \end{bmatrix}= -\begin{bmatrix} e^{bx_1} & ae^{bx_1}x_1 \\ e^{bx_2} & ae^{bx_2}x_2 \\ \vdots & \vdots \\ e^{bx_m} & ae^{bx_m}x_m \\ \end{bmatrix}$

对比代码中的实现：

// 计算 Jacobian 矩阵 (N x 2): ∂r_i/∂a, ∂r_i/∂b
MatrixXd computeJacobian(const vector<double>& 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[i];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;
}

另外，除了迭代过程中的初始条件和迭代停止条件，控制步长的学习率也需要注意。设置的学习率过小，迭代次数就会很长导致收敛很慢；而设置的学习率过大，就容易略过最优解导致结果不问题。因此，在实际的工程应用中，通常不会使用原始的梯度下降法，而是根据需求使用不同优化版本的梯度下降法。关于这一点，有机会的话会在后续的文章中进一步论述。

查看全文

http://www.dtcms.com/a/540569.html