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');
结果如下：