LM算法求解最小二乘问题

Module VarLMImplicit none!//Nparas为参数个数，Ndata为数据个数,Niters为最大迭代次数Integer     ,Parameter :: Nparas = 2,Ndata = 9,Niters = 50,Fileid = 11!//临时变量，order=1，继续迭代，否则停止迭代Integer                :: order!//循环变量，it为迭代次数循环变量Integer                :: it,i,ii!//x_1为观测数据自变量，y_1为观测数据因变量，两者均为已知数据Real(kind=8),Parameter :: x_1(Ndata) = [0.25,0.5,1.0,1.5,2.0,3.0,4.0,6.0,8.0]Real(kind=8),Parameter :: y_1(Ndata)  = [19.306,18.182,16.126,14.302,12.685,9.978,7.849,4.857,3.005]   !//a0与b0分别为参数猜测初始值,RealA、RealB为真值Real(kind=8),Parameter :: a0 = 0.,b0 = 0.,RealA = 20.5,RealB = 0.24!//误差限度Real(kind=4)           :: eps = 1e-8!//y_est为估计值，d为估计值和实际值y_1之间的误差，rd为数组d的变形Real(kind=8)           :: y_est(Ndata) = 0.,d(Ndata) = 0.,rd(Ndata,1) = 0.!//J为雅可比矩阵，JT为J的转置矩阵，H为海森矩阵Real(kind=8)           :: J(Ndata,Nparas) = 0.,JT(Nparas,Ndata) = 0.,H(Nparas,Nparas) = 0.!//Inv_H_1m为H_1m的逆矩阵Real(kind=8)           :: H_Lm(Nparas,Nparas)=0.,Inv_H_Lm(Nparas,Nparas)=0.!//eye为单位矩阵，delta为步长矩阵Real(kind=8)           :: eye(Nparas,Nparas) = 0.0,delta(Nparas,1) = 0.0!//a_est，b_est分别为反演参数Real(kind=8)           :: a_est = 0.,b_est = 0.,e = 0.!//计算新对应的值Real(kind=8)           :: y_est_Lm(Ndata) = 0.,d_Lm(Ndata) = 0.,e_Lm = 0.!//临时变量，lamda为LM算法的阻尼系数初始值Real(kind=8)           :: a_Lm,b_Lm,lamda = 0.01,v = 10.d0,tempEnd Module VarLMProgram LMInclude 'link_fnl_shared.h'Use VarLMUse LINRG_INTImplicit None!//生成单位矩阵Forall( i = 1:Nparas, ii = 1:Nparas, i==ii ) eye(i,ii) = 1.d0!//第一步：变量赋值order = 1a_est = a0b_est = b0!//第二步：迭代Write(*,"('------------------------------------------------')") loop1: Do it = 1,NitersIf ( order==1 ) Then  !//根据当前估计值，计算雅可比矩阵y_est = a_est*Exp( -b_est*x_1 ) !//根据当前a_est，b_est及x_1，得到函数值y_estd = y_1 - y_est   !//计算已知值y_1与y_est的误差loop2: Do i = 1,Ndata  !//计算雅可比矩阵。dy/da = Exp(-b*x)，dy/db = -a*x*Exp(-b*x)J(i,1) = Exp( -b_est*x_1(i) ) J(i,2) = -a_est*x_1(i)*Exp( -b_est*x_1(i) )End Do loop2    JT=Transpose(J)   H = Matmul( JT,J )  !//计算海森矩阵End IfIf ( it==1 ) e = Dot_Product(d,d)  !//若是第一次迭代，计算误差epsilonH_Lm = H + lamda*eye   !//根据阻尼系数lamda混合得到H矩阵!// 縧amda*I縇evenberg靠縧amda*diag(A'A)縈arquardt靠!//计算步长delta，并根据步长计算新的参数估计值Call LINRG( H_Lm,Inv_H_Lm )   !//使用imsl函数库，计算H_Lm的逆矩阵Inv_H_Lmrd = Reshape( d,[Ndata,1] )  !//为了满足内部函数Matmul的计算法则，对d的数组形状进行改变delta = Matmul( Inv_H_Lm,matmul( JT,rd ) )  !//delta为增量a_Lm = a_est + delta(1,1)b_Lm = b_est + delta(2,1)!//如果||delta||<1e-8，终止迭代If ( Dot_Product(delta(:,1),delta(:,1))<eps ) Exit!//计算新的可能估计值对应的y和计算残差ey_est_Lm = a_Lm*Exp( -b_Lm*x_1 )d_Lm = y_1 - y_est_Lme_Lm = Dot_Product( d_Lm,d_Lm )  !//e_LM等于||y_1 - y_est_LM||!//根据误差，决定如何更新参数和阻尼系数!//迭代成功时将lamda减小，否则增大lamdaIf ( e_Lm<e ) Thenlamda = lamda/va_est = a_Lmb_est = b_Lme = e_Lmorder = 1Elseorder = 0lamda = lamda*vEnd IfWrite(*,"('a_est=',g0,2X,'b_est=',g0)") a_est,b_estEnd Do loop1Write(*,"('------------------------------------------------')") Write(*,"('停止迭代，总共迭代',g0,'次')") it - 1!//输出正反演数据以及百分比误差输出到文件,总共100个数据点,计算区域为[0-50]Open(Fileid,file='原始数据与反演数据.dat',status='unknown')Do i = 0,100temp = i/2.d0Write(Fileid,'(f7.3,3f15.8)') temp,RealA*Exp( -RealB*temp ),a_est*Exp( -b_est*temp ),&Abs( a_est*Exp( -b_est*temp )-RealA*Exp( -RealB*temp ) )/RealA/Exp( -RealB*ii )*100End DoClose(Fileid)!//输出反演参数a_est，b_est以及原始数据和拟合数据Write(*,"(/,'反演参数为：')")Write(*,"('a_est=',g0)") a_estWrite(*,"('b_est=',g0)") b_estWrite(*,"(/,'原始数据为：')") Write(*,'(3f9.4/)') y_1Write(*,"('拟合数据为：')") Write(*,'(3f9.4/)') a_est*Exp( -b_est*x_1 )Write(*,"('------------------------------------------------')") End Program LM

