本文將詳解最小二乘法的非線性擬合，高斯牛頓迭代法。

1.原理

高斯—牛頓迭代法的基本思想是使用泰勒級(jí)數(shù)展開(kāi)式去近似地代替非線性回歸模型，然后通過(guò)多次迭代，多次修正回歸系數(shù)，使回歸系數(shù)不斷逼近非線性回歸模型的最佳回歸系數(shù)，最后使原模型的殘差平方和達(dá)到最小。

①已知m個(gè)點(diǎn)：

②函數(shù)原型：

其中：（m>=n）

③目的是找到最優(yōu)解β，使得殘差平方和最小：

殘差：

④要求最小值，即S的對(duì)β偏導(dǎo)數(shù)等于0：

⑤在非線性系統(tǒng)中，是變量和參數(shù)的函數(shù)，沒(méi)有close解。因此，我們給定一個(gè)初始值，用迭代法逼近解：

其中k是迭代次數(shù)，是迭代矢量。

⑥而每次迭代函數(shù)是線性的，在處用泰勒級(jí)數(shù)展開(kāi)：

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

⑦此時(shí)殘差表示為：

⑧帶入公式④有：

化解得：

⑨寫(xiě)成矩陣形式：

⑩所以最終迭代公式為：

其中，Jf是函數(shù)f=（x，β）對(duì)β的雅可比矩陣。

2.代碼

用java代碼實(shí)現(xiàn)，解維基百科的例子：

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

①已知數(shù)據(jù)：

②函數(shù)模型：

③殘差公式：

④對(duì)β求偏導(dǎo)數(shù)：

⑤代碼如下：

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];//系數(shù)β矩陣int m = xData.length;int n = bMatrix.length;int iterations = 7;//迭代次數(shù)//迭代公式求解，即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);//求轉(zhuǎn)置矩陣JTdouble[][] invertedPart = MatrixMath.mult(JT, J);//矩陣JT與J相乘//矩陣invertedPart行列式的值：|JT*J|double result = MatrixMath.mathDeterminantCalculation(invertedPart);//求矩陣invertedPart的逆矩陣:(JT*J)^-1double[][] reversedPart = MatrixMath.getInverseMatrix(invertedPart, result);//求方差r(β)矩陣: ri = yi - f(xi, b)double[][] residuals = new double[m][1];for(int i = 0; i < m; i++) {residuals[i][0] = yData[i] - (bMatrix[0][0] * xData[i]) / (bMatrix[1][0] + xData[i]);}//求矩陣積reversedPart*JT*residuals: (JT*J)^-1*JT*rdouble[][] product = MatrixMath.mult(MatrixMath.mult(reversedPart, JT), residuals);//迭代公式, 即公式⑩double[][] newB = MatrixMath.plus(bMatrix, product);bMatrix = newB; } //顯示系數(shù)值System.out.println("b1: " + bMatrix[0][0] + "\nb2: " + bMatrix[1][0]); }//2中公式④private static double secondDerivative(double x, double b1, double b2, int index){switch(index) {case 0: return x / (b2 + x);//對(duì)系數(shù)bi求導(dǎo)case 1: return -1 * (b1 * x) / Math.pow((b2+x), 2);//對(duì)系數(shù)b2求導(dǎo)}return 0;}public static void main(String[] args) {GaussNewton app = new GaussNewton();app.magic(); } }

運(yùn)行，輸出得到：

b1: 0.3618366954234483?
b2: 0.5562654497238557

⑥其中用到的矩陣運(yùn)算代碼如下：

