1. 引言
在之前的两篇文章《最小二乘问题详解4：非线性最小二乘》、《最小二乘问题详解5：非线性最小二乘求解实例》中，笔者介绍了非线性最小二乘问题，并使用Gauss-Newton方法来进行求解。不过，求解非线性最小二乘问题还有另外一种方法——梯度下降法。
2. 背景
梯度下降法在人工智能的机器学习中使用的非常多，因为机器学习的训练过程通常被形式化为经验风险最小化问题（Empirical Risk Minimization, ERM）：即在训练数据上最小化损失函数。而最小二乘问题其实也是经验风险最小化问题的一种，甚至机器学习的某些任务（比如回归）本身就是最小二乘问题。经验风险最小化问题是一种通用的函数拟合框架，不过损失函数有所不同，通常使用梯度下降法来进行求解。
那么为什么机器学习中使用梯度下降法来求解，而计算机视觉（SLAM、SfM、相机标定、BA）中使用Gauss-Newton/Levenberg-Marquardt来进行求解呢？这是因为机器学习的问题可以只用关心局部的“好解”，而不用像计算机视觉问题那样需要求解全局的“精确最小值”；另外，机器学习问题规模巨大、结构复杂，使用梯度下降法要简单、健壮、高效的多。
3. 求解
接下来就来介绍一下使用梯度下降法求解非线性最小二乘问题。还是先看非线性最小二乘问题的定义：
其中：
(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)的梯度。根据求导的链式法则：
其中：
(J(theta))：雅可比矩阵（Jacobian），大小为 (m times n)
即目标函数的梯度是：(nabla S(theta) = 2 J(theta)^T mathbf{r}(theta))。
另一方面，在每次梯度下降之后，需要更新参数向量：
其中：
(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})：
#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 到 6 double 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.0 cout << "初始猜测: 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残差平方和tt参数 att参数 bn"; cout << "----t----------tt------tt------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 << "tt" << theta(0) << "tt" << 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 求解的关键还是在于计算雅可比矩阵，对于问题模型函数(f(x; boldsymbol{theta}) = a e ^{bx})来说，雅可比矩阵应该是：
对比代码中的实现：
// 计算 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; } 另外，除了迭代过程中的初始条件和迭代停止条件，控制步长的学习率也需要注意。设置的学习率过小，迭代次数就会很长导致收敛很慢；而设置的学习率过大，就容易略过最优解导致结果不问题。因此，在实际的工程应用中，通常不会使用原始的梯度下降法，而是根据需求使用不同优化版本的梯度下降法。关于这一点，有机会的话会在后续的文章中进一步论述。
