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');

结果如下：