算法#03--具體解釋最小二乘法原理和代碼

最小二乘法原理

最小二乘法的目標：求誤差的最小平方和，相應有兩種：線性和非線性。線性最小二乘的解是closed-form（例如以下文），而非線性最小二乘沒有closed-form，通經常使用迭代法求解（如高斯牛頓迭代法，本文不作介紹）。

【首先得到線性方程組】

1.概念

最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的平方和尋找數據的最佳函數匹配。

利用最小二乘法能夠簡便地求得未知的數據，並使得這些求得的數據與實際數據之間誤差的平方和為最小。

最小二乘法還可用於曲線擬合。

2.原理

函數原型：

已知：

（x0，y0）。（x1。y1）…（xi，yi）…（xn，yn）個點，n>=k。

偏差平方和：

偏差平方和最小值能夠通過使偏導數等於零得到：

簡化左邊等式有：

寫成矩陣形式：公式①

將這個范德蒙得矩陣化簡后可得到：公式②

也就是說X*A=Y，那么A = (X’*X)-1*X’*Y，便得到了系數矩陣A，同一時候，我們也就得到了擬合曲線。

高斯消元法

【然后解線性方程組，即公式①】

1.概念

數學上，高斯消元法（或譯：高斯消去法）（英語：Gaussian Elimination），是線性代數中的一個算法，可用來為線性方程組求解，求出矩陣的秩，以及求出可逆方陣的逆矩陣。當用於一個矩陣時。高斯消元法會產生出一個“行梯陣式”。

2.原理

3.偽代碼

這個算法和上面談到的有點不同。它由絕對值最大的部分開始做起。這樣能夠改善算法的穩定性。

本算法由左至右地計算。每作出下面三個步驟，才跳到下一列和下一行：

• 定出i列的絕對值最大的一個非0的數，將第i行的值與該行交換，使得該行擁有該列的最大值；
• 將i列的數字除以該數，使得i列i行的數成為1。
• 第(i+1)行下面（包含第(j+1)行）全部元素都轉化為0。

全部步驟完畢后，這個矩陣會變成一個行梯矩陣，再用代入法就能夠求解該方程組。

 i = 1
 j = 1
 while (i ≤ m and j ≤ n) do
   Find pivot in column j, starting in row i    // 從第i行開始。找出第j列中的最大值（i、j值應保持不變） 
   maxi = i
   for k = i+1 to m do
     if abs(A[k,j]) > abs(A[maxi,j]) then
       maxi = k   // 使用交換法找出最大值（絕對值最大）
     end if
   end for
   if A[maxi,j] ≠ 0 then  // 判定找到的絕對值最大值是否為零：若不為零就進行下面操作；若為零則說明該列第(i+1)行下面（包含第(i+1)行）均為零，不須要再處理，直接跳轉至第(j+1)列第(i+1)行
     swap rows i and maxi, but do not change the value of i   // 將第i行與找到的最大值所在行做交換，保持i值不變（i值記錄了本次操作的起始行）
     Now A[i,j] will contain the old value of A[maxi,j].
     divide each entry in row i by A[i,j]    // 將交換后的第i行歸一化（第i行全部元素分別除以A[i,j]）
     Now A[i,j] will have the value 1.
     for u = i+1 to m do    // 第j列中，第(i+1)行下面（包含第(i+1)行）全部元素都減去A[i,j]，直到第j列的i+1行以後元素均為零
       subtract A[u,j] * row i from row u
       Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
     end for
     i = i + 1   
   end if
   j = j + 1  // 第j列中。第(i+1)行下面（包含第(i+1)行）全部元素均為零。
移至第(j+1)列，從第(i+1)行開始反復上述步驟。
 end while

代碼

public class CurveFitting {
     ///<summary>
    ///最小二乘法擬合二元多次曲線
    ///比如y=ax+b
    ///當中MultiLine將返回a。b兩個參數。
 ///a相應MultiLine[1] ///b相應MultiLine[0] ///</summary> ///<param name="arrX">已知點的x坐標集合</param> ///<param name="arrY">已知點的y坐標集合</param> ///<param name="length">已知點的個數</param> ///<param name="dimension">方程的最高次數</param> public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension) { int n = dimension + 1; //dimension次方程須要求 dimension+1個 系數  double[][] Guass = new double[n][n + 1]; for (int i = 0; i < n; i++){ //求矩陣公式① int j; for (j = 0; j < n; j++){ Guass[i][j] = SumArr(arrX, j + i, length);//公式①等號左邊第一個矩陣，即Ax=b中的A } Guass[i][j] = SumArr(arrX, i, arrY, 1, length);//公式①等號右邊的矩陣，即Ax=b中的b } return ComputGauss(Guass, n);//高斯消元法 } //求數組的元素的n次方的和，即矩陣A中的元素 private static double SumArr(double[] arr, int n, int length) { double s = 0; for (int i = 0; i < length; i++){ if (arr[i] != 0 || n != 0){ s = s + Math.pow(arr[i], n); } else{ s = s + 1; } } return s; } //求數組的元素的n次方的和，即矩陣b中的元素 private static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length) { double s = 0; for (int i = 0; i < length; i++) { if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0)) s = s + Math.pow(arr1[i], n1) * Math.pow(arr2[i], n2); else s = s + 1; } return s; } //高斯消元法解線性方程組 private static double[] ComputGauss(double[][] Guass, int n) { int i, j; int k, m; double temp; double max; double s; double[] x = new double[n]; for (i = 0; i < n; i++) { x[i] = 0.0;//初始化 } for (j = 0; j < n; j++) { max = 0; k = j; // 從第i行開始，找出第j列中的最大值（i、j值應保持不變）  for (i = j; i < n; i++) { if (Math.abs(Guass[i][j]) > max){ max = Guass[i][j];// 使用交換法找出最大值（絕對值最大） k = i; } } if (k != j) { //將第j行與找到的最大值所在行做交換。保持i值不變（j值記錄了本次操作的起始行） for (m = j; m < n + 1; m++) { temp = Guass[j][m]; Guass[j][m] = Guass[k][m]; Guass[k][m] = temp; } } if (max == 0) { // "此線性方程為神秘線性方程"  return x; } // 第m列中，第(j+1)行下面（包含第(j+1)行）全部元素都減去Guass[j][m] * s / (Guass[j][j]) //直到第m列的i+1行以後元素均為零 for (i = j + 1; i < n; i++) { s = Guass[i][j]; for (m = j; m < n + 1; m++) { Guass[i][m] = Guass[i][m] - Guass[j][m] * s / (Guass[j][j]); } } }//結束for (j=0;j<n;j++) //回代過程（見公式4.1.5） for (i = n - 1; i >= 0; i--) { s = 0; for (j = i + 1; j < n; j++) { s = s + Guass[i][j] * x[j]; } x[i] = (Guass[i][n] - s) / Guass[i][i]; } return x; }//返回值是函數的系數 public static void main(String[] args) { double[] x = {0, 1, 2, 3, 4, 5, 6, 7}; double[] y = {0, 1, 4, 9, 16, 25, 36, 49}; double[] a = MultiLine(x, y, 8, 2); for(int i =0; i <a.length;i++){ System.out.println(a[i]); } } }

輸出：

0.708333333333342
-0.37500000000000583
1.0416666666666674

取整就得到y=x^2。