最优化——拟牛顿方法matlab程序

% BFGS

function [x, output] = bfgs(fun, dfun, x0, varargin)

% Step 1: initialization

epsi = 1.0e-6;

k = 0;

funcN = 0;

rho = 0.01; l = 0.15; u = 0.85;

x = x0;

f = feval(fun, x, varargin{:});

funcN = funcN + 1;

n = length(x0);

H = eye(n);

% Step 2: check termination condition

g = feval(dfun, x, varargin{:});

while norm(g) > epsi & k <= 100

itercon = true;

d = -H*g;

% Step 3: line search

alpha_0 = 1.0;

gd = g'*d;

[alpha, funcNk, exitflag] = ...

lines(fun, rho, l, u, alpha_0, f, gd, x, d, varargin{:});

funcN = funcN + funcNk;

if exitflag == -1

itercon = false;

restart = true;

H = eye(n);

gold = g;

end

% Step 4: compute new point

if itercon

s = alpha * d;

x = x + s;

f = feval(fun, x, varargin{:});

funcN = funcN + 1;

gold = g;

g = feval(dfun, x, varargin{:});

% Stept 5: update H

y = g - gold;

hy = H*y; sy = s'*y; yhy = y'*hy;

if sy < 0.2*yhy

theta = 0.8*yhy/(yhy-sy);

s = theta*s + (1.0-theta)*hy;

sy = 0.2*yhy;

end

v = sqrt(yhy)*(s/sy-hy/yhy);

H = H + s*s'/sy-hy*hy'/yhy+v*v';

end

k = k + 1;

end

% Step 6: output

output.fval = f;

output.iteration = k;

output.funcount = funcN;

output.gnorm = norm(g);

end

测试

% f1.m

function f = f1(x)

% objective function

f = (x(2)-x(1)^2)^2+(1-x(1))^2;

end

% df1.m

function g = df1(x)

% grads function

g = [4.0*(x(1)^3-x(1)*x(2))+2*x(1)-2;2.0*(x(2)-x(1)^2)];

end

% Command

>> x0=[-1.9;2];

>> [x, output]=bfgs('f1', 'df1', x0)

结果

x =

1.0000

1.0000

output =

fval: 6.3098e-017

iteration: 12

funcount: 28

gnorm: 3.8353e-008