1. 高斯牛顿法发展历程

从上到下为高斯牛顿法的前世今生以及未来的演化：

• 最速下降法（一阶梯度法）

• 牛顿法（二阶梯度法）

• 高斯牛顿法

• 列文伯格法

• 马夸尔特法

2、问题的引出

先来考虑如下优化目标函数：

我们的目标：

对变量x 进行优化，即寻找一组合适的x使得优化目标函数F(x) 最小。

最直观的方法：

然后基于这些这些组x求出F ( x ) 的全部极值，然后通过比较极值的大小就可以找到F ( x ) 的最小值和其所对应的x 取值。

直观方法的问题：

根据上述方法，我们可以找到F ( x ) 全局最小值，但这里面有一个严重的问题。当F ( x ) 为复杂的非线性函数时，其关于x 的一阶导数是十分复杂的，导致 dF/dx = 0 难以求解，或者说求得其全部解是十分困难的，即我们难以获得F ( x ) 全局性质，自然后续的求极值比大小就无法进行，这个方法就走进了死胡同。

解决直观方法问题的思路：

对于不方便直接求解的最小二乘问题，我们可以用迭代的方式，从一个初始值触发，不断地更新当前的优化变量，使目标函数的值下降。

也就是说，虽然我们难以求解 dF/dx = 0，但我们给定一组x，求 dF/dx 的值是很容易的，就是把给定的x 代入 dF/dx 算，批量处理这一过程是计算机最擅长的事情，即我们可以很容易获得F ( x ) 局部性质。那么我们为了解决这样的一个优化问题，我们可以给定一个初值，然后利用F ( x ) 在初值附近的局部性质（即在初值附近x 如何变化可以使得F ( x ) 更小）得到初值的变化增量 ，通过迭代，当Δ x (DeltaX) 足够小时，我们就找到了F(x) 比较小的一个值和其对应的x 取值。与x相对应，这里Δ x( DeltaX)也是 n 维列向量。对应算法流程为：

这让求解导函数为零的问题，变成了一个不断寻找下降增量 的问题。我们将看到，由于可以对 f 进行线性化，增量的计算将简单的多。当函数下降直到增量非常小的时候，就认为算法收敛，目标函数达到了一个极小值。在这个过程中，问题在于如何找到每次迭代点的增量，而这是一个局部的问题，我们只需要关心 f 在迭代值处的局部性质而非全局性质。这类方法在最优化、机器学习等领域广泛应用。

接下来，我们考察如何寻找这个增量 。 这部分知识实际属于数值优化的领域。

这里有几点值得注意：

• 基于局部性质找到的F ( x ) 较小值并不一定是F ( x ) 的最小值，因为我们只考察了F ( x )
  在初值附近的局部情况，没有对F ( x ) 全局性质进行考察。

• 当 Δx (Deltax) 足够小时，算法结束。这个可以这样理解：此时x 只有在变化一点点的情况下才会使F(x) 变小，这说明F(x) 的函数值已经到了一个谷底附近，x 稍微变大一点都会使得F(x)从谷底向上走从而变大。之所以我们不把算法结束条件设置成Δ x = 0 （Delta x = 0） 这是因为我们恰好让F ( x ) 函数值在谷底点几乎不可能，会导致算法一直在谷底点附近小范围内震荡。

• 局部方法的核心在于如何获得增量Δ x （Delta x ）

3、高斯牛顿法的前世

由于原理的相似性，最速下降法和牛顿法可以统称为一阶，二阶梯度法。这两种方法的原理分别是在局部用一次函数，二次函数对非线性函数F ( x ) 进行近似，然后利用近似函数的最小值推测非线性函数F ( x )的极小值。

3.1、一阶，二阶梯度法共有原理

3.2、最速下降法（一阶梯度法）

1）最速下降法原理

最速下降法（一阶梯度法）就是保留泰勒展开的一阶项用来近似非线性函数F ( x )，即：

2）最速下降法缺点

该方法过于贪心，容易走出锯齿线，反而增加迭代次数。

3.3、牛顿法（二阶梯度法）

1）牛顿法的原理

牛顿法（二阶梯度法）就是保留泰勒展开的一阶项和二阶项用来近似非线性函数F(x)，即：

解得：

2）牛顿法的缺点
海塞矩阵 计算量太大了。

3）阻尼牛顿法
阻尼牛顿法就是在使用牛顿法获得增量方向后，进一步对最优步长进行搜索：

4、高斯牛顿法

4.1 高斯牛顿法的思想

高斯牛顿法针对最小二乘问题，采用一定的方法对牛顿法中的海塞矩阵 进行近似，从而简化了计算量。

注意：只有最小二乘问题才能使用高斯牛顿法

4.2 最小二乘问题

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。最小二乘问题可以表述为：

4.3 高斯牛顿法原理

4.4 高斯牛顿法算法流程

4.5 高斯牛顿法缺点

• 在求解增量 时候，我们需要计算 ，但是 只有半正定性质，在实际计算中可能会遇见 为奇异矩阵或者病态矩阵，从而导致增量不稳定，算法不收敛。出现这种现象的深层原因是 F(x) 在 处的近似不像一个二次函数。这也应该是牛顿法的缺点？

  • 奇异矩阵：不满秩的方阵

  • 病态矩阵：求解方程组时如果对数据进行较小的扰动，则得出的结果具有很大波动，这样的矩阵称为病态矩阵。

• 当我们求得的增量 过大时，我们使用泰勒展开进行局部近似就会不准确。从而导致算法不收敛。

5、曲线拟合示例

5.1 C++代码

