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

2023-10-31 08:40
文章标签 代码 方法 优化 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。如下表所示，已知年份和人口总量，及人口模型方程，求方程中的参数。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

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个参数。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

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优化方法 含代码的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/314136。

相关文章

Nginx设置连接超时并进行测试的方法步骤

Nginx设置连接超时并进行测试的方法步骤

《Nginx设置连接超时并进行测试的方法步骤》在高并发场景下,如果客户端与服务器的连接长时间未响应,会占用大量的系统资源,影响其他正常请求的处理效率,为了解决这个问题,可以通过设置Nginx的连接... 目录设置连接超时目的操作步骤测试连接超时测试方法：总结：设置连接超时目的设置客户端与服务器之间的连接
阅读更多...
Java判断多个时间段是否重合的方法小结

Java判断多个时间段是否重合的方法小结

《Java判断多个时间段是否重合的方法小结》这篇文章主要为大家详细介绍了Java中判断多个时间段是否重合的方法,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录判断多个时间段是否有间隔判断时间段集合是否与某时间段重合判断多个时间段是否有间隔实体类内容public class D
阅读更多...
Python使用国内镜像加速pip安装的方法讲解

Python使用国内镜像加速pip安装的方法讲解

《Python使用国内镜像加速pip安装的方法讲解》在Python开发中,pip是一个非常重要的工具,用于安装和管理Python的第三方库,然而,在国内使用pip安装依赖时,往往会因为网络问题而导致速... 目录一、pip 工具简介1. 什么是 pip？2. 什么是 -i 参数？二、国内镜像源的选择三、如何
阅读更多...
IDEA编译报错“java: 常量字符串过长”的原因及解决方法

IDEA编译报错“java: 常量字符串过长”的原因及解决方法

《IDEA编译报错“java:常量字符串过长”的原因及解决方法》今天在开发过程中,由于尝试将一个文件的Base64字符串设置为常量,结果导致IDEA编译的时候出现了如下报错java:常量字符串过长,... 目录一、问题描述二、问题原因2.1 理论角度2.2 源码角度三、解决方案解决方案①：StringBui
阅读更多...
Linux使用nload监控网络流量的方法

Linux使用nload监控网络流量的方法

《Linux使用nload监控网络流量的方法》Linux中的nload命令是一个用于实时监控网络流量的工具,它提供了传入和传出流量的可视化表示,帮助用户一目了然地了解网络活动,本文给大家介绍了Linu... 目录简介安装示例用法基础用法指定网络接口限制显示特定流量类型指定刷新率设置流量速率的显示单位监控多个
阅读更多...
Java覆盖第三方jar包中的某一个类的实现方法

Java覆盖第三方jar包中的某一个类的实现方法

《Java覆盖第三方jar包中的某一个类的实现方法》在我们日常的开发中,经常需要使用第三方的jar包,有时候我们会发现第三方的jar包中的某一个类有问题,或者我们需要定制化修改其中的逻辑,那么应该如何... 目录一、需求描述二、示例描述三、操作步骤四、验证结果五、实现原理一、需求描述需求描述如下：需要在
阅读更多...
JavaScript中的reduce方法执行过程、使用场景及进阶用法

JavaScript中的reduce方法执行过程、使用场景及进阶用法

《JavaScript中的reduce方法执行过程、使用场景及进阶用法》：本文主要介绍JavaScript中的reduce方法执行过程、使用场景及进阶用法的相关资料,reduce是JavaScri... 目录1. 什么是reduce2. reduce语法2.1 语法2.2 参数说明3. reduce执行过程
阅读更多...
C#中读取XML文件的四种常用方法

C#中读取XML文件的四种常用方法

《C#中读取XML文件的四种常用方法》Xml是Internet环境中跨平台的,依赖于内容的技术,是当前处理结构化文档信息的有力工具,下面我们就来看看C#中读取XML文件的方法都有哪些吧... 目录XML简介格式C#读取XML文件方法使用XmlDocument使用XmlTextReader/XmlTextWr
阅读更多...
Java调用DeepSeek API的最佳实践及详细代码示例

Java调用DeepSeek API的最佳实践及详细代码示例

《Java调用DeepSeekAPI的最佳实践及详细代码示例》：本文主要介绍如何使用Java调用DeepSeekAPI,包括获取API密钥、添加HTTP客户端依赖、创建HTTP请求、处理响应、... 目录1. 获取API密钥2. 添加HTTP客户端依赖3. 创建HTTP请求4. 处理响应5. 错误处理6.
阅读更多...
C++初始化数组的几种常见方法(简单易懂)

C++初始化数组的几种常见方法(简单易懂)

《C++初始化数组的几种常见方法(简单易懂)》本文介绍了C++中数组的初始化方法,包括一维数组和二维数组的初始化,以及用new动态初始化数组,在C++11及以上版本中,还提供了使用std::array... 目录1、初始化一维数组1.1、使用列表初始化（推荐方式）1.2、初始化部分列表1.3、使用std::
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2