Levenberg-Marquardt快速入门教程
本文附的源程序是MATLAB代碼,總共不到80行,實(shí)現(xiàn)了 求雅克比矩陣的解析解,演示了Levenberg-Marquardt最優(yōu)化迭代過(guò)程,演示了如何求解擬合問(wèn)題。本文用圖文介紹了LM算法。轉(zhuǎn)帖請(qǐng)注明來(lái)自蜜蜂電腦,謝謝!
作者:沈樂(lè)君(www.shenlejun.cn)
什么是最優(yōu)化,可分為幾大類?
答:Levenberg-Marquardt算法是最優(yōu)化算法中的一種。最優(yōu)化是尋找使得函數(shù)值最小的參數(shù)向量。它的應(yīng)用領(lǐng)域非常廣泛,如:經(jīng)濟(jì)學(xué)、管理優(yōu)化、網(wǎng)絡(luò)分析 、最優(yōu)設(shè)計(jì)、機(jī)械或電子設(shè)計(jì)等等。
根據(jù)求導(dǎo)數(shù)的方法,可分為2大類。第一類,若f具有解析函數(shù)形式,知道x后求導(dǎo)數(shù)速度快。第二類,使用數(shù)值差分來(lái)求導(dǎo)數(shù)。
根據(jù) 使用模型不同,分為非約束最優(yōu)化、約束最優(yōu)化、最小二乘最優(yōu)化。
什么是Levenberg-Marquardt算法?
它是使用最廣泛的非線性最小二乘算法,中文為列文伯格-馬夸爾特法。它是利用梯度求最大(小)值的算法,形象的說(shuō),屬于“爬山”法的一種。它同時(shí)具有梯度 法和牛頓法的優(yōu)點(diǎn)。當(dāng)λ很小時(shí),步長(zhǎng)等于牛頓法步長(zhǎng),當(dāng)λ很大時(shí),步長(zhǎng)約等于梯度下降法的步長(zhǎng)。在作者的科研項(xiàng)目中曾經(jīng)使用過(guò)多次。圖1顯示了算法從起 點(diǎn),根據(jù)函數(shù)梯度信息,不斷爬升直到最高點(diǎn)(最大值)的迭代過(guò)程。共進(jìn)行了12步。(備注:圖1中綠色線條為迭代過(guò)程)。?
圖1 LM算法迭代過(guò)程形象描述
圖1中,算法從山腳開(kāi)始不斷迭代。可以看到,它的尋優(yōu)速度是比較快的,在山腰部分直接利用梯度大幅度提升(參見(jiàn)后文例子程序中l(wèi)amda較小時(shí)),快到山頂時(shí)經(jīng)過(guò)幾次嘗試(lamda較大時(shí)),最后達(dá)到頂峰(最大值點(diǎn)),算法終止。
?
如何快速學(xué)習(xí)LM算法?
學(xué) 習(xí)該算法的主要困難是入門難。 要么國(guó)內(nèi)中文教材太艱澀難懂,要么太抽象例子太少。目前,我看到的最好的英文入門教程是K. Madsen等人的《Methods for non-linear least squares problems》本來(lái)想把原文翻譯一下,貼到這里。請(qǐng)讓我偷個(gè)懶吧。能找到這里的讀者,應(yīng)該都是E文好手,我翻譯得不清不楚,反而事倍功半了。
可在 下面的鏈接中找到
http://www2.imm.dtu.dk/pubdb/public/publications.php? year=&pubtype=7&pubsubtype=§ion=1&cmd=full_view&lastndays=&order=author
或者直接下載pdf原文:
http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
例子程序(MATLAB源程序)
本程序不到100行,實(shí)現(xiàn)了 求雅克比矩陣的解析解,Levenberg-Marquardt最優(yōu)化迭代,演示了如何求解擬合問(wèn)題。采用《數(shù)學(xué)試驗(yàn)》(第二版)中p190例2來(lái)演示。在MATLAB中可直接運(yùn)行得到最優(yōu)解。
| % 計(jì)算函數(shù)f的雅克比矩陣,是解析式 syms a b y x real; f=a*exp(-b*x); Jsym=jacobian(f,[a b]) ? ? % 擬合用數(shù)據(jù)。參見(jiàn)《數(shù)學(xué)試驗(yàn)》,p190,例2 data_1=[0.25 0.5 1 1.5 2 3 4 6 8]; obs_1=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; ? % 2. LM算法 % 初始猜測(cè)s a0=10; b0=0.5; y_init = a0*exp(-b0*data_1); % 數(shù)據(jù)個(gè)數(shù) Ndata=length(obs_1); % 參數(shù)維數(shù) Nparams=2; % 迭代最大次數(shù) n_iters=50; % LM算法的阻尼系數(shù)初值 lamda=0.01; ? % step1: 變量賦值 updateJ=1; a_est=a0; b_est=b0; ? % step2: 迭代 for it=1:n_iters ??? if updateJ==1 ??????? % 根據(jù)當(dāng)前估計(jì)值,計(jì)算雅克比矩陣 ??????? ? J=zeros(Ndata,Nparams); ??????? for i=1:length(data_1) ??????????? ? J(i,:)=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))]; ?????? ?end ??????? % 根據(jù)當(dāng)前參數(shù),得到函數(shù)值 ??????? y_est = ? a_est*exp(-b_est*data_1); ??????? % 計(jì)算誤差 ??????? ? d=obs_1-y_est; ??????? % 計(jì)算(擬)海塞矩陣 ??????? H=J'*J; ??????? % 若是第一次迭代,計(jì)算誤差 ??????? if it==1 ??????????? ? e=dot(d,d); ??????? end ??? end ? ??? % 根據(jù)阻尼系數(shù)lamda混合得到H矩陣 ? ??H_lm=H+(lamda*eye(Nparams,Nparams)); ??? % 計(jì)算步長(zhǎng)dp,并根據(jù)步長(zhǎng)計(jì)算新的可能的\參數(shù)估計(jì)值 ??? ? dp=inv(H_lm)*(J'*d(:)); ??? g = J'*d(:); ??? ? a_lm=a_est+dp(1); ??? ? b_lm=b_est+dp(2); ??? % 計(jì)算新的可能估計(jì)值對(duì)應(yīng)的y和計(jì)算殘差e ??? y_est_lm = ? a_lm*exp(-b_lm*data_1); ??? ? d_lm=obs_1-y_est_lm; ??? e_lm=dot(d_lm,d_lm); ??? % 根據(jù)誤差,決定如何更新參數(shù)和阻尼系數(shù) ??? if e_lm<e ??????? ? lamda=lamda/10; ??????? ? a_est=a_lm; ??????? ? b_est=b_lm; ??????? e=e_lm; ??????? disp(e); ??????? updateJ=1; ??? else ??????? updateJ=0; ??????? ? lamda=lamda*10; ??? end end %顯示優(yōu)化的結(jié)果 a_est b_est |
演示程序求解的問(wèn)題是《數(shù)學(xué)試驗(yàn)》(蕭樹(shù)鐵,第二版,高等教育出版社)中p190例2。為了方便讀者,提供該書(shū)籍的數(shù)據(jù)和目標(biāo)函數(shù)照片(2012年4月1日)。本文來(lái)自:shenlejun.cn.
C++源碼
?
轉(zhuǎn)載于:https://www.cnblogs.com/rotten_potato/p/3275316.html
總結(jié)
以上是生活随笔為你收集整理的Levenberg-Marquardt快速入门教程的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇: 【转】IOS动画的实现,其实很简单
- 下一篇: CORE ANIMATION的学习备忘录