LM算法求解最小二乘问题相关推荐

非线性优化：徒手实现LM算法
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达本文由知乎作者David LEE授权转载,不得擅自二次转载. 原文链接:https://zhuanla ...
lm opencv 算法_Levenberg–Marquardt算法学习（和matlab的LM算法对比）
回顾高斯牛顿算法,引入LM算法惩罚因子的计算(迭代步子的计算) 完整的算法流程及代码样例 1. 回顾高斯牛顿,引入LM算法根据之前的博文:Gauss-Newton算法学习假设我们研究如 ...
相机校正、张氏标定法、极大似然估计/极大似然参数估计、牛顿法、高斯牛顿法、LM算法、sin/cos/tan/cot
日萌社人工智能AI:Keras PyTorch MXNet TensorFlow PaddlePaddle 深度学习实战(不定时更新) CNN:RCNN.SPPNet.Fast RCNN.Faste ...
梯度下降、牛顿法、高斯牛顿L-M算法比较
本文梳理一下常见的几种优化算法:梯度下降法,牛顿法.高斯-牛顿法和L-M算法,优化算法根据阶次可以分为一阶方法(梯度下降法),和二阶方法(牛顿法等等),不同优化算法的收敛速度,鲁棒性都有所不同.一般来 ...
Levenberg-Marquardt算法求解单应性矩阵
A. Levenberg-Marquardt算法待估计的模型参数 x=[x1,x2,⋯,xn]T\mathbf{x}=[x_1, x_2, \cdots,x_n]^Tx=[x1,x2,⋯,xn ...
都是被逼的... ，LM算法的具体实现python和C++代码
L-M方法全称Levenberg-Marquardt方法,是一种非线性最小二乘优化算法,它通过同时利用高斯牛顿方法和梯度下降方法来解决非线性最小二乘问题.其核心思想是在每次迭代中,根据当前参数估计计算 ...
奇异值分解（SVD）方法求解最小二乘问题
奇异值分解(SVD)方法求解最小二乘问题 1 奇异值分解(SVD)原理 1.1 回顾特征值和特征向量 1.2 SVD的定义 1.3 求出SVD分解后的U,Σ,V矩阵 1.4 SVD的一些性质 2 线性 ...
【lssvm预测】基于天鹰算法优化最小二乘支持向量机lssvm实现数据回归预测附matlab代码
1 简介短时交通流预测是实现智能交通控制与管理,交通流状态辨识和实时交通流诱导的前提及关键,也是智能化交通管理的客观需要.到目前为止,它的研究结果都不尽如人意.现有的以精确数学模型为基础的传统预测方 ...
数学建模——差分算法(求解偏微分方程)
差分算法(求解偏微分方程) 差分算法是数学建模比赛中的一种十分常见的代码,在2018A题和2020A中均用到一维热传导模型,模型的求解用的就是差分算法,具体如何解可以自己去查看相关论文. 定义差分方 ...
Levenberg-Marquardt(LM算法)的理解
Levenberg-Marquardt LM算法的理解 1. convex optimization 1.1 convex set 1.2 convex function 1.3 optimizat ...

LM算法求解最小二乘问题

LM算法求解最小二乘问题相关推荐

最新文章

热门文章