這篇博客是對(duì)[1]中不詳細(xì)的地方進(jìn)行細(xì)節(jié)上的闡述,
并且每句代碼都加了注釋,使得更加容易理解

下面的論述(包括偽代碼和算法)特指被最小化的目標(biāo)函數(shù)是MSE的時(shí)候
需要注意,如果不是MSE為目標(biāo)函數(shù),那么下面的第二張截圖開(kāi)始的偽代碼都是需要更換的,這里的偽代碼僅僅針對(duì)目標(biāo)函數(shù)恰好為MSE的時(shí)候
本文的代碼來(lái)自[2],因?yàn)閇1]中的代碼參數(shù)名字有點(diǎn)奇怪,所以就不予采用
該代碼的運(yùn)行需要安裝opencv,安裝流程使用[3]

###################偽代碼#######################################

這里解釋下:
上面的$▽f$為啥是:
$$2\sum _{i=1}^m{r}_i\frac{?r_i}{?x_j}$$
因?yàn)檫@里的$f$其實(shí)是MSE,也就是:

所以求導(dǎo)后就有了上面$▽f$的模樣

################具體案例與代碼解釋#############################

例子1，根據(jù)美國(guó)1815年至1885年數(shù)據(jù)，估計(jì)人口模型中的參數(shù)A和B。如下表所示，已知年份和人口總量，及人口模型方程，求方程中的參數(shù)(為什么Gauss-Newton用來(lái)求解參數(shù)在本文后面有解釋)。

代碼main.cpp如下:

#include <cstdio> #include <vector> #include <opencv2/core.hpp> #include<iostream> using namespace std; using namespace cv;const double DERIV_STEP = 1e-5; const int MAX_ITER = 100;void GaussNewton(double(*Func)(const Mat &input, const Mat params),const Mat &inputs, const Mat &outputs, Mat params); //算法頂層函數(shù)聲明 double Deriv(double(*Func)(const Mat &input, const Mat params),const Mat &input, const Mat params, int n); //導(dǎo)數(shù)函數(shù)聲明 // The user defines their function here double Func(const Mat &input, const Mat params);//下面是被調(diào)用的函數(shù)//這個(gè)函數(shù)的作用就是f=A·e^(Bt) double Func(const Mat &input, const Mat params) { // 這里的params樣例如下: // params=[7.000153882130696;0.262076597545448] // 也就是說(shuō),這里的params類似一個(gè)列表. //所以,其實(shí)這里的params=[A,B],input=t// Assumes input is a single row matrix// Assumes params is a column matrixdouble A = params.at<double>(0, 0);//一個(gè)浮點(diǎn)數(shù)//這里的(0,0)代表獲取上述列表的第0個(gè)元素double B = params.at<double>(1, 0);//一個(gè)浮點(diǎn)數(shù)//這里的(0,1)代表獲取上述列表的第1個(gè)元素double x = input.at<double>(0, 0);//x就是上面函數(shù)表達(dá)式中的treturn A*exp(x*B); }//calc the n-th params' partial derivation ， the params are our final target double Deriv(double(*Func)(const Mat &input, const Mat params), const Mat &input, const Mat params, int n) {// Assumes input is a single row matrix// Returns the derivative of the nth parameterMat params1 = params.clone();Mat params2 = params.clone();// Use central difference to get derivativeparams1.at<double>(n, 0) -= DERIV_STEP;//對(duì)A,B兩個(gè)系數(shù)縮減params2.at<double>(n, 0) += DERIV_STEP;//對(duì)A,B兩個(gè)系數(shù)遞增double p1 = Func(input, params1);double p2 = Func(input, params2);double d = (p2 - p1) / (2 * DERIV_STEP);//這里是在計(jì)算導(dǎo)數(shù)return d; }void GaussNewton(double(*Func)(const Mat &input, const Mat params),const Mat &inputs, const Mat &outputs, Mat params) {int m = inputs.rows;//行m=8,表示8條數(shù)據(jù)int n = inputs.cols;//列n=1int num_params = params.rows;//nfum_params=2,表示目標(biāo)函數(shù)表達(dá)式中的未知參數(shù)的個(gè)數(shù)Mat r(m, 1, CV_64F); // residual matrix殘差矩陣Mat Jf(m, num_params, CV_64F); // Jacobian of Func()雅各比矩陣Mat input(1, n, CV_64F); // single row inputdouble last_mse = 0;for (int i = 0; i < MAX_ITER; i++){double mse = 0;//擬合效果的指標(biāo)for (int j = 0; j < m; j++)//遍歷行, ?Xj{for (int k = 0; k < n; k++)//遍歷列,?Xk{//copy Independent variable vector, the yearinput.at<double>(0, k) = inputs.at<double>(j, k);}r.at<double>(j, 0) = outputs.at<double>(j, 0) - Func(input, params);//擬合值與實(shí)際值之間的差//之所以是(j,0)是因?yàn)檩敵鲋悼隙ㄖ挥幸涣衜se += r.at<double>(j, 0)*r.at<double>(j, 0);//殘差矩陣的平方,這里之所以要平方是根據(jù)MSE的定義來(lái)的for (int k = 0; k < num_params; k++)//遍歷列{Jf.at<double>(j, k) = Deriv(Func, input, params, k);//對(duì)第k個(gè)元素求偏導(dǎo)//雅各比矩陣中的某個(gè)元素是求導(dǎo)值}}mse /= m;//MSE的定義中需要除以整體元素的個(gè)數(shù)// The difference in mse is very small, so quitif (fabs(mse - last_mse) < 1e-8)//如果MSE不再變化,就認(rèn)為擬合結(jié)束{break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;//計(jì)算△params += delta;//printf("%d: mse=%f\n", i, mse);printf("%d %f\n", i, mse);last_mse = mse;} }//------------------下面是頂層文件----------------------------------------------int main() {// For this demo we're going to try and fit to the function// F = A*exp(t*B)， There are 2 parameters: A Bint num_params = 2;// Generate random data using these parametersint total_data = 8;Mat inputs(total_data, 1, CV_64F);//這里的CV_64F代表一種單通道的矩陣類型Mat outputs(total_data, 1, CV_64F);//-------------------------------下面是采樣數(shù)據(jù)----------------------------------------for (int i = 0; i < total_data; i++){inputs.at<double>(i, 0) = i + 1; //load year}//load America populationoutputs.at<double>(0, 0) = 8.3;outputs.at<double>(1, 0) = 11.0;outputs.at<double>(2, 0) = 14.7;outputs.at<double>(3, 0) = 19.7;outputs.at<double>(4, 0) = 26.7;outputs.at<double>(5, 0) = 35.2;outputs.at<double>(6, 0) = 44.4;outputs.at<double>(7, 0) = 55.9;//-------------------------------------上面是采樣數(shù)據(jù)-----------------------------------------// Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!Mat params(num_params, 1, CV_64F);//下面是初始值設(shè)置params.at<double>(0, 0) = 6;params.at<double>(1, 0) = 0.3;GaussNewton(Func, inputs, outputs, params);cout<<"最終params="<<params<<endl;printf("Parameters from GaussNewton: %f %f\n", params.at<double>(0, 0), params.at<double>(1, 0));return 0; }

