LM算法的C++实现

本文提供了一个使用Levenberg-Marquardt (LM) 算法进行数据拟合的实例,通过C++与OpenCV实现,详细展示了算法的迭代过程及关键步骤。此示例有助于理解LM算法的工作原理。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

这是一个数据拟合的例子,并没有采用面向对象的设计方法是使能更好的理解LM算法的流程,简约而不简单。算法详细过程不多介绍。程序中用到opencv库中的矩阵类Mat。

例:

#pragma  once
#include <stdio.h>
#include "opencv2\core\core.hpp"

#pragma comment(lib,"opencv_core248d.lib")

const int  MAXTIME = 50;

using namespace cv;
FileStorage fs;

Mat jacobin(const Mat& pk/*[a,b]*/, const Mat& x); //f = a*exp(-b*x)
Mat yEstimate(const Mat& p, const Mat& x);
inline void outData(FileStorage& fs, Mat & m, char* filename)
{
	fs.open(filename, FileStorage::FORMAT_XML | FileStorage::WRITE);
	char *temp = new char[10]; 
	strcpy_s(temp, 10,filename); 
	*strchr(temp, '.') = '\0';
	fs << temp << m; 
	fs.release();
	delete[] temp;
}

void LM(double* p0, int pN, double* x, int xN, double* y, double lamda, double step, double ep = 0.0001)
{
	int iters = 0;
	int updateJ = 1;
	double ek = 0.0, ekk = 0.0;//估计误差
	Mat_<double> xM(xN, 1, x), yM(xN, 1, y), pM(pN, 1, p0), JM, yEM, yEMM, dM, gM, dMM, dpM;//至少需要JM,gM,dpM,pM
	for (; iters < MAXTIME; iters++)
	{
		if (updateJ == 1)
		{
			JM = jacobin(pM, xM);
			//outData(fs, JM, "J.xml");
			yEM = yEstimate(pM, xM);
			dM = yM - yEM;
			gM = JM.t()*dM;
			if (iters == 0)
				ek = dM.dot(dM);
			//outData(fs, dM, "d.xml");
		}
		Mat_<double> NM = JM.t()*JM + lamda*(Mat::eye(pN, pN, CV_64F));
		if (solve(NM, gM, dpM))
		{
			Mat_<double> pMM = pM + dpM;
			yEMM = yEstimate(pMM, xM);
			dMM = yM - yEMM;
			ekk = dMM.dot(dMM);
			//outData(fs, dMM, "dlm.xml");
			//outData(fs, dpM, "dp.xml");
			if (ekk < ek)//成功则更新向量与估计误差
			{
				printf("the %d iterator result\n", iters);
				if (dpM.dot(dpM) < ep)
				{
					outData(fs, pM, "p.xml");
					return;
				}
				else
				{
					pM = pMM;
					ek = ekk;
					lamda = lamda / step;
					updateJ = 1;
					continue;
				}
			}
			else
			{
				lamda = lamda*step;
				updateJ = 0;
			}
		}
		else
		{
			outData(fs, JM, "badJ.xml");
			//return;
		}
	}
}

Mat jacobin(const Mat& pk/*[a,b]*/, const Mat& x)
{
	Mat_<double> J(x.rows, pk.rows), da, db;
	exp(-pk.at<double>(1)*x, da);
	exp(-pk.at<double>(1)*x, db);
	db = -pk.at<double>(0)*x.mul(db);
	//outData(fs, da, "da.xml");
	//outData(fs, db, "db.xml");
	da.copyTo(J(Rect(0, 0, 1, J.rows)));
	db.copyTo(J(Rect(1, 0, 1, J.rows)));
	return J;
}
Mat yEstimate(const Mat& p, const Mat& x)
{
	Mat_<double> Y(x.rows, x.cols);
	exp(-p.at<double>(1)*x, Y);
	Y = p.at<double>(0)*Y;
	return Y;
}


#include "LMM.h"

int main()
{
	double data[] = { 0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8 };
	double obs[] = { 19.21, 18.15, 15.36, 14.10, 12.89, 9.32, 7.45, 5.24, 3.01 };
	double p0[] = { 1, 1 };
	LM(p0, 2, data, 9, obs, 0.01, 10);
}



LM算法,全称为Levenberg-Marquard算法,它可用于解决非线性最小二乘问题,多用于曲线拟合等场合。 LM算法实现并不算难,它的关键是用模型函数 f 对待估参数向量 p 在其邻域内做线性近似,忽略掉二阶以上的导数项,从而转化为线性最小二乘问题,它具有收敛速度快等优点。LM算法属于一种“信赖域法”——所谓的信赖域法,此处稍微解释一下:在最优化算法中,都是要求一个函数的极小值,每一步迭代中,都要求目标函数值是下降的,而信赖域法,顾名思义,就是从初始点开始,先假设一个可以信赖的最大位移 s ,然后在以当前点为中心,以 s 为半径的区域内,通过寻找目标函数的一个近似函数(二次的)的最优点,来求解得到真正的位移。在得到了位移之后,再计算目标函数值,如果其使目标函数值的下降满足了一定条件,那么就说明这个位移是可靠的,则继续按此规则迭代计算下去;如果其不能使目标函数值的下降满足一定的条件,则应减小信赖域的范围,再重新求解。 事实上,你从所有可以找到的资料里看到的LM算法的说明,都可以找到类似于“如果目标函数值增大,则调整某系数再继续求解;如果目标函数值减小,则调整某系数再继续求解”的迭代过程,这种过程与上面所说的信赖域法是非常相似的,所以说LM算法是一种信赖域法。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值