function [a,b]=unrolling(s,y,b0)

% "Unrolling the BFGS formula" 
%
% This code is an implementation of Procedure 7.6 found in 
% J. Nocedal, S.J. Wright, Numerical Optimization, Springer-Verlag, New
% York, second edition, 2006.  

% Inputs: S, Y contain the L-BFGS updates 
% where B_0 = b0*I is the initial L-BFGS matrix

% The most recently computed s_i and y_i vectors are stored in the last
% columns of S and Y, respectively.


[m,n] = size(s);
a = zeros(m,n); 
b = zeros(m,n); 

for i = 1:n
  b(:,i) = y(:,i)/(sqrt(y(:,i)'*s(:,i)));
  a(:,i) = b0*s(:,i);
  for j = 1: (i-1)
    a(:,i) = a(:,i) + b(:,j)'*s(:,i)*b(:,j) - a(:,j)'*s(:,i)*a(:,j);
  end
  a(:,i) = a(:,i)/(sqrt(s(:,i)'*a(:,i)));
end


