高斯牛顿迭代法的原理及实现（经典例子，附C和C++代码，含运行结果）

可帮助理解的一些文章：

梯度，雅克比矩阵和海森矩阵

海森矩阵和牛顿法

如何通俗易懂地讲解牛顿迭代法？

最优化方法：牛顿迭代法和拟牛顿迭代法

经典举例：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>#define MAXN 7
typedef struct {                            //n*m矩阵int m,n;double data[MAXN][MAXN];
}mat;mat CreateMatrix(int p,int q);
double arry2matrix(mat *a,mat *b,mat *r,double *x,double *y,double *beta,int m,int n);
mat trans(mat a);
mat mul(mat a,mat b);
double *SolveMatrix(mat c,mat d);
void arryans(double *beta,mat c,mat d);int main(){double x[]={0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740};double y[]={0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317};mat r = {0,0,0};mat a = {0,0,0};mat b = {0,0,0};mat c = {0,0,0};mat d = {0,0,0};double beta[2] = {0.9,0.2};double rss = 0;int T;for(T=0;T<10;++T){rss = arry2matrix(&a,&b,&r,x,y,beta,7,2);printf("%d %f %f %f\n",T,beta[0],beta[1],rss);c = mul(b,a);d = mul(b,r);arryans(beta,c,d);}system("pause");
}// Input x[],y[],m,n
// Output a,b,r[]
double arry2matrix(mat *a,mat *b,mat *r,double *x,double *y,double *beta,int m,int n){int i;double rss = 0;(*a).m = m;        // 雅克比矩阵a(*a).n = n;(*b).m = n;       // a的转置矩阵b(*b).n = m;(*r).m = m;(*r).n = 1;for(i=0;i<m;++i){(*a).data[i][0] = -x[i]/(beta[1]+x[i]);(*a).data[i][1] = beta[0]*x[i]/(beta[1]+x[i])/(beta[1]+x[i]);(*b).data[0][i] = (*a).data[i][0];(*b).data[1][i] = (*a).data[i][1];(*r).data[i][0] = y[i]-beta[0]*x[i]/(beta[1]+x[i]);rss += (*r).data[i][0]*(*r).data[i][0];  }return rss;
}/********************************************************************
*Fuction:
*Input：  a,b
*Output:  c
********************************************************************/
mat mul(mat a,mat b){int i,j,k;mat c = {0,0,0};c.m = a.m;c.n = b.n;for(i=0;i<c.m;++i){for(j=0;j<c.n;++j){c.data[i][j] = 0;for(k=0;k<a.n;k++){c.data[i][j] += a.data[i][k]*b.data[k][j];}}}return c;
}/********************************************************************
*Fuction:
*Input：  c,d,beta
*Output:  new_beta
********************************************************************/
void arryans(double *beta,mat c,mat d){                 //解一个n阶的线性方程组int i, j, k, mi;double mx, tmp;double beta_tmp[2];beta_tmp[0] = beta[0];beta_tmp[1] = beta[1];for (i = 0; i < c.m-1; i++) {mi = i;mx = fabs(c.data[i][i]);for (j = i + 1; j < c.m; j++)           //找主元素if (fabs(c.data[j][i]) > mx) {mi = j;                         //记录矩阵下三角一列中绝对值最大值的行号 mx = fabs(c.data[j][i]);     //记录矩阵下三角一列中绝对值最大值 }if (i < mi) {                            //交换两行tmp = d.data[i][0];d.data[i][0] = d.data[mi][0];d.data[mi][0] = tmp;for (j = i; j < c.m; j++) {tmp = c.data[i][j];c.data[i][j] = c.data[mi][j];c.data[mi][j] = tmp;}}//---高斯消元---//for (j = i + 1; j < c.m; j++) {             //第i行，第i+1、i+2、…、n-1、n列tmp = -c.data[j][i] / c.data[i][i];        //第j行除以第i行d.data[j][0] += d.data[i][0] * tmp;for (k = i; k < c.m; k++)              //将第i行的tmp倍加到第j行c.data[j][k] += c.data[i][k] * tmp;       }}beta[c.m - 1] = d.data[c.m - 1][0] / c.data[c.m - 1][c.m - 1];       //求解方程for (i = c.m - 2; i >= 0; i--) {beta[i] = d.data[i][0];for (j = i + 1; j < c.m; j++)beta[i] -= c.data[i][j] * beta[j];beta[i] /= c.data[i][i];}beta[0] = beta_tmp[0]-beta[0];         beta[1] = beta_tmp[1]-beta[1];
}

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
using namespace std;int n;
double x[10],y[10],beta[10];struct matrix{double a[10][10],tmp[10][10];int n,m;void clear(){                            // 清零memset(a,0,sizeof(a));memset(tmp,0,sizeof(tmp));} void cpy(matrix &b){                 // 复制矩阵n=b.n;m=b.m;for (int i=1;i<=n;i++)    for (int j=1;j<=m;j++)a[i][j]=b.a[i][j];}void mul(matrix &b){for (int i=0;i<=n;i++) for (int j=0;j<=b.m;j++) tmp[i][j]=0;for (int i=0;i<=n;i++)for (int k=0;k<=m;k++)if (a[i][k]){for (int j=0;j<=b.m;j++)tmp[i][j]+=a[i][k]*b.a[k][j];a[i][k]=0;}for (int i=0;i<=n;i++)for (int j=0;j<=b.m;j++)a[i][j]=tmp[i][j];m=b.m;}void getinv(){for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) a[i][n+j]=0;for (int i=1;i<=n;i++) a[i][n+i]=1;for (int i=1;i<=n;i++){int po;double maxi=0;for (int j=i;j<=n;j++){if (fabs(a[j][i])>maxi){maxi=fabs(a[j][i]);po=j;}}for (int j=1;j<=2*n;j++){double t=a[i][j];a[i][j]=a[po][j];a[po][j]=t;}if (fabs(maxi)==0) continue;for (int j=i+1;j<=n;j++){double tim=a[j][i]/a[i][i];for (int k=i;k<=2*n;k++) a[j][k]-=a[i][k]*tim;}}for (int i=n;i>=1;i--){for (int j=i+1;j<=n;j++){for (int k=n+1;k<=2*n;k++)a[i][k]-=a[i][j]*a[j][k];a[i][j]=0;          }for (int j=n+1;j<=2*n;j++)a[i][j]/=a[i][i];a[i][i]=1;  }for (int i=1;i<=n;i++)for (int j=1;j<=n;j++)a[i][j]=a[i][j+n];for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) a[i][j+n]=0;}void trav(){for (int i=1;i<=m;i++)for (int j=1;j<=n;j++)tmp[i][j]=a[j][i],a[j][i]=0;for (int i=1;i<=m;i++)for (int j=1;j<=n;j++)a[i][j]=tmp[i][j];swap(n,m);}
}a,b,c,d;int main(){      n = 7;double x[]={0,0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740};double y[]={0,0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317};double r = 0;double rss = 0;beta[1]=0.9;beta[2]=0.2;for (int i=1;i<=n;i++){r = y[i]-beta[1]*x[i]/(beta[2]+x[i]);rss += r*r;}printf("0 %f %f %f\n",beta[1],beta[2],rss);for (int T=1;T<=9;T++){a.clear();b.clear();c.clear();d.clear();    a.n=n;a.m=2;for (int i=1;i<=n;i++){a.a[i][1]=-x[i]/(beta[2]+x[i]);a.a[i][2]=beta[1]*x[i]/(beta[2]+x[i])/(beta[2]+x[i]);}b.cpy(a);c.cpy(a);// 残差矩阵d.n=n;d.m=1;for (int i=1;i<=n;i++){d.a[i][1]=y[i]-beta[1]*x[i]/(beta[2]+x[i]);}a.trav();           // a转置a.mul(b);         // a.getinv();c.trav();         // a.mul(c);a.mul(d);           //beta[1]=beta[1]-a.a[1][1];beta[2]=beta[2]-a.a[2][1];double rss=0;for (int i=1;i<=n;i++)rss+=(y[i]-beta[1]*x[i]/(beta[2]+x[i]))*(y[i]-beta[1]*x[i]/(beta[2]+x[i]));printf("%d %f %f %f\n",T,beta[1],beta[2],rss);}system("pause");
}

运行结果与C语言运行结果一致。

可参考博客：

高斯-牛顿迭代

详解高斯牛顿迭代法原理和代码

高斯牛顿迭代法的原理及实现（经典例子，附C和C++代码，含运行结果）相关推荐

梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现
---------------------梯度下降法------------------- 梯度的一般解释: f(x)在x0的梯度:就是f(x)变化最快的方向.梯度下降法是一个最优化算法,通常也称为最 ...
高斯牛顿迭代法matlab代码,优化算法--牛顿迭代法
简书同步更新牛顿法给出了任意方程求根的数值解法,而最优化问题一般会转换为求函数之间在"赋范线性空间"的距离最小点,所以,利用牛顿法去求解任意目标函数的极值点是个不错的思路. 方程 ...
经典重读 | 用高斯牛顿的方法来进行IKF的更新步骤
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达作者丨泡泡编辑组来源丨泡泡机器人SLAM 点击进入->3D视觉工坊学习交流群标题:The ...
艾特肯法方程解matlab程序,牛顿迭代法matlab代码
牛顿法迭代公式: x(k1) xk [2 f (x(k) )]1f (x(k) ) Matlab 代码: function [x1,k] =newton(x1,eps) hs=inline('(x ...
Matlab利用牛顿迭代法求解非线性方程组
我们得首先了解牛顿迭代法的原理是什么: 在这里,我以二阶非线性方程组为例: f1(x,y)=0 f2(x,y)=0,求解x,y 原理假设方程组的一组近似解为(x0,y0),将方程f1(x,y)=0与 ...
100个python算法超详细讲解：牛顿迭代法求方程根
1．问题描述编写用牛顿迭代法求方程根的函数.方程为ax 3 +bx 2 +cx+d=0,系数a. b.c.d由主函数输入,求x在1附近的一个实根.求出根后,由主函数输出. 2．问题分析牛顿迭代法是 ...
用c语言编制牛顿法程序,求解试用newton法求函数,YTU 2405: C语言习题牛顿迭代法求根...
2405: C语言习题牛顿迭代法求根时间限制: 1 Sec 内存限制: 128 MB 提交: 562 解决: 317 题目描述用牛顿迭代法求根.方程为ax3+bx2+cx+d=0.系数a,b ...
牛顿迭代法python_python 牛顿迭代法
使用牛顿迭代法求方程在x附近的一个实根. 赋值X,即迭代初值:用初值x代入方程中计算此时的f(x)=(a * x * x * x + b * x * x + c * x + d)和f'(x)=(3 ...
【清橙A1094】【牛顿迭代法】牛顿迭代法求方程的根
问题描述给定三次函数f(x)=ax3+bx2+cx+d的4个系数a,b,c,d,以及一个数z,请用牛顿迭代法求出函数f(x)=0在z附近的根,并给出迭代所需要次数. 牛顿迭代法的原理如下(参考下图) ...
java 牛顿迭代算术平方根,牛顿迭代法求n方根
一.简单推导二.使用借助上述公式,理论上可以求任意次方根,假设要求a(假设非负)的n次方根,则有xn=a,令f(x)=xn-a,则只需求f(x)=0时x的值即可.由上述简单推导知,当f(x)=0时 ...

高斯牛顿迭代法的原理及实现（经典例子，附C和C++代码，含运行结果）

高斯牛顿迭代法的原理及实现（经典例子，附C和C++代码，含运行结果）相关推荐

最新文章

热门文章