【转】Gauss-Newton优化方法含代码

本文主要是介绍【转】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，则牛顿迭代公式为：

注：高斯牛顿舍弃了黑森矩阵的计算，用H=JtJl来实现。用一阶导数信息逼近二阶信息项

另外：Hdx=B; dx=H.inverse()*B;H为黑森矩阵，B是什么？B是结果项，观测项。dx求出的为最优方向。或者最优梯度。

看到这里大家都明白高斯牛顿和牛顿法的差异了吧，就在这迭代项上。经典高斯牛顿算法迭代步长λ为1.

那回过头来，高斯牛顿法里为啥要舍弃黑森矩阵的二阶偏导数呢？主要问题是因为牛顿法中Hessian矩阵中的二阶信息项通常难以计算或者花费的工作量很大，而利用整个H的割线近似也不可取，因为在计算梯度时已经得到J(x)，这样H中的一阶信息项JTJ几乎是现成的。鉴于此，为了简化计算，获得有效算法，我们可用一阶导数信息逼近二阶信息项。注意这么干的前提是，残差r接近于零或者接近线性函数从而接近与零时，二阶信息项才可以忽略。通常称为“小残量问题”，否则高斯牛顿法不收敛。

3. 举例

接下来的代码里并没有保证算法收敛的机制，在例子2的自嗨中可以看到劣势。关于自变量维数，代码可以支持多元，但两个例子都是一维的，比如例子1中只有年份t，其实可以增加其他因素的，不必在意。

例子1，根据美国1815年至1885年数据，估计人口模型中的参数A和B。如下表所示，已知年份和人口总量，及人口模型方程，求方程中的参数。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

// A simple demo of Gauss-Newton algorithm on a user defined function

#include <cstdio>

#include <vector>

#include <opencv2/core/core.hpp>

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 ¶ms), // function pointer

const Mat &inputs, const Mat &outputs, Mat ¶ms);

double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer

const Mat &input, const Mat ¶ms, int n);

// The user defines their function here

double Func(const Mat &input, const Mat ¶ms);

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 B

int num_params = 2;

// Generate random data using these parameters

int total_data = 8;

Mat inputs(total_data, 1, CV_64F);

Mat outputs(total_data, 1, CV_64F);

//load observation data

for(int i=0; i < total_data; i++) {

inputs.at<double>(i,0) = i+1; //load year

}

//load America population

outputs.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);

//init guess

params.at<double>(0,0) = 6;

params.at<double>(1,0) = 0.3;

GaussNewton(Func, inputs, outputs, params);

printf("Parameters from GaussNewton: %f %f\n", params.at<double>(0,0), params.at<double>(1,0));

return 0;

}

double Func(const Mat &input, const Mat ¶ms)

{

// Assumes input is a single row matrix

// Assumes params is a column matrix

double A = params.at<double>(0,0);

double B = params.at<double>(1,0);

double x = input.at<double>(0,0);

return 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 ¶ms), const Mat &input, const Mat ¶ms, int n)

{

// Assumes input is a single row matrix

// Returns the derivative of the nth parameter

Mat params1 = params.clone();

Mat params2 = params.clone();

// Use central difference to get derivative

params1.at<double>(n,0) -= DERIV_STEP;

params2.at<double>(n,0) += DERIV_STEP;

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 ¶ms),

const Mat &inputs, const Mat &outputs, Mat ¶ms)

{

int m = inputs.rows;

int n = inputs.cols;

int num_params = params.rows;

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 input

double last_mse = 0;

for(int i=0; i < MAX_ITER; i++) {

double mse = 0;

for(int j=0; j < m; j++) {

for(int k=0; k < n; k++) {//copy Independent variable vector, the year

input.at<double>(0,k) = inputs.at<double>(j,k);

}

r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);//diff between estimate and observation population

mse += r.at<double>(j,0)*r.at<double>(j,0);

for(int k=0; k < num_params; k++) {

Jf.at<double>(j,k) = Deriv(Func, input, params, k);

}

mse /= m;

// The difference in mse is very small, so quit

if(fabs(mse - last_mse) < 1e-8) {

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;

}

　　运行结果：

A=7.0,B=0.26 (初始值，A=6,B=0.3)，100次迭代到第4次就收敛了。

若初始值A=1,B=1，则要迭代14次收敛。

下图为根据上面得到的A、B系数，利用matlab拟合的人口模型曲线

例子2：我想要拟合如下模型，

由于缺乏观测量，就自导自演，假设4个参数已知A=5,B=1,C=10,D=2，构造100个随机数作为x的观测值，计算相应的函数观测值。然后，利用这些观测值，反推4个参数。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

// A simple demo of Gauss-Newton algorithm on a user defined function

#include <cstdio>

#include <vector>

#include <opencv2/core/core.hpp>

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 ¶ms), // function pointer

const Mat &inputs, const Mat &outputs, Mat ¶ms);

