1. 程式人生 > >數值優化之高斯-牛頓法(Gauss-Newton)

數值優化之高斯-牛頓法(Gauss-Newton)

一、基本概念定義

1.非線性方程定義及最優化方法簡述

   指因變數與自變數之間的關係不是線性的關係,比如平方關係、對數關係、指數關係、三角函式關係等等。對於此類方程,求解n元實函式f在整個n維向量空間Rn上的最優值點往往很難得到精確解,經常需要求近似解問題。

   求解該最優化問題的方法大多是逐次一維搜尋的迭代演算法,基本思想是在一個近似點處選定一個有利於搜尋方向,沿這個方向進行一維搜尋,得到新的近似點。如此反覆迭代,知道滿足預定的精度要求為止。根據搜尋方向的取法不同,這類迭代演算法可分為兩類:

解析法:需要用目標函式的到函式,

梯度法:又稱最速下降法,是早期的解析法,收斂速度較慢

牛頓法:收斂速度快,但不穩定,計算也較困難。高斯牛頓法基於其改進,但目標作用不同

共軛梯度法:收斂較快,效果好

變尺度法:效率較高,常用DFP法(Davidon Fletcher Powell)

直接法:不涉及導數,只用到函式值。有交替方向法(又稱座標輪換法)、模式搜尋法、旋轉方向法、鮑威爾共軛方向法和單純形加速法等。


2.非線性最小二乘問題

   非線性最小二乘問題來自於非線性迴歸,即通過觀察自變數和因變數資料,求非線性目標函式的係數引數,使得函式模型與觀測量儘量相似。

   高斯牛頓法解決非線性最小二乘問題的最基本方法,並且它只能處理二次函式。(使用時必須將目標函式轉化為二次的)

   Unlike Newton'smethod, the Gauss–Newton algorithm can only be used to minimize a sum ofsquared function values


3.基本數學表達

a.梯度gradient,由多元函式的各個偏導陣列成的向量

以二元函式為例,其梯度為:


b.黑森矩陣Hessian matrix,由多元函式的二階偏導陣列成的方陣,描述函式的局部曲率,以二元函式為例,


c.雅可比矩陣 Jacobian matrix,是多元函式一階偏導數以一定方式排列成的矩陣,體現了一個可微方程與給出點的最優線性逼近。以二元函式為例,


如果擴充套件多維的話F: Rn-> Rm,則雅可比矩陣是一個m行n列的矩陣:


雅可比矩陣作用,如果P是Rn中的一點,F在P點可微分,那麼在這一點的導數由JF(P)給出,在此情況下,由F(P)描述的線性運算元即接近點P的F的最優線性逼近:


d.殘差 residual,表示實際觀測值與估計值(擬合值)之間的差


二、牛頓法

牛頓法的基本思想是採用多項式函式來逼近給定的函式值,然後求出極小點的估計值,重複操作,直到達到一定精度為止。

1.考慮如下一維無約束的極小化問題:


因此,一維牛頓法的計算步驟如下:



需要注意的是,牛頓法在求極值的時候,如果初始點選取不好,則可能不收斂於極小點


2.下面給出多維無約束極值的情形:

若非線性目標函式f(x)具有二階連續偏導,在x(k)為其極小點的某一近似,在這一點取f(x)的二階泰勒展開,即:


  如果f(x)是二次函式,則其黑森矩陣H為常數,式(1)是精確的(等於號),在這種情況下,從任意一點處罰,用式(2)只要一步可求出f(x)的極小點(假設黑森矩陣正定,所有特徵值大於0)

  如果f(x)不是二次函式,式(1)僅是一個近似表示式,此時,按式(2)求得的極小點,只是f(x)的近似極小點。在這種情況下,常按照下面選取搜尋方向:


牛頓法收斂的速度很快,當f(x)的二階導數及其黑森矩陣的逆矩陣便於計算時,這一方法非常有效。【但通常黑森矩陣很不好求】

3.下面給出一個實際計算例子。


例:試用牛頓法求的極小值

解:



【f(x)是二次函式,H矩陣為常數,只要任意點出發,只要一步即可求出極小點】


三、牛頓高斯法

1.      gauss-newton是如何由上述派生的

