功能：拟牛顿法BFGS实例源码。
function x = Opt_BFGS(x0, Iter_max, eps)
% BFGS法确定多个变量的最优解
%
=======================================================================
% 输入：
% x0：当前点；
% Iter_max：最大迭代次数；
% eps：黄金分割法的精度
% 输出：
% x：最优解；
%
=======================================================================
ite = 1; % 循环次数
gk = BFGS_Gradient(x0); % 初始梯度
H = eye(length(x0)); % 初始H阵
d = -gk*H'; % 初始搜索方向
alpha = Advance_Retreat_Gold(x0, d, 0,
0.01, eps); % 初始搜索步长
while norm(gk) > eps && ite
< Iter_max
x0 = x0 + alpha*d; % 新的循环点
gkk = BFGS_Gradient(x0); % 新的梯度
y = gkk - gk; % 梯度差
s = alpha*d; % 迭代点差
H = H + (1 +
(y*H*y')/(s*y'))*(s'*s)/(s*y') - (s'*y*H + H*y'*s)/(s*y'); %
H阵更新
d = -gkk*H'; % 新的搜索方向
alpha = Advance_Retreat_Gold(x0, d, 0,
0.01, 1e-4); % 新的搜索步长
gk = gkk; % 更新梯度
ite = ite + 1; % 循环次数加1
end
x = x0 + alpha*d; % 极小值
end
% 梯度矩阵
function g = BFGS_Gradient(x0)
g = [2*x0(1) - 2, 2*x0(2) - 4, 2*x0(3) -
6, 2*x0(4) - 2, 2*x0(5) - 2, 2*x0(6) - 2];
end
% 一维搜索
function result = Advance_Retreat_Gold(x0,
d, t0, step, eps)
% 进退法和黄金分割法确定步长
%
=======================================================================
% 输入：
% x0：当前点；
% d：搜索方向；
% t0：进退法端点值；
% step：进退法步长；
% eps：黄金分割法的精度
% 输出：
% result：步长；
% 作者：王飞飞
% QQ：1105936347
%
=======================================================================
% 进退法确定大致区间
t1= t0 + step;
ft0 = BFGS_fun(x0, t0, d); % 目标函数
ft1 = BFGS_fun(x0, t1, d);
if ft1 <= ft0
step = 2*step;
t2 = t1 + step;
ft2 = BFGS_fun(x0, t2, d);
while ft1 > ft2
t1 = t2;
step = 2*step;
t2 = t1 + step;
ft1 = BFGS_fun(x0, t1, d);
ft2 = BFGS_fun(x0, 2, d);
end
else
step = step/2;
t2 = t1;
t1 = t2 - step;
ft1 = BFGS_fun(x0, t1, d);
while(ft1 > ft0)
step = step/2;
t2 = t1;
t1 = t2 - step;
ft1 = BFGS_fun(x0, t2, d);
end
end
a = 0;
b = t2;
% 黄金分割法确定精确解
a1 = a + 0.382*(b - a);
a2 = a + 0.618*(b - a);
f1 = BFGS_fun(x0, a1, d);
f2 = BFGS_fun(x0, a2, d);
while abs(b - a) >= eps
if f1 < f2
b = a2;
a2 = a1;
f2 = f1;
a1 = a + 0.382*(b - a);
f1 = BFGS_fun(x0, a1, d);
else
a = a1;
a1 = a2;
f1 = f2;
a2 = a + 0.618*(b - a);
f2 = BFGS_fun(x0, a2, d);
end
end;
result = 0.5*(a + b);
end
% 目标函数
function f = BFGS_fun(x0, r, d)
f = (x0(1) + r*d(1) - 1)^2 + (x0(2) +
r*d(2) - 2)^2 + (x0(3) + r*d(3) - 3)^2 + ...
(x0(4) + r*d(4) - 1)^2 + (x0(5) + r*d(5)
- 1)^2 + (x0(6) + r*d(6) - 1)^2;
end