function [L, dw, x_out, energy_evals] = back_lbfgs(energy,loss,x0,w,N,m,stepsize)

% evaluate the gradient of a function implicitly defined in terms of 
% an LBFGS optimization
%
% see paper "generic methods for optimization based modeling", AISTATS 2012
%
% author: justin domke
% contact (my first name).(my last name)@nicta.com.au
% history: Written in late 2011
%          Published Jan 31, 2014

% by default, run for 10 iterations 
if nargin < 5
    N = 10;
end

% by default, store 100 gradients
if nargin < 6
    m = min(N,100);
end

energy_evals = 0;

% this is a list of step sizes
if nargin < 7
    stepsize = .1;
end
lambdas = stepsize*ones(N,1);

% LBFGS constant
gamma0 = 1;

D = length(x0);

x = zeros(D,N+1);
g = zeros(D,N);
r = zeros(D,N);
s = zeros(D,N);
y = zeros(D,N);
p = zeros(1,N);
a = zeros(1,N);
q = zeros(D,N);

dx = zeros(D,N+1);
dg = zeros(D,N);
dr = zeros(D,N);
ds = zeros(D,N);
dy = zeros(D,N);
dp = zeros(1,N);
da = zeros(1,N);
dq = zeros(D,N);

x(:,1)=x0;

for k=1:N
    [~, g(:,k)] = energy(x(:,k),w);
    
    energy_evals = energy_evals + 1;
    if k>1 && norm(g(:,k))<1e-10
        N=k-1;
        break;
    end
    
    mp = min(m,k-1);
    
    q(:,k) = g(:,k);
    for i=k-1:-1:k-mp
        s(:,i) = x(:,i+1) - x(:,i);
        y(:,i) = g(:,i+1) - g(:,i);
        p(i)   = 1/(y(:,i)'*s(:,i));
        a(i)   = p(i)*(s(:,i)'*q(:,i+1));
        q(:,i) = q(:,i+1) - a(i)*y(:,i);
    end
    
    if k > 1,
        gamma = (s(:,k-1)'*y(:,k-1))/(y(:,k-1)'*y(:,k-1));
    else
        gamma = gamma0;
    end
    
    r(:,k-mp) = gamma*q(:,k-mp);
    for i=k-mp:k-1
        r(:,i+1) = r(:,i) + s(:,i)*(a(i)-p(i)*(y(:,i)'*r(:,i)));
    end
    
    x(:,k+1) = x(:,k) - lambdas(k)*r(:,k);
end

dw = 0*w;

if isempty(loss), L=0; dw=0*w; x_out = x(:,N+1); return; end;
[L dx(:,N+1)] = loss(x(:,N+1));
if nargout==1, return; end

for k=N:-1:1
    mp = min(m,k-1);
    
    dr(:,k) = -lambdas(k)*dx(:,k+1);
    dx(:,k) = dx(:,k) + dx(:,k+1);
    
    q(:,k) = g(:,k);
    for i=k-1:-1:k-mp
        s(:,i) = x(:,i+1) - x(:,i);
        y(:,i) = g(:,i+1) - g(:,i);
        p(i)   = 1/(y(:,i)'*s(:,i));
        a(i)   = p(i)*(s(:,i)'*q(:,i+1));
        q(:,i) = q(:,i+1) - a(i)*y(:,i);
    end
    
    if k > 1,
        gamma = (s(:,k-1)'*y(:,k-1))/(y(:,k-1)'*y(:,k-1));
    else
        gamma = gamma0;
    end
    
    r(:,k-mp) = gamma*q(:,k-mp);
    for i=k-mp:k-1
        r(:,i+1) = r(:,i) + s(:,i)*(a(i)-p(i)*(y(:,i)'*r(:,i)));
    end
    for i=k-1:-1:k-mp
        s_dr = s(:,i)'*dr(:,i+1);
        y_r  = y(:,i)'*r(:,i);
        dr(:,i) = (1-p(i)*(y(:,i)'*s(:,i)))*dr(:,i+1);
        dr(:,i) = dr(:,i+1)-p(i)*s_dr*y(:,i);
        ds(:,i) = (a(i)-p(:,i)*y_r) * dr(:,i+1);
        da(:,i) =  s(:,i)' * dr(:,i+1);
        dy(:,i) = -p(i) * s_dr * r(:,i);
        dp(:,i) = -y_r * s_dr;

    end
    dq(:,k-mp) = gamma*dr(:,k-mp);
    
    if k>1
        dgamma = q(:,k-mp)'*dr(:,k-mp);
        y_y = y(:,k-1)'*y(:,k-1);
        s_y = s(:,k-1)'*y(:,k-1);
        ds(:,k-1) = ds(:,k-1) + y(:,k-1)*(dgamma/y_y);
        dy(:,k-1) = dy(:,k-1) + s(:,k-1)*(dgamma/y_y) - y(:,k-1)*(2*dgamma*s_y/y_y^2);
    end
    
    for i=k-mp:k-1
        y_s = y(:,i)'*s(:,i);
        dq(:,i+1) = dq(:,i);
        da(i)     = da(i)   - (y(:,i)'*dq(:,i));
        dy(:,i)   = dy(:,i) - a(i) *dq(:,i);
        dp(i)     = dp(i)   + (s(:,i)'*q(:,i+1))*da(i);
        ds(:,i)   = ds(:,i)   + p(i)*da(:,i)*q(:,i+1);
        dq(:,i+1) = dq(:,i+1) + p(i)*da(:,i)*s(:,i);
        dy(:,i)   = dy(:,i)   - s(:,i)*(dp(i)/y_s^2);
        ds(:,i)   = ds(:,i)   - y(:,i)*(dp(i)/y_s^2);
        dg(:,i)   = dg(:,i)   - dy(:,i);
        dg(:,i+1) = dg(:,i+1) + dy(:,i);
        dx(:,i)   = dx(:,i)   - ds(:,i);
        dx(:,i+1) = dx(:,i+1) + ds(:,i);
    end
    dg(:,k) = dg(:,k) + dq(:,k);
    
    % get stepsize
    e = (eps)^(1/3)*(1+norm(x(:,k),inf))/norm(dg(:,k),inf);
    
    % perform perturbations
    x2 = x(:,k) + e*dg(:,k);
    x3 = x(:,k) - e*dg(:,k);
    [~,dEdx2,dEdw2] = energy(x2  ,w);
    [~,dEdx3,dEdw3] = energy(x3  ,w);
    energy_evals = energy_evals + 2;
    
    dx(:,k) = dx(:,k) + (1/e/2)*(dEdx2-dEdx3);
    dw      = dw      + (1/e/2)*(dEdw2-dEdw3);
end

x_out = x(:,end);

end