有時候為了擬合數據,比如根據重投影誤差求相機位姿(R,T為方程係數),常常將求解模型轉化為非線性最小二乘問題。高斯牛頓法正是用於解決非線性最小二乘問題,達到資料擬合、引數估計和函式估計的目的。

假設我們研究如下形式的非線性最小二乘問題:



這兩個位置間殘差(重投影誤差):

如果有大量觀測點(多維),我們可以通過選擇合理的T使得殘差的平方和最小求得兩個相機之間的位姿。機器視覺這塊暫時不擴充套件,接著說怎麼求非線性最小二乘問題。

若用牛頓法求式3,則牛頓迭代公式為:


看到這裡大家都明白高斯牛頓和牛頓法的差異了吧,就在這迭代項上。經典高斯牛頓演算法迭代步長λ為1.

那回過頭來,高斯牛頓法裡為啥要捨棄黑森矩陣的二階偏導數呢?主要問題是因為牛頓法中Hessian矩陣中的二階資訊項通常難以計算或者花費的工作量很大,而利用整個H的割線近似也不可取,因為在計算梯度時已經得到J(x),這樣H中的一階資訊項JTJ幾乎是現成的。鑑於此,為了簡化計算,獲得有效演算法,我們可用一階導數資訊逼近二階資訊項。注意這麼幹的前提是,殘差r接近於零或者接近線性函式從而接近與零時,二階資訊項才可以忽略。通常稱為小殘量問題,否則高斯牛頓法不收斂。


3.  舉例

接下來的程式碼裡並沒有保證演算法收斂的機制,在例子2的自嗨中可以看到劣勢。關於自變數維數,程式碼可以支援多元,但兩個例子都是一維的,比如例子1中只有年份t,其實可以增加其他因素的,不必在意。

