L-M算法

L-M算法

编程基础 • 2024-12-21 23:51 • 阅读 14

最近深入学习了L-M（Levenberg-Marquardt）算法，感觉还是挺有意思的。帮我解决了三阶指数拟合的问题。在METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS这本书里有详细的介绍。理论上该方法应该可以解决大多数非线性拟合的问题。下面是算法截图

部分代码如下：

 cal_func(x,t,y,func); // f=y-f(x,t) cal_jacob(x,t,jacob); // j=f'(x,t) cal_hessian(jacob,hessian); // A=j(T) * j cal_grad(jacob,func,grad); // g=j(T) * f bFound = (cal_max(grad) <= e1); //结束标志一般初值为0 除非初始解为最优化解 miu = cal_max2(hessian) * 1e-6; //下降长度的初值 for(n=0; n<500 && !bFound; n++) { cal_hess_miu(hessian,miu,hess_miu); // A'=A+miu * I gauss_jordan(hess_miu,grad,h_lm); // h_lm=h^-1 * grad if( cal_norm(h_lm) <= e2*(e2+cal_norm(x)) ) bFound = true; else { for(i=0; i<N; i++) x_new[i]=x[i]+h_lm[i]; // x_new=x+h_lm objFi=cal_objFi(func); //φ=∑f ^ 2 cal_func(x_new,t,y,fun_new); // f_new=y-f(x_new,t) objFi_new=cal_objFi(fun_new); //φ_new = ∑f_new ^ 2 rou = (objFi-objFi_new)/cal_l(h_lm,grad,miu); if(rou>0) { for(i=0; i<N; i++) x[i]=x_new[i]; cal_func(x,t,y,func); // f=y-f(x,t) cal_jacob(x,t,jacob); // j=f'(x,t) cal_hessian(jacob,hessian); // A=j(T) * j cal_grad(jacob,func,grad); // g=j(T) * f bFound = cal_max(grad) <= e1; miu *= (1.0f/3.0f) > (1-pow(2*rou-1,3)) ? (1.0f/3.0f) : (1-pow(2*rou-1,3)); v = 2; cout<<n<<"\t"<<miu<<"\t"<<cal_max(grad)<<"\t"<<objFi<<endl; } else { miu *= v; v *= 2; cout<<n<<"\t"<<miu<<"\t"<<cal_max(grad)<<"\t"<<objFi<<endl; } } }

以下是函数声明

int gauss_jordan(double c[N][N], double d[N], double e[N]); int cal_jacob(double x[N],int t[M],double jacob[M][N]); //计算雅克比矩阵jacob(77*N) int cal_hessian(double jacob[M][N],double hessian[N][N]); //计算(拟)海赛矩阵hessian(N*N) int cal_func(double x[N],int t[M],double y[M],double func[M]); //计算函数f(77*1) double cal_objFi(double func[M]); //计算目标函数objFi(1*1) int cal_grad(double jacob[M][N],double func[M],double grad[N]); //计算梯度gradient(N*1) int cal_hess_miu(double hessian[N][N],double miu,double hess_miu[N][N]); //计算混合矩阵(N*N) double cal_max(double ch[N]); //计算一维最大值 double cal_max2(double ch[N][N]); //计算二维最大值 double cal_norm(double ch[N]); //计算向量的模 double cal_l(double h_lm[N],double grad[N],double miu); //计算分母L(0)-L(h_lm)

需要根据拟合的函数自行编写。

仅供参考，欢迎指正、讨论。

今天的文章 L-M算法分享到此就结束了，感谢您的阅读。

相关推荐