function lbfgs(n) % L-BFGS, Algorithm 7.5 in N&W % n = problem dimension m = 5; x0 = zeros(n,1); xE = ones(n,1); % Exact Hk0 = eye(n); [fk, g0] = rosenbrock(x0); % linesearch parameters alpha_ini = 1.0; alpha_min = 1.0E-08; alpha_max = 1.0E+02; c1 = 1.0E-04; c2 = 0.9E-01; xtol = 1.0E-08; nfuneval = 20; quiet = 1; % Vector of vectors to store sk and yk s = zeros(m,n); y = zeros(m,n); Nmax = 1000; xk = x0; gk = g0; for k = 1:Nmax if norm(gk) <= 1.0E-12*norm(g0) % converged! break; end pk = twolooprec(Hk0, gk, m, s, y); % strong Wolfe linesearch [~,~,~,~,alpha,info,~] = minpack_cvsrch(@rosenbrock,xk, ... fk,gk,pk,alpha_ini, ... c1,c2,xtol, ... alpha_min,alpha_max, ... nfuneval,quiet); if info ~= 1 fprintf('Warning: linesearch routine returned error %d', ... info); end xk = xk + alpha*pk; [fk, gk1] = rosenbrock(xk); fprintf('%3d: f=%10.3e, |df|=%10.3e, |x-x*| = %10.3e\n', ... k, fk, norm(gk1), norm(xk-xE)); % Update storage s = [s(2:m,:); zeros(1,n)]; % Remove last and move others upwards y = [y(2:m,:); zeros(1,n)]; s(m,:) = alpha*pk'; % Store new values y(m,:) = gk1 - gk; gamma = (s(m,:)*y(m,:)') / (y(m,:)*y(m,:)'); Hk0 = gamma*eye(n); gk = gk1; end end function p = twolooprec(Hk0, gk, m, s, y) % Algorithm 7.4 q = gk; alpha = zeros(m,1); for i=m:-1:1 rho_inv = s(i,:)*y(i,:)'; if (rho_inv > 0) alpha(i) = (s(i,:)*q) / rho_inv; q = q - alpha(i)*y(i,:)'; end end r = Hk0*q; for i=1:m rho_inv = s(i,:)*y(i,:)'; if (rho_inv > 0) beta = (y(i,:)*r) / rho_inv; r = r + s(i,:)'*(alpha(i) - beta); end end p = -r; end % Backtracking algorithm function alpha=armijo(x,p) alpha = 1.0; rho = 0.5; c1 = 1.0E-03; [f0, df0] = rosenbrock(x); df0 = c1*df0'*p; while rosenbrock(x+alpha*p) > f0+alpha*df0 alpha = alpha*rho; end end