Gauss-Newton算法代码详细解释(转载+自己注释)

这篇博客是对[1]中不详细的地方进行细节上的阐述,
并且每句代码都加了注释,使得更加容易理解

下面的论述(包括伪代码和算法)特指被最小化的目标函数是MSE的时候
需要注意,如果不是MSE为目标函数,那么下面的第二张截图开始的伪代码都是需要更换的,这里的伪代码仅仅针对目标函数恰好为MSE的时候
本文的代码来自[2],因为[1]中的代码参数名字有点奇怪,所以就不予采用
该代码的运行需要安装opencv,安装流程使用[3]

###################伪代码#######################################

这里解释下:
上面的▽f\triangledown f▽f为啥是:
2∑i=1mri∂ri∂xj2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial x_j}2i=1∑mri∂xj∂ri
因为这里的fff其实是MSE,也就是:

所以求导后就有了上面▽f\triangledown f▽f的模样

################具体案例与代码解释#############################

例子1，根据美国1815年至1885年数据，估计人口模型中的参数A和B。如下表所示，已知年份和人口总量，及人口模型方程，求方程中的参数(为什么Gauss-Newton用来求解参数在本文后面有解释)。

代码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);
//算法顶层函数声明
double Deriv(double(*Func)(const Mat &input, const Mat params),const Mat &input, const Mat params, int n);
//导数函数声明
// The user defines their function here
double Func(const Mat &input, const Mat params);//下面是被调用的函数//这个函数的作用就是f=A·e^(Bt)
double Func(const Mat &input, const Mat params)
{
//    这里的params样例如下:
//    params=[7.000153882130696;0.262076597545448]
//    也就是说,这里的params类似一个列表.
//所以,其实这里的params=[A,B],input=t// Assumes input is a single row matrix// Assumes params is a column matrixdouble A = params.at<double>(0, 0);//一个浮点数//这里的(0,0)代表获取上述列表的第0个元素double B = params.at<double>(1, 0);//一个浮点数//这里的(0,1)代表获取上述列表的第1个元素double x = input.at<double>(0, 0);//x就是上面函数表达式中的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;//对A,B两个系数缩减params2.at<double>(n, 0) += DERIV_STEP;//对A,B两个系数递增double p1 = Func(input, params1);double p2 = Func(input, params2);double d = (p2 - p1) / (2 * DERIV_STEP);//这里是在计算导数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条数据int n = inputs.cols;//列n=1int num_params = params.rows;//nfum_params=2,表示目标函数表达式中的未知参数的个数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;//拟合效果的指标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);//拟合值与实际值之间的差//之所以是(j,0)是因为输出值肯定只有一列mse += r.at<double>(j, 0)*r.at<double>(j, 0);//残差矩阵的平方,这里之所以要平方是根据MSE的定义来的for (int k = 0; k < num_params; k++)//遍历列{Jf.at<double>(j, k) = Deriv(Func, input, params, k);//对第k个元素求偏导//雅各比矩阵中的某个元素是求导值}}mse /= m;//MSE的定义中需要除以整体元素的个数// The difference in mse is very small, so quitif (fabs(mse - last_mse) < 1e-8)//如果MSE不再变化,就认为拟合结束{break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;//计算△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);//-------------------------------下面是采样数据----------------------------------------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;//-------------------------------------上面是采样数据-----------------------------------------// 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);//下面是初始值设置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是(没这个玩意儿还真运行不了):

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

关于这个代码的一个问题:
NewtonGauss算法明明是为了计算最小值而存在的,为什么到了代码里变成了求参数值A和B?
答:
因为这个代码中,我们常见的Jacobian是对于x1,x2求偏导,其实在代码里面是对应于A和B求偏导,
也就是说,代码其实是在求解对于A和B的Jacobian矩阵.
所以这里关注的"目标函数迭代到最小值"其实是MSE,并不是fff
算法的目的是,当A和B为何值时,目标函数的值最小.
这样就符合了Newton-Gauss算法的本意:
求解目标函数MSE的最低值以及最低值对应的变量A,B的具体数值

代码和伪代码的对应关系:

代码	伪代码
num_params=A,B	x1,x2x_1,x_2x1,x2
m	m
Jf.at(j, k)函数在某个点的导数	Jr=αfαxkJ_r=\frac{\alpha f}{\alpha {x_k}}Jr=αxkαf
点坐标(A,B,outputs.at(i,0))	点坐标(x1,x2,f(x1,x2)x_1,x_2,f(x_1,x_2)x1,x2,f(x1,x2))

Reference:
[1]Gauss-Newton算法学习
[2]梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现
[3]opencv4.0.1配合contrib在linux下面安装编译全过程

Gauss-Newton算法代码详细解释(转载+自己注释)相关推荐

吴恩达机器学习神经网络作业1(用已经求好的权重进行手写数字分类) Python实现代码详细解释
整个项目的github:https://github.com/RobinLuoNanjing/MachineLearning_Ng_Python 里面可以下载进行代码实现的数据集题目介绍: In t ...
吴恩达机器学习逻辑回归作业3(手写数字分类) Python实现代码详细解释
整个项目的github:https://github.com/RobinLuoNanjing/MachineLearning_Ng_Python 里面可以下载进行代码实现的数据集题目介绍: In t ...
吴恩达机器学习逻辑回归作业2(芯片预测) Python实现代码详细解释
整个项目的github:https://github.com/RobinLuoNanjing/MachineLearning_Ng_Python 里面可以下载进行代码实现的数据集题目介绍: In t ...
Apollo Planning决策规划算法代码详细解析 (1)：Scenario选择
本文重点讲解Apollo代码中怎样配置Scenario以及选择当前Scenario,Scenario场景决策是Apollo规划算法的第一步,本文会对代码进行详细解析,也会梳理整个决策流程,码字不易,喜 ...
PointNet代码详细解释(Pytorch版本)
pointnet.pytorch的代码详细解释 1. PointNet的Pytorch版本代码解析链接 2. 代码解释 2.1 代码结构思维导图 2.2 代码注释 2.2.1 build.sh 2.2 ...
关于 IMCRA+OMLSA 语音降噪算法的详细解释
关于 IMCRA+OMLSA 语音降噪算法的详细解释概述 OMLSA 算法 IMCRA 算法概述 IMCRA+OMLSA 算法及其一些变体是目前语音降噪中常用的算法.很多文献在解释这两种算法的时候 ...
SORT跟踪算法的详细解释，不容错过
转载自:https://blog.csdn.net/HaoBBNuanMM/article/details/85555547 SORT - SIMPLE ONLINE AND REALTIME TRA ...
din算法代码_DIN算法代码详细解读
首先给出论文的地址:Deep Interest Network for Click-Through Rate Prediction 然后给出两篇对论文进行了详细介绍的文章: 王喆:推荐系统中的注意力机 ...
Apollo Planning决策规划算法代码详细解析 (5)：规划算法流程介绍
之前的章节介绍了planning模块的整体框架,经过scenario与stage的选择,便进入了具体的task任务,由一系列配置好的task组成了具体的规划算法,本章以apollo中的PublicRo ...

Gauss-Newton算法代码详细解释(转载+自己注释)

Gauss-Newton算法代码详细解释(转载+自己注释)相关推荐

最新文章

热门文章