LM擬合算法

一、  Levenberg-Marquardt算法

（1）y=a*e.^(-b*x)形式擬合

clear all
% 計算函數f的雅克比矩陣，是解析式
syms a b y x real;
f=a*exp(-b*x);
Jsym=jacobian(f,[a b]);
% 擬合用數據。參見《數學試驗》，p190，例2
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
load fla
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
data_1=1:26;
obs_1=fla3_1;%y軸的值

% 2. LM算法
% 初始猜測s
a0=1e+09; b0=0;
y_init = a0*exp(-b0*data_1);
% 數據個數
Ndata=length(obs_1);
% 參數維數
Nparams=2;
% 迭代最大次數
n_iters=50;
% LM算法的阻尼系數初值
lamda=0.01;
% step1: 變量賦值
updateJ=1;
a_est=a0;
b_est=b0;
%% 
% step2: 迭代
for it=1:n_iters %迭代次數
   if updateJ==1
       % 根據當前估計值，計算雅克比矩陣
       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))];%雅克比矩陣的生成
                                                                                 % 對系數a和b求導
       end
       % 根據當前參數，得到函數值
       y_est = a_est*exp(-b_est*data_1);
       % 計算誤差
       d=obs_1-y_est;
       % 計算（擬）海塞矩陣
       H=J'*J;
       % 若是第一次迭代，計算誤差
       if it==1
           e=dot(d,d);%向量點乘
       end
   end
   % 根據阻尼系數lamda混合得到H矩陣
   H_lm=H+(lamda*eye(Nparams,Nparams)); % J'*J+lamda*I
   % 計算步長dp，並根據步長計算新的可能的參數估計值
   dp=inv(H_lm)*(J'*d(:)); %deltp=[J'*J+lamda*I]^(-1)*(J'*εk)矩陣求逆
   g = J'*d(:);
   a_lm=a_est+dp(1);
   b_lm=b_est+dp(2);
   % 計算新的可能估計值對應的y和計算殘差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);
   % 根據誤差，決定如何更新參數和阻尼系數
   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
%% 
%顯示優化的結果
a_est
b_est
%從圖形上觀察擬合程度
plot(data_1,obs_1,'*')
hold on
x=1:54;
y=a_est*exp(-b_est*x);
plot(x,y,'o-')

（2）y=a*x.^b+c形式

注意公式變化后，求導的變化 

clear all
% 計算函數f的雅克比矩陣，是解析式
syms a b c x ;
f=a*x.^b+c;
Jsym=jacobian(f,[a b c]);
% 擬合用數據。參見《數學試驗》，p190，例2
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
load fla
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
data_1=1:length(fla3_1);
obs_1=fla3_1;%y軸的值

