function [old_dirs,old_stps,Hdiag,Bcompact] = lbfgsUpdate(y,s,corrections,debug,old_dirs,old_stps,Hdiag)

%B0 = eye(length(y))/Hdiag;
S = old_dirs(:,2:end);
Y = old_stps(:,2:end);
k = size(Y,2);
L = zeros(k);
for j = 1:k
    for i = j+1:k
        L(i,j) = S(:,i)'*Y(:,j);
    end
end
D = diag(diag(S'*Y));
N = [S/Hdiag Y];
M = [S'*S/Hdiag L;L' -D];

ys = y'*s;
Bs = s/Hdiag - N*(M\(N'*s)); % Product B*s
sBs = s'*Bs;

eta = .02;
if ys < eta*sBs
    if debug
        fprintf('Damped Update\n');
    end
    theta = min(max(0,((1-eta)*sBs)/(sBs - ys)),1);
    y = theta*y + (1-theta)*Bs;
end


numCorrections = size(old_dirs,2);
if numCorrections < corrections
    % Full Update
    old_dirs(:,numCorrections+1) = s;
    old_stps(:,numCorrections+1) = y;
else
    % Limited-Memory Update
    old_dirs = [old_dirs(:,2:corrections) s];
    old_stps = [old_stps(:,2:corrections) y];
end

% Update scale of initial Hessian approximation
Hdiag = (y'*s)/(y'*y);