L-M算法

最近深入学习了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)

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

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