最小二乘問(wen)題(ti)詳(xiang)解6：梯度下降法
1. 引言
在之(zhi)前的兩篇(pian)文章、中，筆(bi)者介紹(shao)了非(fei)線(xian)(xian)性最小(xiao)二乘問題(ti)，并使用Gauss-Newton方法(fa)來(lai)進行(xing)求解。不過，求解非(fei)線(xian)(xian)性最小(xiao)二乘問題(ti)還有(you)另外一種(zhong)方法(fa)——梯度下降法(fa)。
2. 背景
梯度下降法在人工智能的機器學習中使用的非常多，因為機器學習的訓練過程通常被形式化為經驗風險最小化問題（Empirical Risk Minimization, ERM）：即(ji)在(zai)訓練數據上(shang)最(zui)(zui)(zui)小(xiao)化損失函數。而(er)最(zui)(zui)(zui)小(xiao)二乘(cheng)問(wen)(wen)題其實也(ye)是經(jing)(jing)驗風(feng)險最(zui)(zui)(zui)小(xiao)化問(wen)(wen)題的(de)一(yi)種，甚至機器學習的(de)某些(xie)任務（比如回歸）本身就是最(zui)(zui)(zui)小(xiao)二乘(cheng)問(wen)(wen)題。經(jing)(jing)驗風(feng)險最(zui)(zui)(zui)小(xiao)化問(wen)(wen)題是一(yi)種通用(yong)的(de)函數擬(ni)合(he)框(kuang)架，不過(guo)損失函數有(you)所(suo)不同，通常使用(yong)梯度下降(jiang)法(fa)來進行求解。
那(nei)(nei)么(me)(me)為什么(me)(me)機器學(xue)習中(zhong)使用(yong)梯度下(xia)降法(fa)來求(qiu)解(jie)，而(er)計算機視(shi)覺(jue)（SLAM、SfM、相機標定、BA）中(zhong)使用(yong)Gauss-Newton/Levenberg-Marquardt來進行求(qiu)解(jie)呢？這是(shi)因為機器學(xue)習的問題(ti)可(ke)以(yi)只(zhi)用(yong)關心局部的“好(hao)解(jie)”，而(er)不用(yong)像計算機視(shi)覺(jue)問題(ti)那(nei)(nei)樣需要求(qiu)解(jie)全(quan)局的“精確(que)最小值(zhi)”；另外，機器學(xue)習問題(ti)規模巨(ju)大、結構復雜，使用(yong)梯度下(xia)降法(fa)要簡單、健壯、高效的多。
3. 求解
接下(xia)來(lai)就來(lai)介(jie)紹一下(xia)使(shi)用梯度下(xia)降法求解非(fei)線性最(zui)(zui)小二乘問題(ti)。還是先看非(fei)線性最(zui)(zui)小二乘問題(ti)的(de)定(ding)義：
其中：
\(\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)\)：目標函數（損(sun)失函數），是我們要最小(xiao)化的殘差平方和
梯度下降法的核心思想是：在當前點，沿著目標函數下降最快的方向走一步，然后重復。而這個“最快下降方向”就是負梯度方向\(-\nabla S(\theta)\)。因此問題的關鍵在于計算目標函數\(S(\theta) = \ \mathbf{r}(\theta)\ ^2\)的(de)梯度。根據求導的(de)鏈(lian)式法(fa)則：
其中：
\(J(\theta)\)：雅可比矩陣（Jacobian），大小為 \(m \times n\)
即目標函數的梯度是：\(\nabla S(\theta) = 2 J(\theta)^T \mathbf{r}(\theta)\)。
另一方面，在每次梯度下降之后，需要更(geng)新參(can)數向量：
其中：
\(\theta_k\)：第 \(k\) 次迭代的參數
\(\alpha > 0\)：學習率（step size），控制步長
\(J_k = J(\theta_k)\)，\(\mathbf{r}_k = \mathbf{r}(\theta_k)\)
因此(ci)，將梯度下降方法完整的流程總結(jie)如下：
- 初始化：選一個初始猜測 θ?
- 設置學習率 α（例如 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殘差平方和\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
求解的關鍵還是在于計算雅可比矩陣，對于問題模型函數\(f(x; \boldsymbol{\theta}) = a e ^{bx}\)來說，雅可比矩(ju)陣應該(gai)是：
對比代碼中的實現：
// 計算 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;
}
另外(wai)，除了迭代過程中(zhong)的(de)(de)(de)初(chu)始條件和迭代停止條件，控制步長(chang)的(de)(de)(de)學習率也需要注意。設(she)置的(de)(de)(de)學習率過小，迭代次(ci)數就(jiu)會很長(chang)導(dao)致收(shou)斂很慢；而設(she)置的(de)(de)(de)學習率過大，就(jiu)容(rong)易略過最優解導(dao)致結(jie)果不問題。因此，在實際的(de)(de)(de)工(gong)程應用中(zhong)，通常不會使(shi)用原始的(de)(de)(de)梯度(du)下(xia)降法(fa)，而是(shi)根據需求使(shi)用不同優化版(ban)本的(de)(de)(de)梯度(du)下(xia)降法(fa)。關于(yu)這一(yi)點，有(you)機會的(de)(de)(de)話會在后續(xu)的(de)(de)(de)文章(zhang)中(zhong)進一(yi)步論述。
