matlab練習程式(高斯牛頓法最優化)

Dsp Tian發表於2019-01-03

計算步驟如下:

圖片來自《視覺slam十四講》6.2.2節。

下面使用書中的練習y=exp(a*x^2+b*x+c)+w這個模型驗證一下,其中w為噪聲,a、b、c為待解算係數。

程式碼如下:

clear all;
close all;
clc;

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

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);      %步驟1
for i=1:1000
    
    f = exp(pre(1)*x.^2+pre(2)*x+pre(3));
    g = y-f;                    %步驟2中的誤差 
    
    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];             %步驟2中的雅克比矩陣
    
    delta = inv(J'*J)*J'* g;    %步驟3,inv(J'*J)*J'為H的逆
    
    pcur = pre+delta;           %步驟4
    if norm(delta) <1e-16
        break;
    end
    pre = pcur;
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'

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

相關文章