double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer

const Mat &input, const Mat ¶ms, int n);

// The user defines their function here

double Func(const Mat &input, const Mat ¶ms);

int main()

{

// For this demo we're going to try and fit to the function

// F = A*sin(Bx) + C*cos(Dx)

// There are 4 parameters: A, B, C, D

int num_params = 4;

// Generate random data using these parameters

int total_data = 100;

double A = 5;

double B = 1;

double C = 10;

double D = 2;

Mat inputs(total_data, 1, CV_64F);

Mat outputs(total_data, 1, CV_64F);

for(int i=0; i < total_data; i++) {

double x = -10.0 + 20.0* rand() / (1.0 + RAND_MAX); // random between [-10 and 10]

double y = A*sin(B*x) + C*cos(D*x);

// Add some noise

// y += -1.0 + 2.0*rand() / (1.0 + RAND_MAX);

inputs.at<double>(i,0) = x;

outputs.at<double>(i,0) = y;

}

// 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) = 1;

params.at<double>(1,0) = 1;

params.at<double>(2,0) = 8; // changing to 1 will cause it not to find the solution, too far away

params.at<double>(3,0) = 1;

GaussNewton(Func, inputs, outputs, params);

printf("True parameters: %f %f %f %f\n", A, B, C, D);

printf("Parameters from GaussNewton: %f %f %f %f\n", params.at<double>(0,0), params.at<double>(1,0),

params.at<double>(2,0), params.at<double>(3,0));

return 0;

}

double Func(const Mat &input, const Mat ¶ms)

{

// Assumes input is a single row matrix

// Assumes params is a column matrix

double A = params.at<double>(0,0);

double B = params.at<double>(1,0);

double C = params.at<double>(2,0);

double D = params.at<double>(3,0);

double x = input.at<double>(0,0);

return A*sin(B*x) + C*cos(D*x);

}

double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), const Mat &input, const Mat ¶ms, int n)

{

// Assumes input is a single row matrix

// Returns the derivative of the nth parameter

Mat params1 = params.clone();

Mat params2 = params.clone();

// Use central difference to get derivative

params1.at<double>(n,0) -= DERIV_STEP;

params2.at<double>(n,0) += DERIV_STEP;

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 ¶ms),

const Mat &inputs, const Mat &outputs, Mat ¶ms)

{

int m = inputs.rows;

int n = inputs.cols;

int num_params = params.rows;

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 input

double last_mse = 0;

for(int i=0; i < MAX_ITER; i++) {

double mse = 0;

for(int j=0; j < m; j++) {

for(int k=0; k < n; k++) {

input.at<double>(0,k) = inputs.at<double>(j,k);

}

r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);

mse += r.at<double>(j,0)*r.at<double>(j,0);

for(int k=0; k < num_params; k++) {

Jf.at<double>(j,k) = Deriv(Func, input, params, k);

}

mse /= m;

// The difference in mse is very small, so quit

if(fabs(mse - last_mse) < 1e-8) {

break;

}

Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;

params += delta;

//printf("%d: mse=%f\n", i, mse);

printf("%f\n",mse);

last_mse = mse;

}

　　运行结果，得到的参数并不够理想，50次后收敛了

下图中，每次迭代残差并没有持续减少，有反复

4．优缺点分析

优点：

对于零残量问题，即r=0，有局部二阶收敛速度

对于小残量问题，即r较小或接近线性，有快的局部收敛速度

对于线性最小二乘问题，一步达到极小点

缺点：

对于不是很严重的大残量问题，有较慢的局部收敛速度

对于残量很大的问题或r的非线性程度很大的问题，不收敛

不一定总体收敛

如果J不满秩，则方法无定义

对于它的缺点，我们通过增加线性搜索策略，保证目标函数每一步下降，对于几乎所有非线性最小二乘问题，它都具有局部收敛性及总体收敛，即所谓的阻尼高斯牛顿法。

这篇关于【转】Gauss-Newton优化方法含代码的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

【转】Gauss-Newton优化方法含代码

相关文章

Vue3 的 shallowRef 和 shallowReactive：优化性能

HDFS—存储优化（纠删码）

使用opencv优化图片（画面变清晰）

【C++】_list常用方法解析及模拟实现

活用c4d官方开发文档查询代码

浅谈主机加固，六种有效的主机加固方法

webm怎么转换成mp4？这几种方法超多人在用！

poj 1258 Agri-Net（最小生成树模板代码）

MySQL高性能优化规范

透彻！驯服大型语言模型（LLMs）的五种方法，及具体方法选择思路

【转】Gauss-Newton优化方法 含代码

相关文章

【转】Gauss-Newton优化方法含代码