算法#04--详解高斯牛顿迭代法原理和代码

最小二乘法–高斯牛顿迭代法

上一篇文章详解了最小二乘法的线性拟合。本文将详解最小二乘法的非线性拟合，高斯牛顿迭代法。

1.原理

高斯—牛顿迭代法的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。

①已知m个点：

②函数原型：

其中：（m>=n）

③目的是找到最优解β，使得残差平方和最小：

残差：

④要求最小值，即S的对β偏导数等于0：

⑤在非线性系统中，是变量和参数的函数，没有close解。因此，我们给定一个初始值，用迭代法逼近解：

其中k是迭代次数，是迭代矢量。

⑥而每次迭代函数是线性的，在处用泰勒级数展开：

其中：J是已知的矩阵，为了方便迭代，令。

⑦此时残差表示为：

⑧带入公式④有：

化解得：

⑨写成矩阵形式：

⑩所以最终迭代公式为：

其中，Jf是函数f=（x，β）对β的雅可比矩阵。

2.代码

用java代码实现，解维基百科的例子：

https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm

①已知数据：

②函数模型：

③残差公式：

④对β求偏导数：

⑤代码如下：

public class GaussNewton {
    double[] xData = new double[]{
  0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740};
    double[] yData = new double[]{
  0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317};

    double[][] bMatrix = new double[2][1];//系数β矩阵
    int m = xData.length;
    int n = bMatrix.length;
    int iterations = 7;//迭代次数

    //迭代公式求解，即1中公式⑩
    private void magic(){
        //β1,β2迭代初值
        bMatrix[0][0] = 0.9;
        bMatrix[1][0] = 0.2;

        //求J矩阵
        for(int k = 0; k < iterations; k++) { 
            double[][] J = new double[m][n];
            double[][] JT = new double[n][m];
            for(int i = 0; i < m; i++){
                for(int j = 0; j < n; j++) {
                    J[i][j] = secondDerivative(xData[i], bMatrix[0][0], bMatrix[1][0], j);
                }
            }

            JT = MatrixMath.invert(J);//求转置矩阵JT
            double[][] invertedPart = MatrixMath.mult(JT, J);//矩阵JT与J相乘

            //矩阵invertedPart行列式的值：|JT*J|
            double result = MatrixMath.mathDeterminantCalculation(invertedPart);

            //求矩阵invertedPart的逆矩阵:(JT*J)^-1
            double[][] reversedPart = MatrixMath.getIn