% 2. LM算法
% 初始猜測s
a0=1e+09; b0=0;c0=1.4e+09;
y_init = a0*data_1.^b0+c0;
% 數據個數
Ndata=length(obs_1);
% 參數維數
Nparams=3;
% 迭代最大次數
n_iters=50;
% LM算法的阻尼系數初值
lamda=0.01;
% step1: 變量賦值
updateJ=1;
a_est=a0;
b_est=b0;
c_est=c0;
%% 
% step2: 迭代
for it=1:n_iters %迭代次數
   if updateJ==1
       % 根據當前估計值，計算雅克比矩陣
       J=zeros(Ndata,Nparams);%數據的維度，系數的維度
       for i=1:length(data_1)
           J(i,:)=[data_1(i).^b_est a_est*log(data_1(i))*data_1(i).^b_est 0];%雅克比矩陣的生成
                                                                                 % 對系數a和b求導
       end
       % 根據當前參數，得到函數值
       y_est =a_est*data_1.^b_est+c_est ;
       % 計算誤差
       d=obs_1-y_est;
       % 計算（擬）海塞矩陣
       H=J'*J;
       % 若是第一次迭代，計算誤差
       if it==1
           e=dot(d,d);%向量點乘
       end
   end
   % 根據阻尼系數lamda混合得到H矩陣
   H_lm=H+(lamda*eye(Nparams,Nparams)); % J'*J+lamda*I
   % 計算步長dp，並根據步長計算新的可能的參數估計值
   dp=inv(H_lm)*(J'*d(:)); %deltp=[J'*J+lamda*I]^(-1)*(J'*εk)矩陣求逆
   g = J'*d(:);
   a_lm=a_est+dp(1);
   b_lm=b_est+dp(2);
   c_lm=c_est+dp(3);
   % 計算新的可能估計值對應的y和計算殘差e
   y_est_lm = a_lm*data_1.^b_lm+c_lm;
   d_lm=obs_1-y_est_lm;
   e_lm=dot(d_lm,d_lm);
   % 根據誤差，決定如何更新參數和阻尼系數
   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
%% 
%顯示優化的結果
a_est
b_est
c_est
%從圖形上觀察擬合程度
plot(data_1,obs_1,'*')
hold on
x=1:length(data_1);
y=a_est*data_1.^b_est+c_est;
plot(x,y,'o-')

 

擬合不是很好

更改初始參數c0=1.42e+09后，變為如下圖，擬合效果好很多，所以參數調整很重要。

LM函數

function [a_est,b_est,c_est]=LM(data)
    data_1=1:length(data);
    obs_1=data;%y軸的值

    % 2. LM算法
    % 初始猜測s
    a0=1; b0=0;c0=1.42e+09;
    y_init = a0*data_1.^b0+c0;
    % 數據個數
    Ndata=length(obs_1);
    % 參數維數
    Nparams=3;
    % 迭代最大次數
    n_iters=500;
    % LM算法的阻尼系數初值
    lamda=0.01;
    % step1: 變量賦值
    updateJ=1;
    a_est=a0;
    b_est=b0;
    c_est=c0;
    %% 
    % step2: 迭代
    for it=1:n_iters %迭代次數
       if updateJ==1
           % 根據當前估計值，計算雅克比矩陣
           J=zeros(Ndata,Nparams);%數據的維度，系數的維度
           for i=1:length(data_1)
               J(i,:)=[data_1(i).^b_est a_est*log(data_1(i))*data_1(i).^b_est 0];%雅克比矩陣的生成
                                                                                     % 對系數a和b求導
           end
           % 根據當前參數，得到函數值
           y_est =a_est*data_1.^b_est+c_est ;
           % 計算誤差
           d=obs_1-y_est;
           % 計算（擬）海塞矩陣
           H=J'*J;
           % 若是第一次迭代，計算誤差
           if it==1
               e=dot(d,d);%向量點乘
           end
       end
       % 根據阻尼系數lamda混合得到H矩陣
       H_lm=H+(lamda*eye(Nparams,Nparams)); % J'*J+lamda*I
       % 計算步長dp，並根據步長計算新的可能的參數估計值
       dp=inv(H_lm)*(J'*d(:)); %deltp=[J'*J+lamda*I]^(-1)*(J'*εk)矩陣求逆
       g = J'*d(:);
       a_lm=a_est+dp(1);
       b_lm=b_est+dp(2);
       c_lm=c_est+dp(3);
       % 計算新的可能估計值對應的y和計算殘差e
       y_est_lm = a_lm*data_1.^b_lm+c_lm;
       d_lm=obs_1-y_est_lm;
       e_lm=dot(d_lm,d_lm);
       % 根據誤差，決定如何更新參數和阻尼系數
       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

　

 

 二、nlinfit擬合

clear;clc
x=[36 39 45 48 53 61 73 97 108 121 140 152];
y=[10000 7000 6500 5500 4500 3500 2000 300 750 1000 250 100];
a0=[0 0 0];
a=nlinfit(x,y,'gx',a0);
s=30:0.01:160;
plot(x,y,'*',x,gx(a,x),'--r')
a=vpa(a,5)

function q=gx(a,x)
q=a(1)*x.^a(2)+a(3);
end

　

 三、最小二乘法

%最小二乘法
function [a,b]=Least_sqr(data)%%
    x=1:length(data);
    y=data;
    s_xy=sum(x.*y);
    s_xx=sum(x.*x);
    s_sxy=sum(x)*sum(y)/length(x);
    s_sx=sum(x).*sum(x)/length(x);
    ave_x=sum(x)/length(x);
    ave_y=sum(y)/length(y);
    a=(s_xy-s_sxy)/(s_xx-s_sx);
    b=ave_y-a*ave_x;
end

　

 四、對LM算法的介紹

LM算法的實現並不算難，它的關鍵是用模型函數 f 對待估參數向量p在其鄰域內做線性近似，忽略掉二階以上的導數項，從而轉化為線性最小二乘問題，它具有收斂速度快等優點。LM算法屬於一種“信賴域法”——所謂的信賴域法，此處稍微解釋一下：在最優化算法中，都是要求一個函數的極小值，每一步迭代中，都要求目標函數值是下降的，而信賴域法，顧名思義，就是從初始點開始，先假設一個可以信賴的最大位移s，然后在以當前點為中心，以s為半徑的區域內，通過尋找目標函數的一個近似函數（二次的）的最優點，來求解得到真正的位移。在得到了位移之后，再計算目標函數值，如果其使目標函數值的下降滿足了一定條件，那么就說明這個位移是可靠的，則繼續按此規則迭代計算下去；如果其不能使目標函數值的下降滿足一定的條件，則應減小信賴域的范圍，再重新求解。

事實上，你從所有可以找到的資料里看到的LM算法的說明，都可以找到類似於“如果目標函數值增大，則調整某系數再繼續求解；如果目標函數值減小，則調整某系數再繼續求解”的迭代過程，這種過程與上面所說的信賴域法是非常相似的，所以說LM算法是一種信賴域法。

 

LM算法需要對每一個待估參數求偏導，所以，如果你的目標函數f非常復雜，或者待估參數相當地多，那么可能不適合使用LM算法，而可以選擇Powell算法——Powell算法不需要求導。

 

至於這個求導過程是如何實現的，我還不能給出建議，我使用過的方法是拿到函數的方程，然后手工計算出其偏導數方程，進而在函數中直接使用，這樣做是最直接，求導誤差也最小的方式。不過，在你不知道函數的形式之前，你當然就不能這樣做了——例如，你提供給了用戶在界面上輸入數學函數式的機會，然后在程序中解析其輸入的函數，再做后面的處理。在這種情況下，我猜是需要使用數值求導算法的，但我沒有親自試驗過這樣做的效率，因為一些優秀的求導算法——例如Ridders算法——在一次求導數值過程中，需要計算的函數值次數也會達到5次以上。這樣的話，它當然要比手工求出導函數（只需計算一次，就可以得到導數值）效率要差得多了。不過，我個人估計（沒有任何依據的，只是猜的）：依賴於LM算法的高效，就算添加了一個數值求導的“拖油瓶”，整個最優化過程下來，它仍然會優於Powell等方法。

關於偏導數的求取

個人認為：在條件允許、對速度和精度任何以方面都有一定要求的前提下，如果待求解的函數形式是顯式的，應當盡量自己計算目標函數的偏導數方程。原因在於，在使用數值法估計偏導數值時，盡管我們可以控制每一步偏導數值的精度，但是由於求解過程需要進行多次迭代，特別是收斂過程比較慢的求解過程，需要進行很多次的求解，每一次求解的誤差偏差都會在上一步偏差的基礎上不斷累積。盡管在最后依然可以收斂，但是得到的解已經離可以接受的解偏離比較遠了。因此，在求解函數形式比較簡單、偏導數函數比較容易求取時，還是盡量手動計算偏導數，得到的結果誤差相對更小一些。

 

在這篇解釋信賴域算法的文章中，我們已經知道了LM算法的數學模型：

可以證明，此模型可以通過解方程組(Gk+μI)s=−gk確定sk來表征。
即：LM算法要確定一個μ≥0，使得Gk+μI正定，並解線性方程組(Gk+μI)sk=−gk求出sk。
下面來看看LM算法的基本步驟：

·從初始點x0，μ0>0開始迭代

·到第k步時，計算xk和μk

·分解矩陣Gk+μkI，若不正定，令μk=4μk並重復到正定為止

·解線性方程組(Gk+μkI)sk=−gk求出sk並計算rk

·若rk<0.25，令μk+1=4μk；若rk>0.75，令μk+1=μk2；若0.25≤rk≤0.75，令μk+1=μk

·若rk≤0，說明函數值是向着上升而非下降的趨勢變化了（與最優化的目標相反），這說明這一步走錯了，而且錯得“離譜”，此時，不應該走到下一點，而應“原地踏步”，即xk+1=xk，並且和上面rk<0.25的情況一樣對μk進行處理。反之，在rk>0的情況下，都可以走到下一點，即xk+1=xk+sk

·        迭代的終止條件：∥gk∥<ε，其中ε是一個指定的小正數（大家可以想像一下二維平面上的尋優過程（函數圖像類似於拋物線），當接近極小值點時，迭代點的梯度趨於0）

從上面的步驟可見，LM求解過程中需要用到求解線性方程組的算法，一般我們使用高斯約當消元法，因為它非常穩定——雖然它不是最快最好的算法。
同時，上面的算法步驟也包含對矩陣進行分解的子步驟。為什么要先分解矩陣，再解線性方程組？貌似是這樣的（數學不好的人再次淚奔）：不分解矩陣使之正定，就無法確定那個線性方程組是有解的。矩陣分解有很多算法，例如LU分解等，這方面我沒有看。

https://www.cnblogs.com/shhu1993/p/4878992.html 

三，程序實現

推導過程的博客 https://www.cnblogs.com/monoSLAM/p/5249339.html

clear all
% 計算函數f的雅克比矩陣，是解析式
syms a b c x ;
f=a*x.^b+c;%函數形式
Jsym=jacobian(f,[a b c]);
data=[1425204748.0, 1425204757.0, 1425204759.0, 1425204766.0,...
       1425204768.0, 1425206629.0, 1425209138.0, 1425209139.0,...
       1425209159.0, 1425209167.0, 1425219290.0, 1425219298.0,...
       1425219307.0, 1425219307.0, 1425219381.0, 1425219382.0,...
       1425219390.0,1425223403.0, 1425223411.0, 1425223420.0,...
       1425225096.0, 1425479571.0, 1425479580.0, 1427109036.0,...
       1427109038.0, 1427115143.0, 1427449736.0, 1427449755.0,...
       1427449763.0, 1427449774.0, 1427449775.0, 1427717884.0];
data_1=1:length(data);
obs_1=data;%y軸的值

% 2. LM算法
% 初始猜測s
a0=1; b0=0;c0=data(1);  %設置初始值
y_init = a0*data_1.^b0+c0;
% 數據個數
Ndata=length(obs_1);
% 參數維數
Nparams=3;
% 迭代最大次數
n_iters=200;
% LM算法的阻尼系數初值
lamda=0.01;
% step1: 變量賦值
updateJ=1
a_est=a0;
b_est=b0;
c_est=c0;
%% 
% step2: 迭代
for it=1:n_iters %迭代次數
    if updateJ==1
       % 根據當前估計值，計算雅克比矩陣
       J=zeros(Ndata,Nparams);%數據的維度，系數的維度
       for i=1:length(data_1)
           J(i,:)=[data_1(i).^b_est a_est*log(data_1(i))*data_1(i).^b_est 0];%雅克比矩陣的生成,32*3矩陣
                                                                                 % 對系數a、b、c求導
       end
       % 根據當前參數，得到函數值
       y_est =a_est*data_1.^b_est+c_est ;%32個預測值
       % 計算誤差
       d=obs_1-y_est;
       % 計算（擬）海塞矩陣，是一個多元函數的二階偏導數構成的方陣，描述了函數的局部曲率，對稱矩陣
       H=J'*J;
       % 若是第一次迭代，計算誤差
       if it==1
           e=dot(d,d);%向量點乘
       end
    end
    
   % 根據阻尼系數lamda混合得到H矩陣
   H_lm=H+(lamda*eye(Nparams,Nparams)); % J'*J+lamda*I
   % 計算步長dp，並根據步長計算新的可能的參數估計值
   dp=inv(H_lm)*(J'*d'); %deltp=[J'*J+lamda*I]^(-1)*(J'*εk)矩陣求逆
   g = J'*d(:);
   
   a_lm=a_est+dp(1);%新的系數值
   b_lm=b_est+dp(2);
   c_lm=c_est+dp(3);
   
   % 計算新的可能估計值對應的y和計算殘差e
   y_est_lm = a_lm*data_1.^b_lm+c_lm;
   d_lm=obs_1-y_est_lm;
   e_lm=dot(d_lm,d_lm);
   
   % 根據誤差，決定如何更新參數和阻尼系數
   if e_lm<e  %誤差變小
       lamda=lamda/10;
       a_est=a_lm;
       b_est=b_lm;
       c_est=c_lm;
       e=e_lm;
       disp(e);
       updateJ=1;
   else
       updateJ=0;
       lamda=lamda*10;
   end
end
%% 
%顯示優化的結果
a_est
b_est
c_est
%從圖形上觀察擬合程度
plot(data_1,obs_1,'*')
hold on
x=1:length(data);
y=a_est*data_1.^b_est+c_est;
plot(x,y,'o-')