Clion2018.3中的CMakeLists.txt是(沒(méi)這個(gè)玩意兒還真運(yùn)行不了):

cmake_minimum_required(VERSION 3.12) project(GaussNewton)set(CMAKE_CXX_STANDARD 14) include_directories($ENV{CMAKE_INCLUDE_PATH}) set(CMAKE_CXX_STANDARD 14)#C++ standard set(OpenCV_DIR /home/appleyuchi/opencv/opencv_install/lib/cmake/opencv4) find_package( OpenCV REQUIRED ) # locate OpenCV in system include_directories( ${OpenCV_INCLUDE_DIRS} ) # provide library headers add_executable(GaussNewton main.cpp) target_link_libraries(GaussNewton ${OpenCV_LIBS} /home/appleyuchi/opencv/opencv_install/lib/libopencv_highgui.so) # link OpenCV libraries , hightgui.so not found by cmake so this hack MESSAGE("OpenCV_LIBS: " ${OpenCV_LIBS} ) #display opencv libs found

使用的opencv版本是4.0.1

關(guān)于這個(gè)代碼的一個(gè)問(wèn)題:
NewtonGauss算法明明是為了計(jì)算最小值而存在的,為什么到了代碼里變成了求參數(shù)值A(chǔ)和B?
答:
因?yàn)檫@個(gè)代碼中,我們常見(jiàn)的Jacobian是對(duì)于x1,x2求偏導(dǎo),其實(shí)在代碼里面是對(duì)應(yīng)于A和B求偏導(dǎo),
也就是說(shuō),代碼其實(shí)是在求解對(duì)于A和B的Jacobian矩陣.
所以這里關(guān)注的"目標(biāo)函數(shù)迭代到最小值"其實(shí)是MSE,并不是$f$
算法的目的是,當(dāng)A和B為何值時(shí),目標(biāo)函數(shù)的值最小.
這樣就符合了Newton-Gauss算法的本意:
求解目標(biāo)函數(shù)MSE的最低值以及最低值對(duì)應(yīng)的變量A,B的具體數(shù)值

代碼和偽代碼的對(duì)應(yīng)關(guān)系:

代碼偽代碼

num_params=A,B	$x_1,x_2$
m	m
Jf.at(j, k)函數(shù)在某個(gè)點(diǎn)的導(dǎo)數(shù)	$J_r=\frac{\alpha f}{\alpha x_k}$
點(diǎn)坐標(biāo)(A,B,outputs.at(i,0))	點(diǎn)坐標(biāo)($x_1,x_2,f(x_1,x_2)$)

Reference:
[1]Gauss-Newton算法學(xué)習(xí)
[2]梯度下降法，牛頓法，高斯-牛頓迭代法，附代碼實(shí)現(xiàn)
[3]opencv4.0.1配合contrib在linux下面安裝編譯全過(guò)程

總結(jié)

以上是生活随笔為你收集整理的Gauss-Newton算法代码详细解释(转载+自己注释)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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