采用Armjio准则求步长的BFGS/DFP拟牛顿方法--MATLAB实现
Posted Z.Q.Feng
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了采用Armjio准则求步长的BFGS/DFP拟牛顿方法--MATLAB实现相关的知识,希望对你有一定的参考价值。
实现原理
具体数学实现原理可参考这篇文章:拟牛顿法/Quasi-Newton,DFP算法/Davidon-Fletcher-Powell,及BFGS算法/Broyden-Fletcher-Goldfarb-Shanno
代码实现
拟牛顿方法
算法伪码如下:
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
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文章合集
有帮助可以点赞哦,谢谢大家的支持~
以上是关于采用Armjio准则求步长的BFGS/DFP拟牛顿方法--MATLAB实现的主要内容,如果未能解决你的问题,请参考以下文章
采用非精确线搜索求步长的Armjio准则--MATLAB实现
采用非精确线搜索求步长的Armjio准则--MATLAB实现
采用Armjio非精确线搜索准则的最速下降法--MATLAB实现
采用Armjio非精确线搜索准则的最速下降法--MATLAB实现