www.pudn.com > starter.zip > lbfgs.m, change:2011-01-04,size:924b

function [d] = lbfgs(g,s,y,Hdiag) % BFGS Search Direction % % This function returns the (L-BFGS) approximate inverse Hessian, % multiplied by the gradient % % If you pass in all previous directions/sizes, it will be the same as full BFGS % If you truncate to the k most recent directions/sizes, it will be L-BFGS % % s - previous search directions (p by k) % y - previous step sizes (p by k) % g - gradient (p by 1) % Hdiag - value of initial Hessian diagonal elements (scalar) [p,k] = size(s); for i = 1:k ro(i,1) = 1/(y(:,i)'*s(:,i)); end q = zeros(p,k+1); r = zeros(p,k+1); al =zeros(k,1); be =zeros(k,1); q(:,k+1) = g; for i = k:-1:1 al(i) = ro(i)*s(:,i)'*q(:,i+1); q(:,i) = q(:,i+1)-al(i)*y(:,i); end % Multiply by Initial Hessian r(:,1) = Hdiag*q(:,1); for i = 1:k be(i) = ro(i)*y(:,i)'*r(:,i); r(:,i+1) = r(:,i) + s(:,i)*(al(i)-be(i)); end d=r(:,k+1);