matlab練習程式（Levenberg-Marquardt法最優化）

Dsp Tian發表於2019-01-04
Matlab優化

上一篇部落格中介紹的高斯牛頓演算法可能會有J'*J為奇異矩陣的情況，這時高斯牛頓法穩定性較差，可能導致演算法不收斂。比如當係數都為7或更大的時候，演算法無法給出正確的結果。

Levenberg-Marquardt法一定程度上修正了這個問題。

計算迭代係數deltaX公式如下：

當lambda很小的時候，H佔主要地位，公式變為高斯牛頓法，當lambda很大的時候，H可以忽略，公式變為最速下降法。該方法提供了更穩定的deltaX。

演算法步驟如下：

1.給定初始係數，以及初始優化半徑u。

2.計算使用當前係數的模型得到的結果與測量結果差值e。

3.使用迭代公式更新帶解算係數。

4.計算更新後係數的模型得到的結果與測量結果差值ecur。

5.如果ecur>e，則u=2*u；否則u=u/2，並且更新模型係數x(k+1)=x(k)+deltaX。

6.判斷演算法是否收斂，不收斂返回2，否則結束。

程式碼如下：

clear all;
close all;
clc;
warning off all;

a=7;b=7;c=7;              %待求解的係數

x=(0:0.01:1)';
w=rand(length(x),1)*2-1;   %生成噪聲
y=exp(a*x.^2+b*x+c)+w;     %帶噪聲的模型 
plot(x,y,'.')

pre=rand(3,1);             
update=1;
u=0.1;
for i=1:100    
    if update==1
        f = exp(pre(1)*x.^2+pre(2)*x+pre(3));
        g = y-f;                                        %計算誤差 

        p1 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x.^2;    %對a求偏導
        p2 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x;       %對b求偏導
        p3 = exp(pre(1)*x.^2+pre(2)*x+pre(3));          %對c求偏導
        J = [p1 p2 p3];                                 %計算雅克比矩陣
        H=J'*J;
        if i==1
            e=dot(g,g);
        end          
    end
    
    delta = inv(H+u*eye(length(H)))*J'* g;      
    pcur = pre+delta;                           %迭代
    fcur = exp(pcur(1)*x.^2+pcur(2)*x+pcur(3)); 
    ecur = dot(y-fcur,y-fcur);
    
    if ecur<e                                   %比較兩次差值，新模型好則使用
        if norm(pre-pcur)<1e-10
           break; 
        end
        u=u/2;  
        pre=pcur;
        e=ecur;
        update=1;    
    else
        u=u*2;        
        update=0;
    end 
end

hold on;
plot(x,exp(a*x.^2+b*x+c),'r');
plot(x,exp(pre(1)*x.^2+pre(2)*x+pre(3)),'g');

%比較一下
[a b c]
pre'

迭代結果，其中散點為帶噪聲資料，紅線為原始模型，綠線為解算模型

參考：

《視覺slam十四講》

http://www.docin.com/p-63281100.html

相關文章