现在有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$ 。不断迭代直到满足结束条件。

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

首先将误差函数进行一阶泰勒展开，然后取平方（因为线性是没最小值的）。 $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)$ 。
对一阶泰勒的平方以da为变量求导。当导数为零时候da是最优解。求导并令导数等于0后的结果为： $J(a)^Tg(a)+J(a)^TJ(a)da=0$ .于是我们得到了一个线性方程： $J(a)^TJ(a)da=-J(a)^Tg(a)$ 。其中da是变量其他都是常量。所以我们可以解出最优的da。
更新参数a=a+da。
直到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/

实例通俗解释高斯牛顿法求解曲线拟合的最小二乘法问题与c++编程实践教程 要解决的问题是： 怎么解决这个问题 更详细的讲讲怎么用高斯牛顿法最小化误差函数 c++编程实践

怎么解决这个问题

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

c++编程实践

相关推荐

实例通俗解释高斯牛顿法求解曲线拟合的最小二乘法问题与c++编程实践教程要解决的问题是：怎么解决这个问题更详细的讲讲怎么用高斯牛顿法最小化误差函数 c++编程实践