#include <iostream>
#include <chrono>
#include <opencv2/opencv.hpp>
#include <Eigen/Core>
#include <Eigen/Dense>



using namespace std;
using namespace Eigen;

int Test001LeastSquares();
int LeastSquare002();

int main(int argc, char **argv)
{  
  //在实际的项目中，不会有预设的真值，只有直接获取到的观测值，目的是为了依据观测值估算出真值。对于最小二乘法，是通过计算误差平方和最小来估算真值
  double a_real = 1.0, b_real = 2.0, c_real = 3.0;   // 真实参数值: 这里所谓真值，是为了模拟结合随机数用于构造实际观测值。
  double a_esti = 2.0, b_esti = -1.0, c_esti = 5.0;   // 估计参数值

  std::string str_init_value = std::to_string(a_esti) + ", " + std::to_string(b_esti) + ", " + std::to_string(c_esti) + "";

  int point_num = 100;                                     // 数据点
  double w_sigma = 1.0;                                    // 噪声Sigma(西格玛)值
  double inv_sigma = 1.0 / w_sigma;                        // 噪声Sigma(西格玛)的倒数???
  cv::RNG rng;                                             // OpenCV随机数产生器

  //创建模拟实际观测点
  vector<double> x_data, y_data;
  {
    for (int index_point = 0; index_point < point_num; index_point++)
    {
      double rng_value = rng.gaussian(w_sigma * w_sigma);

      //模拟实际观测点
      double x_temp = (double)index_point / (double)point_num;   //进行 point_num(100) 等分
      double y_temp = exp(a_real * x_temp * x_temp + b_real * x_temp + c_real) + rng_value;

      x_data.push_back(x_temp);
      y_data.push_back(y_temp);
    }
  }


  // 开始Gauss-Newton迭代
  int iterations = 100;           // 迭代次数
  double cost = 0, lastCost = 0;  // 本次迭代的cost和上一次迭代的cost

  chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
  for (int index_iter = 0; index_iter < iterations; index_iter++)
  {
    Matrix3d hessian = Matrix3d::Zero();  //海森矩阵：计算所有点的累加和 Hessian = J^T W^{-1} J in Gauss-Newton
    Vector3d bias = Vector3d::Zero();     //bias: 计算所有点的累加和 bias 偏导?  斜的,偏动的
    cost = 0;                             //误差平方和: 计算所有点的累加和
    {
      //针对海森矩阵(hessian)、bias、误差平方和(cost) ,计算所有点的累加和
      for (int index_points = 0; index_points < point_num; index_points++)
      {
        double x_temp = x_data[index_points]; // 第i个数据点(实际观测点)
        double y_temp = y_data[index_points]; // 第i个数据点(实际观测点)

        Vector3d jacobian; //雅可比矩阵: 针对每个点，计算对应的雅可比矩阵
        {
          jacobian[0] = -x_temp * x_temp * exp(a_esti * x_temp * x_temp + b_esti * x_temp + c_esti);  //求解雅可比矩阵: de/da 对系数a求导, 参考: P133, (6.40)
          jacobian[1] = -x_temp * exp(a_esti * x_temp * x_temp + b_esti * x_temp + c_esti);           //求解雅可比矩阵: de/db 对系数b求导, 参考: P133, (6.40)
          jacobian[2] = -exp(a_esti * x_temp * x_temp + b_esti * x_temp + c_esti);                    //求解雅可比矩阵: de/dc 对系数c求导, 参考: P133, (6.40)
        }

        //海森矩阵：计算所有点的累加和
        hessian += inv_sigma * inv_sigma * jacobian * jacobian.transpose();                           //计算累加和: 高斯牛顿法的增量方程? P133, (6.41) 公式的左侧 ?

        //bias: 计算所有点的累加和
        double error = y_temp - exp(a_esti * x_temp * x_temp + b_esti * x_temp + c_esti);          //误差值
        bias += -inv_sigma * inv_sigma * error * jacobian;                                            //计算累加和: 高斯牛顿法的增量方程? P133, (6.41) 公式的右侧 ?

        //误差平方和: 计算所有点的累加和
        cost += error * error;  //求解最小二乘问题(误差平方和)，估计曲线参数 ? P133, (6.38)
     
      } //end for : for (int index_points = 0; index_points < point_num; index_points++)
      
    }

    if (index_iter > 0 && cost >= lastCost)
    {
      cout << "cost: " << cost << " > last cost: " << lastCost << ", break11." << endl;
      break;
    }


    // 接下来，便是求解 (6.41) 公式 中的 X (delta_X_k) :  求解线性方程 Hx=b
    // 求解线性方程 Hx=b  P129, (公式 6.33)
    Vector3d deltaX = hessian.ldlt().solve(bias);
    if(std::isnan(deltaX[0]))
    {
      cout << "result is nan!" << endl;
      break;
    }

    a_esti += deltaX[0];
    b_esti += deltaX[1];
    c_esti += deltaX[2];

    lastCost = cost;

    cout << "iterations[" << index_iter << "/" << iterations << "] init = " << str_init_value << ", estimated: " << a_esti << ", " << b_esti << ", " << c_esti
         << "\t deltaX(has update estimate): " << deltaX.transpose() << ", "
         << "\t total cost: " << cost << endl;
         
  }//end for : for (int index_iter = 0; index_iter < iterations; index_iter++)

  chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
  chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
  cout << "solve time cost = " << time_used.count() << " seconds. " << endl;

  cout << "estimated abc = " << a_esti << ", " << b_esti << ", " << c_esti << endl;
  
  return 0;

}

5.2 运行结果

5.3 图形化显示

光滑曲线为拟合后的曲线。锯齿形的拟合前观测点。