public class MatrixMath {/*** 矩陣基本運(yùn)算：加、減、乘、除、轉(zhuǎn)置*/public final static int OPERATION_ADD = 1; public final static int OPERATION_SUB = 2; public final static int OPERATION_MUL = 4; /*** 矩陣加法* @param a 加數(shù)* @param b 被加數(shù)* @return 和*/public static double[][] plus(double[][] a, double[][] b){if(legalOperation(a, b, OPERATION_ADD)) { for(int i=0; i<a.length; i++) { for(int j=0; j<a[0].length; j++) { a[i][j] = a[i][j] + b[i][j]; }} }return a;}/*** 矩陣減法* @param a 減數(shù)* @param b 被減數(shù)* @return 差*/public static double[][] substract(double[][] a, double[][] b){if(legalOperation(a, b, OPERATION_SUB)) { for(int i=0; i<a.length; i++) { for(int j=0; j<a[0].length; j++) { a[i][j] = a[i][j] - b[i][j]; }} }return a;} /*** 判斷矩陣行列是否符合運(yùn)算* @param a 進(jìn)行運(yùn)算的矩陣* @param b 進(jìn)行運(yùn)算的矩陣* @param type 運(yùn)算類(lèi)型* @return 符合/不符合運(yùn)算*/private static boolean legalOperation(double[][] a, double[][] b, int type) { boolean legal = true; if(type == OPERATION_ADD || type == OPERATION_SUB) { if(a.length != b.length || a[0].length != b[0].length) { legal = false; } } else if(type == OPERATION_MUL) { if(a[0].length != b.length) { legal = false; } } return legal; } /*** 矩陣乘法* @param a 乘數(shù)* @param b 被乘數(shù)* @return 積*/public static double[][] mult(double[][] a, double[][] b){if(legalOperation(a, b, OPERATION_MUL)) {double[][] result = new double[a.length][b[0].length];for(int i=0; i< a.length; i++) { for(int j=0; j< b[0].length; j++) { result[i][j] = calculateSingleResult(a, b, i, j); }}return result;}else{return null;} }/*** 矩陣乘法* @param a 乘數(shù)* @param b 被乘數(shù)* @return 積*/public static double[][] mult(double[][] a, int b) { for(int i=0; i<a.length; i++) { for(int j=0; j<a[0].length; j++) { a[i][j] = a[i][j] * b; } } return a; }/*** 乘數(shù)矩陣的行元素與被乘數(shù)矩陣列元素積的和* @param a 乘數(shù)矩陣* @param b 被乘數(shù)矩陣* @param row 行* @param column 列* @return 值*/private static double calculateSingleResult(double[][] a, double[][] b, int row, int column) { double result = 0.0; for(int k = 0; k< a[0].length; k++) { result += a[row][k] * b[k][column]; } return result; } /*** 矩陣的轉(zhuǎn)置* @param a 要轉(zhuǎn)置的矩陣* @return 轉(zhuǎn)置后的矩陣*/public static double[][] invert(double[][] a){double[][] result = new double[a[0].length][a.length];for(int i=0;i<a.length;i++){for(int j=0;j<a[0].length;j++){ result[j][i]=a[i][j]; } } return result;} /** * 求可逆矩陣（使用代數(shù)余子式的形式） */ /** * 求傳入的矩陣的逆矩陣 * @param value 需要轉(zhuǎn)換的矩陣 * @return 逆矩陣 */ public static double[][] getInverseMatrix(double[][] value,double result){ double[][] transferMatrix = new double[value.length][value[0].length]; //計(jì)算代數(shù)余子式，并賦值給|A| for (int i = 0; i < value.length; i++) { for (int j = 0; j < value[i].length; j++) { transferMatrix[j][i] = mathDeterminantCalculation(getNewMatrix(i, j, value)); if ((i+j)%2!=0) { transferMatrix[j][i] = -transferMatrix[j][i]; } transferMatrix[j][i] /= result; } } return transferMatrix; } /*** * 求行列式的值* @param value 需要算的行列式 * @return 計(jì)算的結(jié)果 */ public static double mathDeterminantCalculation(double[][] value){ if (value.length == 1) { //當(dāng)行列式為1階的時(shí)候就直接返回本身 return value[0][0]; }if (value.length == 2) { //如果行列式為二階的時(shí)候直接進(jìn)行計(jì)算 return value[0][0]*value[1][1]-value[0][1]*value[1][0]; } //當(dāng)行列式的階數(shù)大于2時(shí) double result = 1; for (int i = 0; i < value.length; i++) { //檢查數(shù)組對(duì)角線位置的數(shù)值是否是0，如果是零則對(duì)該數(shù)組進(jìn)行調(diào)換，查找到一行不為0的進(jìn)行調(diào)換 if (value[i][i] == 0) { value = changeDeterminantNoZero(value,i,i); result*=-1; } for (int j = 0; j <i; j++) { //讓開(kāi)始處理的行的首位為0處理為三角形式 //如果要處理的列為0則和自己調(diào)換一下位置，這樣就省去了計(jì)算 if (value[i][j] == 0) { continue; } //如果要是要處理的行是0則和上面的一行進(jìn)行調(diào)換 if (value[j][j]==0) { double[] temp = value[i]; value[i] = value[i-1]; value[i-1] = temp; result*=-1; continue; } double ratio = -(value[i][j]/value[j][j]); value[i] = addValue(value[i],value[j],ratio); } } return mathValue(value,result);} /** * 計(jì)算行列式的結(jié)果 * @param value * @return */ private static double mathValue(double[][] value,double result){ for (int i = 0; i < value.length; i++) { //如果對(duì)角線上有一個(gè)值為0則全部為0，直接返回結(jié)果 if (value[i][i]==0) { return 0; } result *= value[i][i]; } return result; } /*** * 將i行之前的每一行乘以一個(gè)系數(shù)，使得從i行的第i列之前的數(shù)字置換為0 * @param currentRow 當(dāng)前要處理的行 * @param frontRow i行之前的遍歷的行 * @param ratio 要乘以的系數(shù) * @return 將i行i列之前數(shù)字置換為0后的新的行 */ private static double[] addValue(double[] currentRow,double[] frontRow, double ratio){ for (int i = 0; i < currentRow.length; i++) { currentRow[i] += frontRow[i]*ratio; } return currentRow; } /** * 指定列的位置是否為0，查找第一個(gè)不為0的位置的行進(jìn)行位置調(diào)換，如果沒(méi)有則返回原來(lái)的值 * @param determinant 需要處理的行列式 * @param line 要調(diào)換的行 * @param row 要判斷的列 */ private static double[][] changeDeterminantNoZero(double[][] determinant,int line,int column){ for (int i = line; i < determinant.length; i++) { //進(jìn)行行調(diào)換 if (determinant[i][column] != 0) { double[] temp = determinant[line]; determinant[line] = determinant[i]; determinant[i] = temp; return determinant; } } return determinant; } /** * 轉(zhuǎn)換為代數(shù)余子式 * @param row 行 * @param line 列 * @param matrix 要轉(zhuǎn)換的矩陣 * @return 轉(zhuǎn)換的代數(shù)余子式 */ private static double[][] getNewMatrix(int row,int line,double[][] matrix){ double[][] newMatrix = new double[matrix.length-1][matrix[0].length-1]; int rowNum = 0,lineNum = 0; for (int i = 0; i < matrix.length; i++) { if (i == row){ continue; } for (int j = 0; j < matrix[i].length; j++) { if (j == line) { continue; } newMatrix[rowNum][lineNum++%(matrix[0].length-1)] = matrix[i][j]; } rowNum++; } return newMatrix; } public static void main(String[] args) { //double[][] test = {{0,0,0,1,2},{0,0,0,2,3},{1,1,0,0,0},{0,1,1,0,0},{0,0,1,0,0}}; double[][] test = {{3.8067488033632655, -2.894113667134647}, {-2.894113667134647, 3.6978894069779504}};double result; try { double[][] temp = new double[test.length][test[0].length]; for (int i = 0; i < test.length; i++) { for (int j = 0; j < test[i].length; j++) { temp[i][j] = test[i][j]; } } //先計(jì)算矩陣的行列式的值是否等于0，如果不等于0則該矩陣是可逆的 result = mathDeterminantCalculation(temp); if (result == 0) { System.out.println("矩陣不可逆"); }else { System.out.println("矩陣可逆"); //求出逆矩陣 double[][] result11 = getInverseMatrix(test,result); //打印逆矩陣 for (int i = 0; i < result11.length; i++) { for (int j = 0; j < result11[i].length; j++) { System.out.print(result11[i][j]+" "); } System.out.println(); } } } catch (Exception e) { e.printStackTrace(); System.out.println("不是正確的行列式！！"); } } }

?

總結(jié)

以上是生活随笔為你收集整理的高斯牛顿迭代法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。