例子1,根據美國1815年至1885年資料,估計人口模型中的引數A和B。如下表所示,已知年份和人口總量,及人口模型方程,求方程中的引數。



  1. // A simple demo of Gauss-Newton algorithm on a user defined function
  2. #include <cstdio>
  3. #include <vector>
  4. #include <opencv2/core/core.hpp>
  5. usingnamespace std;  
  6. usingnamespace cv;  
  7. constdouble DERIV_STEP = 1e-5;  
  8. constint MAX_ITER = 100;  
  9. void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer
  10.                  const Mat &inputs, const Mat &outputs, Mat ¶ms);  
  11. double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer
  12.              const Mat &input, const Mat ¶ms, int n);  
  13. // The user defines their function here
  14. double Func(const Mat &input, const Mat ¶ms);  
  15. int main()  
  16. {  
  17.     // For this demo we're going to try and fit to the function
  18.     // F = A*exp(t*B)
  19.     // There are 2 parameters: A B
  20.     int num_params = 2;  
  21.     // Generate random data using these parameters
  22.     int total_data = 8;  
  23.     Mat inputs(total_data, 1, CV_64F);  
  24.     Mat outputs(total_data, 1, CV_64F);  
  25.     //load observation data
  26.     for(int i=0; i < total_data; i++) {  
  27.         inputs.at<double>(i,0) = i+1;  //load year
  28.     }  
  29.     //load America population
  30.     outputs.at<double>(0,0)= 8.3;  
  31.     outputs.at<double>(1,0)= 11.0;  
  32.     outputs.at<double>(2,0)= 14.7;  
  33.     outputs.at<double>(3,0)= 19.7;  
  34.     outputs.at<double>(4,0)= 26.7;  
  35.     outputs.at<double>(5,0)= 35.2;  
  36.     outputs.at<double>(6,0)= 44.4;  
  37.     outputs.at<double>(7,0)= 55.9;  
  38.     // Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!
  39.     Mat params(num_params, 1, CV_64F);  
  40.     //init guess
  41.     params.at<double>(0,0) = 6;  
  42.     params.at<double>(1,0) = 0.3;  
  43.     GaussNewton(Func, inputs, outputs, params);  
  44.     printf("Parameters from GaussNewton: %f %f\n", params.at<double>(0,0), params.at<double>(1,0));  
  45.     return 0;  
  46. }  
  47. double Func(const Mat &input, const Mat ¶ms)  
  48. {  
  49.     // Assumes input is a single row matrix
  50.     // Assumes params is a column matrix
  51.     double A = params.at<double>(0,0);  
  52.     double B = params.at<double>(1,0);  
  53.     double x = input.at<double>(0,0);  
  54.     return A*exp(x*B);  
  55. }  
  56. //calc the n-th params' partial derivation , the params are our  final target
  57. double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), const Mat &input, const Mat ¶ms, int n)  
  58. {  
  59.     // Assumes input is a single row matrix
  60.     // Returns the derivative of the nth parameter
  61.     Mat params1 = params.clone();  
  62.     Mat params2 = params.clone();  
  63.     // Use central difference  to get derivative
  64.     params1.at<double>(n,0) -= DERIV_STEP;  
  65.     params2.at<double>(n,0) += DERIV_STEP;  
  66.     double p1 = Func(input, params1);  
  67.     double p2 = Func(input, params2);  
  68.     double d = (p2 - p1) / (2*DERIV_STEP);  
  69.     return d;  
  70. }  
  71. void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms),  
  72.                  const Mat &inputs, const Mat &outputs, Mat ¶ms)  
  73. {  
  74.     int m = inputs.rows;  
  75.     int n = inputs.cols;  
  76.     int num_params = params.rows;  
  77.     Mat r(m, 1, CV_64F); // residual matrix
  78.     Mat Jf(m, num_params, CV_64F); // Jacobian of Func()
  79. 相關推薦

    數值優化-牛頓Gauss-Newton

    一、基本概念定義1.非線性方程定義及最優化方法簡述   指因變數與自變數之間的關係不是線性的關係,比如平方關係、對數關係、指數關係、三角函式關係等等。對於此類方程,求解n元實函式f在整個n維向量空間Rn上的最優值點往往很難得到精確解,經常需要求近似解問題。   求解該最優化問

    聚類混合模型Gaussian Mixture Model

    k-means應該是原來級別的聚類方法了,這整理下一個使用後驗概率準確評測其精度的方法—高斯混合模型。 我們談到了用 k-means 進行聚類的方法,這次我們來說一下另一個很流行的演算法:Gaussian Mixture Model (GMM)。事實上,GMM

    聚類混合模型Gaussian Mixture Model【轉】

    k-means應該是原來級別的聚類方法了,這整理下一個使用後驗概率準確評測其精度的方法—高斯混合模型。 我們談到了用 k-means 進行聚類的方法,這次我們來說一下另一個很流行的演算法:Gaussian Mixture Model (GMM)。事實上,GMM 和 k-means 很像,不過 GMM 是學習

    迭代求解最優化問題——最小二乘問題、牛頓

    最小二乘問題 最小二乘問題是應用最廣泛的優化問題,它的一般形式如下: minx||r(x)||2 該問題的損失函式為S(x)=||r(x)||2。其中r(x)為殘差函式,一般表示預測值與實際值的差別。一個最簡單的最小二乘問題就是線性迴歸問題,對於這個問題的

    EM演算法混合模型

    單個高斯模型 如果我們有一堆資料,其分佈屬於一個高斯模型,那麼有 p(X)=N(x|μ,Σ)=1(2π)m|Σ|‾‾‾‾‾‾‾‾√exp[−12(x−μ)TΣ−1(x−μ)](1.1) p(X) = N(x|\mu,\Sigma) = \

    混合模型GMM model以及梯度下降gradient descent更新引數

    關於GMM模型的資料和 EM 引數估算的資料,網上已經有很多了,今天想談的是GMM的協方差矩陣的分析、GMM的引數更新方法 1、GMM協方差矩陣的物理含義 涉及到每個元素,是這樣求算: 用中文來描述就是: 注意後面的那個除以(樣本數-1),就是大括號外面的E求期望 (這叫

    梯度下降牛頓牛頓、LM最優化算

    src tro 分享 image 最優化 ima str img 圖片 1、梯度下降法 2、牛頓法 3、高斯牛頓法 4、LM算法 梯度下降法、牛頓法、高斯牛頓法、LM最優化算法

    樸素貝葉Naive Bayes

    ive log 分布 做了 規模 line clas 獨立 輸入數據 1. 前言 說到樸素貝葉斯算法,首先牽扯到的一個概念是判別式和生成式。 判別式:就是直接學習出特征輸出\(Y\)和特征\(X\)之間的關系,如決策函數\(Y=f(X)\),或者從概率論的角度,求出條件分

    Spark MLlib水塘抽樣算Reservoir Sampling

    抽樣 返回 算法 蓄水池抽樣 seq pack param long nds 1.理解  問題定義可以簡化如下:在不知道文件總行數的情況下,如何從文件中隨機的抽取一行?   首先想到的是我們做過類似的題目嗎?當然,在知道文件行數的情況下,我們可以很容易的用C運行庫的rand

    廣義線性迴歸邏輯諦迴歸 Logistic Regression

    廣義線性模型 邏輯斯諦迴歸概念可以認為是屬於廣義線性迴歸的範疇,但它是用來進行分類的。 線性模型的表示式為: f (

    排序演算法希爾排序c#實現

    希爾排序演算法是將陣列的所有元素按照一定增量d分組,對每組內的資料實行插入排序,之後不斷減小增量,每組內所包含的元素也越多,增量減少至1時,整個陣列被分成一組,插入排序結束後整個陣列的排序便完成。演算法流程圖:操作步驟:初始時,有一個大小為 10 的無序序列。(1)在第一趟排

    Android本地圖片或者網路圖片模糊效果毛玻璃效果圖片模糊效果一行程式碼搞定

    一,實現本地圖片或者網路圖片的毛玻璃效果特別方便,只需要把下面的FastBlurUtil類複製到你的專案中就行 package com.testdemo.blur_image_lib10;   import android.graphics.Bitmap;   import andr

    hdu3364-消元取模

    列n個方程,表示每個燈會被哪些開關控制, 得到一個一個n*m的矩陣,最後一列為所要求的狀態 則對這個(n+1)*(m)的矩陣高斯消元,得到方案數,答案爆int #include <cstdio> #include <cmath> #inclu

    優化演算法牛頓

    一、牛頓法 上述描述的都是隻有一個自變數X的一元情況,如果是多元的,比如x1,x2,x3...,xn 呢? 二、對比分析梯度下降演算法 從本質上去看,牛頓法是二階收斂,梯度下降是一階收斂,所以牛頓法就更快。如果更通俗地說的話,比如你想找一條最短的路徑走到一個盆地的最底部,梯度下降

    復習——消元ssoi

    ring 模板 con 這樣的 using pos 但是 -1 stream 題目: 題目描述 Tom 是個品學兼優的好學生,但由於智商問題,算術學得不是很好,尤其是在解方程這個方面。雖然他解決 2x=2 這樣的方程遊刃有余,但是對於下面這樣的方程組就束手無策了。x+y=

    優化演算法:牛頓Newton

    學習深度學習時遇到二階優化演算法牛頓法,查閱了相關書籍進行記錄。 :函式的梯度向量 :函式的Hessian矩陣,其第i行第j列的元素為. 假設是二階連續可微函式,。最速下降法因為迭代路線呈鋸齒形,固收斂速度慢,僅是線性的。最速下降法本質使用線性函式去近似目標函式。要得到快速的演算法,

    順序消元Python實現

    main python實現 ber seq rev div 順序 inf break # coding: utf8 import numpy as np # 設置矩陣 def getInput(): matrix_a = np.mat([[2, 3, 11,

    UnityShader例項14:螢幕特效模糊Gaussian Blur

    Shader "PengLu/ImageEffect/Unlit/GaussianBlur" { Properties { _MainTex ("Base (RGB)", 2D) = "white" {} } CGINCLUDE #include "UnityCG.cginc" s

    梯度下降牛頓-牛頓迭代,附程式碼實現

    ---------------------梯度下降法------------------- 梯度的一般解釋: f(x)在x0的梯度:就是f(x)變化最快的方向。梯度下降法是一個最優化演算法,通常也稱為最速下降法。 假設f(x)是一座山,站在半山腰,往x方向走1米,高度上

    使用牛頓確定邏輯諦迴歸Logistic Regression最佳迴歸係數

    邏輯斯諦迴歸 在邏輯斯諦迴歸中,因為使用梯度上升(gradient ascent)收斂較慢,固本文采用牛頓法(Newton’s Method)進行引數求解,試驗發現通常迭代10次左右就可達到收斂,而梯度上升法則需要迭代上百甚至上千次,當然實際的迭代次數也要視實際資料而定