实例通俗解释高斯牛顿法求解曲线拟合的最小二乘法问题与c++编程实践教程

要解决的问题是：

现在有N个数据点(x,y)。我们假设这个曲线$y=exp(ax^2+bx+c)+noise$可以拟合那堆数据，其中a,b,c是待求解的参数，noise是噪声。我们要根据那堆数据去算出a,b,c的值。用的方法是高斯牛顿法。为啥有个牛顿？因为它和牛顿法一样都是用泰勒展开，只不过高斯牛顿法是一阶泰勒展开。一阶泰勒展开意味着它是线性方程，所以需要用高斯消元法去解方程。因此名字中的高斯就是这么来的。

怎么解决这个问题

现在我们知道了数据的模型，和数据(x,y)。a,b,c是待求解的参数。那么怎么知道a,b,c是设置的是适合这个数据还是不适合呢？答：计算误差不就可以了么。假设第i个样本数据是$(x_i,y_i)$，那么现在我们给定a,b,c值下的模型误差为：$[exp(ax_i^2+bx_i+c)+noise-y_i]^2$。由于二次方求导会前面有个系数2，为了求导方便我们习惯性在误差前面乘个$\frac 1 2$。这就是我们经常看到的$\frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$. 由于不是只有一个样本。我们当然希望整个样本的误差都很小。所以要将所有样本误差累加起来，以衡量现在我们设定的a,b,c参数值是不是不错。于是就得到了最小二乘法的终极版目标函数$\sum_{i=1}^{N} \frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$

我们的目标是：找到参数a,b,c合适的取值来最小化这个误差$\sum_{i=1}^{N} \frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$。最小二乘法中的最小就是这么来的。

我们来根据那个误差函数找最优的a,b,c的方法叫做高斯牛顿法。接下来我们看看高斯牛顿法怎么做的。
注意：我们是根据误差函数$\sum_{i=1}^{N} \frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$找最优的a,b,c。最小化误差函数。

重复注意：我们是根据误差函数$\sum_{i=1}^{N} \frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$找最优的a,b,c。最小化误差函数。
重复注意：我们是根据误差函数$\sum_{i=1}^{N} \frac 1 2[exp(ax_i^2+bx_i+c)+noise-y_i]^2$找最优的a,b,c。最小化误差函数。

高斯牛顿法的思想就是虽然我不知道a,b,c是多少。但是我可以猜啊。怎么猜？首先第一次猜的时候当然要随便给a,b,c赋值。假如第0次猜a的值为$a=0.2333$，那么接下来我往哪猜比较好？数学家连猜都是认真的。假设我猜下次的a的值为$a=0.2333+da$，现在我们要确定的是$da$到底取多少是比较好的，也就是说$da$现在是个不确定量.数学家想到一个办法能不能在$a=0.2333$这地方对原先函数以$da$为变量进行泰勒展 开来代替原函数。我们知道一阶泰勒展开是酱紫$g(a_0+da)=g(a_0)+g'(a_0)da$。我们可以对误差函数进行泰勒展开，用展开后的式子代替原来的误差函数，我们目标反正是求误差函数的最小值。那么我们$da$的最优取值只需要求$g(a_0+da)=g(a_0)+g'(a_0)da$这个式子的最小值时候$da$的取值即可。

我们看看$g(a_0+da)=g(a_0)+g'(a_0)da$这个式子是不是关于$da$的二次函数？那就非常简单了。因为二次函数的最小值是在它导数等于0的那个点取得的。所以我们只需要对误差函数的一阶泰勒展开以变量$da$来求导，然后令导数等于0.求得最优的$da$。然后就找到下一次猜的$a_{下次}=a_{上次}+da$。不断迭代直到满足结束条件。

更详细的讲讲怎么用高斯牛顿法最小化误差函数

1. 首先将误差函数进行一阶泰勒展开，然后取平方（因为线性是没最小值的）。$\frac 1 2 (\sum_{i=1}^{N}[exp(ax_i^2+bx_i+c)+noise-y_i]^2$+$\sum_{i=1}^{N} [exp(ax_i^2+bx_i+c)+noise-y_i](x_i^2exp(ax_i^2+bx_i+c))da)^2$。为了简化这个公式我们还是把它写成矩阵的形式比较好。我们假设误差函数的泰勒展开为$g(a+da)=g(a)+J(a)da$。然后泰勒展开的平方那就对应矩阵的2范数即$\frac 1 2 ||g(a)+J(a)da||_2^2=\frac 1 2 (g(a)+J(a)da)^T(g(a)+J(a)da)=\frac 1 2 (||g(a)||^2+g(a)^TJ(a)da+da^TJ(a)^Tg(a)+da^TJ(a)^TJ(a)da)$。
2. 对一阶泰勒的平方以da为变量求导。当导数为零时候da是最优解。求导并令导数等于0后的结果为：$J(a)^Tg(a)+J(a)^TJ(a)da=0$.于是我们得到了一个线性方程：$J(a)^TJ(a)da=-J(a)^Tg(a)$。其中da是变量其他都是常量。所以我们可以解出最优的da。
3. 更新参数a=a+da。
4. 直到da足够小，才停止迭代。
   注意：高斯牛顿法存在一些情况是解不出来的。比如$$J(a)^TJ(a)$是不可逆矩阵那就很尴尬。

c++编程实践

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

using namespace std;
using namespace Eigen;

int 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值
    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));
    }

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

    for (int iter = 0; iter < iterations; iter++) {

        Matrix3d H = Matrix3d::Zero();             // Hessian = J^T J in Gauss-Newton
        Vector3d b = Vector3d::Zero();             // bias
        cost = 0;

        for (int i = 0; i < N; i++) {
            double xi = x_data[i], yi = y_data[i];  // 第i个数据点
            // start your code here
            double error = 0;   // 第i个数据点的计算误差
            error = 0; // 填写计算error的表达式
            Vector3d J; // 雅可比矩阵
            J[0] = 0;  // de/da
            J[1] = 0;  // de/db
            J[2] = 0;  // de/dc

            H += J * J.transpose(); // GN近似的H
            b += -error * J;
            // end your code here

            cost += error * error;
        }

        // 求解线性方程 Hx=b，建议用ldlt
 	// start your code here
        Vector3d dx;
	// end your code here

        if (isnan(dx[0])) {
            cout << "result is nan!" << endl;
            break;
        }

        if (iter > 0 && cost > lastCost) {
            // 误差增长了，说明近似的不够好
            cout << "cost: " << cost << ", last cost: " << lastCost << endl;
            break;
        }

        // 更新abc估计值
        ae += dx[0];
        be += dx[1];
        ce += dx[2];

        lastCost = cost;

        cout << "total cost: " << cost << endl;
    }

    cout << "estimated abc = " << ae << ", " << be << ", " << ce << endl;
    return 0;
}

参考文献：
[1] http://fourier.eng.hmc.edu/e176/lectures/NM/node36.html
[2] https://web.stanford.edu/~boyd/cvxbook/