【数学与算法】非线性最小二乘法的解法【最速梯度下降法】 【牛顿法】 【高斯牛顿法

关于非线性优化的问题，可以推荐观看视觉SLAM十四讲视频的第六讲 非线性优化。

如果不明白线性和非线性，可参考这篇博客：线性最小二乘和非线性最小二乘

这篇博客的后面有讲到几种优化方法(最速梯度下降法、牛顿法、高斯牛顿法，LM算法)，很容易记住，不像其他的公式推导那么生硬：Bundle Adjustment—即最小化重投影误差（高翔slam—第七讲）

这篇博客也很棒非线性优化（高翔slam—第六讲 ）

下图摘自非线性优化

非线性最小二乘:

所谓非线性，就是f(x)无法表示为的线性关系，而是某种非线性关系。

这让求解导函数为零的问题变成了一个不断寻找下降增量 ∆xk 的问题。以下就是寻找增量的方法。

下面是对平方展开的：

一阶梯度法(梯度下降法):

不需要求步长(设定固定步长)叫梯度下降法，

每一步迭代都需要求最优步长的叫做最速下降法。

二阶梯度法(又叫牛顿法):

牛顿法与高斯牛顿法的不同之处：

牛顿法是直接对误差的平方∣∣f(x+Δx)∣∣22\color{blue}||f(x+\Delta{x})||^2_2∣∣f(x+Δx)∣∣22​ 在x处进行泰勒展开。

高斯牛顿法是对误差函数f(x+Δx)f(x+\Delta{x})f(x+Δx)在x处进行泰勒展开后得到f(x+∆x)=f(x)+J(x)∆x\color{blue}f(x+∆x) = f(x)+ J(x)∆xf(x+∆x)=f(x)+J(x)∆x，再对f(x+∆x)=f(x)+J(x)∆x\color{blue}f(x+∆x) = f(x)+ J(x)∆xf(x+∆x)=f(x)+J(x)∆x进行平方操作。

评价一阶梯度与二阶梯度：

一阶梯度法(梯度下降法)，太慢了，因为它每次都会找到最陡的方向进行下降。如果是最速下降法(比梯度下降法每一步迭代多了求最优步长)，前后两次迭代的梯度向量方向是正交的，也就是说一直走的是直角，故就算很简单的也会走很多步，即过于贪婪(zigzag 问题)，过于贪心，容易走出锯齿路线，反而增加了迭代次数；

而二阶梯度法(牛顿法)，迭代次数少，对高阶表现良好，但不可避免的我们需要计算目标函数的 H 矩阵，这在问题规模较大时非常困难，我们通常倾向于避免 H 的计算。

高斯牛顿法:

令H=J(x)J(x)T，g=−J(x)f(x)\color{blue}H=J(x)J(x)^T，g=-J(x)f(x)H=J(x)J(x)T，g=−J(x)f(x)，则高斯牛顿方程变为：

H∆x=gH∆x=gH∆x=g

重中之重，牛顿法和高斯牛顿法的对比：对比牛顿法中 H∗∆x=−J\color{blue}H*∆x=-JH∗∆x=−J，高斯牛顿法用J(x)J(x)T\color{blue}J(x)J(x)^TJ(x)J(x)T作为牛顿法中Hessian矩阵的近似，从而省略了H\color{blue}HH的计算。

求解高斯牛顿方程是整个优化问题的核心所在，如果能解出该方程，则高斯牛顿法的步骤如下：

为了求解增量方程，需要解出H−1\color{blue}H^{-1}H−1，这需要H\color{blue}HH矩阵可逆。但是实际上H\color{blue}HH只有半正定，也就是H\color{blue}HH可能会是奇异矩阵或ill-condition的情况，此时增量的稳定性较差，导致算法不收敛。就算H非奇异也非病态，但是如果求出来的步长 ∆x\color{blue}∆x∆x 太大，也无法保证能迭代收敛。

牛顿法与高斯牛顿法的不同之处：

牛顿法是直接对误差的平方∣∣f(x+Δx)∣∣22\color{blue}||f(x+\Delta{x})||^2_2∣∣f(x+Δx)∣∣22​ 在x处进行泰勒展开。

高斯牛顿法是对误差函数f(x+Δx)f(x+\Delta{x})f(x+Δx)在x处进行泰勒展开后得到f(x+∆x)=f(x)+J(x)∆x\color{blue}f(x+∆x) = f(x)+ J(x)∆xf(x+∆x)=f(x)+J(x)∆x，再对f(x+∆x)=f(x)+J(x)∆x\color{blue}f(x+∆x) = f(x)+ J(x)∆xf(x+∆x)=f(x)+J(x)∆x进行平方操作。

列文伯格-马夸尔特(Levenberg-Marquadt):

例题：使用高斯牛顿法进行曲线拟合

#include <iostream>#include <opencv2/opencv.hpp>#include <Eigen/Core>#include <Eigen/Dense>using namespace std;using namespace Eigen;//g++ test.cpp `pkg-config opencv --libs --cflags` -std=c++11 -o testint main(int argc, char **argv) {double ar = 1.0, br = 2.0, cr = 1.0; // 真实参数值double ae = 2.0, be = -1.0, ce = 5.0; // 估计参数值int N = 100; // 数据点double w_sigma = 1.0; // 噪声Sigma值double inv_sigma = 1.0 / w_sigma;cv::RNG rng; // OpenCV随机数产生器vector<double> x_data, y_data;// 数据for (int i = 0; i < N; i++) {double x = i / 100.0;x_data.push_back(x);y_data.push_back(exp(ar * x * x + br * x + cr) + rng.gaussian(w_sigma * w_sigma));}// 开始Gauss-Newton迭代int iterations = 100; // 迭代次数double cost = 0, lastCost = 0; // 本次迭代的cost和上一次迭代的costfor (int iter = 0; iter < iterations; iter++) {Matrix3d H = Matrix3d::Zero(); // Hessian = J^T W^{-1} J in Gauss-NewtonVector3d b = Vector3d::Zero(); // biascost = 0;for (int i = 0; i < N; i++) {double xi = x_data[i], yi = y_data[i]; // 第i个数据点double error = yi - exp(ae * xi * xi + be * xi + ce);Vector3d J; // 雅可比矩阵J[0] = -xi * xi * exp(ae * xi * xi + be * xi + ce); // de/daJ[1] = -xi * exp(ae * xi * xi + be * xi + ce); // de/dbJ[2] = -exp(ae * xi * xi + be * xi + ce); // de/dcH += inv_sigma * inv_sigma * J * J.transpose();b += -inv_sigma * inv_sigma * error * J;cost += error * error;}// 求解线性方程 Hx=bVector3d dx = H.ldlt().solve(b);if (isnan(dx[0])) {cout << "result is nan!" << endl;break;}if (iter > 0 && cost >= lastCost) {cout << "cost: " << cost << ">= last cost: " << lastCost << ", break." << endl;break;}ae += dx[0];be += dx[1];ce += dx[2];lastCost = cost;cout << "total cost: " << cost << ", \t\tupdate: " << dx.transpose() <<"\t\testimated params: " << ae << "," << be << "," << ce << endl;}cout << "estimated abc = " << ae << ", " << be << ", " << ce << endl;return 0;}

【数学与算法】非线性最小二乘法的解法【最速梯度下降法】 【牛顿法】 【高斯牛顿法】 【LM算法】