最优化建模算法理论之BFGS/DFP拟牛顿方法（数学原理及MATLAB实现）

---

一、实现原理

具体数学实现原理可参考这篇文章：拟牛顿法/Quasi-Newton，DFP算法/Davidon-Fletcher-Powell，及BFGS算法/Broyden-Fletcher-Goldfarb-Shanno

---

二、代码实现

1. 拟牛顿方法

算法伪码如下：

2. BFGS算法

采用 BFGS 公式的拟牛顿算法如下：

function [f, xk, k] = BFGS(x0, fun, grid, eps, kmax)
    %
    % function [f, xk, k] = BFGS(x0, fun, grid, eps, kmax)
    % 求出函数fun从初始点x0处以拟牛顿方向为下降方向，
    % 采用 Armjio 准则计算迭代步长，求出函数的极小点
    % -----------------------------------------------------------
    % 输入:
    % 	x0    初始点(列向量)
    % 	fun 	函数文件名称(字符变量）
    %		grid 	梯度函数文件名称(字符变量)
    %		eps   精度要求
    %     kmax  函数的最大迭代次数
    %
    % 输出:
    %		  f	    函数在极小值 xk 处的目标函数值
    %		  xk		函数采用此方法求得的极小点
    %	  	k		  求极小点算法迭代次数
    % -----------------------------------------------------------
    % by Zhi Qiangfeng 
    % 
  k = 0;
  n = length(x0)
  H0 = eye(n); % 初始选取单位阵作为Hessen矩阵的逆的近似阵
  Hk = H0;
  xk = x0;
  gk = feval(grid, xk);
  while k <= kmax
      if norm(gk) < eps
          break;
      end
      dk = -Hk * gk; % 拟牛顿下降方向
      alpha = Armjio(fun, grid, xk, dk);
      x_ = xk; % x_ 保存上一个点坐标
      xk = x_ + alpha * dk; % 更新 xk
      gk_ = gk; % gk_ 保存上一个点的梯度值
      gk = feval(grid, xk); % 更新 gk
      sk = xk - x_; % 记 xk - x_ 为 sk
      yk = gk - gk_; % 记 gk - gk_ 为 yk
      if sk' * yk > 0
          v = yk' * sk;
          % BFGS公式
          Hk = Hk + (1 + (yk' * Hk * yk) / v) * (sk * sk') / v - (sk * yk' * Hk + Hk * yk * sk') / v;
      end
      k = k + 1;
  end
  f = feval(fun, xk);
end

---

3. DFP算法

注：若采用 DFP 算法，则只需将上述代码中的第 40、41、42 行改为如下 DFP 公式即可，拟牛顿方法主体算法一致。

v = Hk * yk;
% DFP公式
Hk = Hk + (sk * sk') / (sk' * yk) - (v * v') / (yk' * v);

---

三、文件化输出

以 Rosenbrock 函数为例，这是优化领域中一个著名的检验函数，其函数与其梯度函数如下：

编写函数文件 Rosenbrock.m 如下：

function f = Rosenbrock(x)
f = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
end

随后是梯度函数文件 grid.m 如下：

function g = grid(x)
g = [-400 * x(1) * x(2) + 400 * x(1)^3 + 2 * x(1) - 2;
    200 * x(2) - 200 * x(1)^2];
end

求步长的 Armjio 函数请看作者之前的这篇博客：采用非精确线搜索求步长的Armjio准则–MATLAB实现，编写函数文件 Armjio.m 如下，与刚才两个文件放在同一文件夹目录下：

function [alpha] = Armjio(fun, grid, x0, dk)
    %
    % Function [alpha, xk, fx, k] = Armjio(fun, grid, x0, dk)
    % 求出函数fun在x0处以dk为下降方向时的步长alpha，同时返回相对应的下
    % 一个下降点xk以及xk处的函数值fx，k为迭代次数
    % -----------------------------------------------------------
    % 输入: 
    % 		fun 	函数名称(字符变量）
    %		grid 	梯度函数名称(字符变量)
    %		x0		迭代点(列向量)
    %		dk		函数在迭代点处的下降方向(列向量)
    %
    % 输出:
    %		alpha	函数在x0处以dk为下降方向时的下降步长
    %		xk		函数在x0处以dk为下降方向，以alpha为步长
    %				求得的下降点
    %		fx		函数在下降点xk处的函数值
    %		k		求步长算法迭代次数
    % -----------------------------------------------------------
    % by Zhi Qiangfeng 
    %
    beta = 0.333; 		% 步长 alpha 的迭代系数，小于 1
    rho = 1e-3; 		% 泰勒展开式补足系数，0 < rho < 1/2
    alpha = 1; 			% 初始步长为 1
    k = 0; 				% 统计迭代次数
    gk = feval(grid, x0);	% x0处的梯度值
    fd = feval(fun, x0 + alpha * dk); 	% 函数在下一个迭代点处的目标函数值
    fk = feval(fun, x0) + alpha * rho * gk' * dk; 	% 函数在下一个迭代点处的泰勒展开值
    while fd > fk
        alpha = beta * alpha;
        fd = feval(fun, x0 + alpha * dk);
        fk = feval(fun, x0) + alpha * rho * gk' * dk;
        k = k + 1;
    end
end

以 Rosenbrock 函数为例，该函数的的极小点为[1; 1]，设置精度 eps = 1e-5，迭代最大次数 kmax = 3000，我们生成一个 2 行 10 列，范围在 [-10, 10] 的矩阵，每次选取其中的一列作为初始点，编写脚本文件如下：

X = randi([-10, 10], 1, 1) * rand(2, 10);
eps = 1e-5;
kmax = 3000;
file = fopen("./BFGSdata.txt", "w");
fprintf(file, "初始点\t\t\t\t\t 极小点\t\t\t\t  目标函数值\t\t 迭代次数\t 运行时间\n");
for i = 1:10
    tic 
    x0 = X(:, i);
    [f, xk, k] = BFGS(x0, "fun", "grid", eps, kmax);
    t = toc;
    fprintf(file, "[%f, %f]\t[%f, %f]\t%f\t\t%d\t\t\t%f\n", x0(1), x0(2), xk(1), xk(2), f, k, t);
end  

得到文本文件 BFGSdata.txt，内容如下：

初始点					极小点				  	目标函数值		迭代次数	 	运行时间
[0.119559, 0.469560]	[1.000000, 1.000000]	0.000000		22			0.002685
[0.706317, 1.642388]	[1.000000, 1.000000]	0.000000		15			0.002144
[0.030807, 0.086048]	[1.000000, 1.000000]	0.000000		27			0.002457
[0.337980, 1.298231]	[1.000000, 1.000000]	0.000000		20			0.001869
[1.463445, 1.295492]	[1.000001, 1.000001]	0.000000		21			0.002422
[0.901847, 1.094018]	[1.000000, 1.000000]	0.000000		14			0.001120
[0.592642, 1.489386]	[1.000000, 1.000000]	0.000000		19			0.001511
[0.377910, 1.373551]	[1.000000, 1.000000]	0.000000		20			0.001629
[0.367022, 0.736969]	[1.000000, 1.000000]	0.000000		18			0.001403
[1.251237, 1.560455]	[1.000000, 1.000000]	0.000000		14			0.001040

将文本文件数据导入 Excel 或 MATLAB 中，便可对该方法进行误差分析。

---

附

最优化相关算法设计数学原理：最优化/Optimization文章合集

---

有帮助可以点赞哦，谢谢大家的支持~