matlab練習程序（BFGS）

BFGS和DFP都是擬牛頓法，和高斯牛頓法不同的地方是不用直接求黑塞矩陣了，而BFGS又比DFP算法有更好的數值穩定性。

算法步驟如下：

1. 給一個待求參數的初始值x(1)。

2. 給定H(1)矩陣為單位陣，並且計算出待優化函數在x(k)處的梯度g(k)。

3. 令d(k) = -H(k)*g(k)，得到搜索方向。

4. 從x(k)出發，沿着d(k)方向搜索，給定一個步長，找到其搜索范圍內比當前參數x(k)更小函數值對應的x值，記為x(k+1)。

5. 計算待優化函數在x(k+1)處的梯度g(k+1)。

6. 計算此時參數位移p = x(k+1) - x(k)和梯度位移q = g(k+1) - g(k)。

7. 最后根據下面迭代公式更新H矩陣即可，當g小於一定閾值時認為迭代結束。

下面兩個公式是對偶的，形式上就是把p和q對換了一下，通常BFGS性能更優。

DFP迭代公式：

BFGS迭代公式：

這里提供一個NIST非線性擬合例題作為示例。

matlab代碼如下：

clear all;close all;clc;
warning off;

data=[109 1;        %待擬合數據
      149 2;
      149 3;
      191 5;
      213 7;
      224 10];

y = data(:,1);
x = data(:,2);
plot(x,y,'ro');

syms b1 b2;
ff = sum((y - (b1*(1-exp(-b2*x)))).^2); %定義待優化函數
dff1 = diff(ff,b1);                     %兩個參數的偏導
dff2 = diff(ff,b2);

f=inline(char(ff),'b1','b2');           %轉為函數
g1=inline(char(dff1),'b1','b2');
g2=inline(char(dff2),'b1','b2');

b = [1;1];                              %初始參數
rho = 0.5;sigma = 0.5;                  %迭代步長
H = eye(2);                             %用來替代hessian矩陣的H矩陣

re=[];
for i=1:200
    g = [g1(b(1),b(2));g2(b(1),b(2))]; 
    dk = -inv(H)*g;
    
    mk=1;   
    for j=1:20              %局部最優搜索
        new=f(b(1)+rho^j*dk(1),b(2)+rho^j*dk(2));   
        old=f(b(1),b(2));
        
        if new < old + sigma*rho^j*g'*dk
            mk = j;
            break;
        end
    end
    
    norm(g)
    if norm(g)<1e-20 || isnan(new)
        break;
    end 
    b = b + rho^mk*dk;      %向局部最優方向移動
    gg=[g1(b(1),b(2));g2(b(1),b(2))];
    
    q = gg - g;             %q = g(k+1)-g(k)
    p = rho^mk*dk;          %p = x(k+1)-x(k)
    H = H - (H*p*p'*H)/(p'*H*p) + (q*q')/(q'*p);    %矩陣更新
end
b

x1 = min(x):0.01:max(x);
y1 = (b(1)*(1-exp(-b(2)*x1)));
hold on;
plot(x1,y1,'b');

結果如下：