%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a demo driver for recursionSolve_qN.m.  
% 
% For details, see
%
% "Limited-Memory BFGS Systems with Diagonal Updates" 
% by Jennifer B. Erway and Roummel F. Marcia
%
% which is available at 
%
% http://www.wfu.edu/~erwayjb/publications.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

%set parameters
M     =    5;   %number of L-BFGS updates
n     = 1000;   %dimension of problem
sigma =  0.1;   %sigma in (B + sigma*I)
epsil = 1e-8;   %lower bound on gamma*sigma

%form the matrices S and Y from quasi-Newton method
S = randn(n,M);
Y = randn(n,M);

%ensure B is positive definite
for i=1:M
  if S(:,i)'*Y(:,i)<0
    S(:,i)=-S(:,i);
  end
end

%set 1/gamma and ensure that gamma*sigma > epsil
gamma_inv = norm(Y(:,M))^2/(S(:,M)'*Y(:,M));
gamma_inv = min(gamma_inv, sigma/epsil);

%right-hand side
z = randn(n,1);

%unroll the L-BFGS formula
[a,b] = unrolling(S,Y,gamma_inv);

tic;
%solve (B + sigma*I)x=z
[x] = recursionSolve_qN(a,b,gamma_inv,sigma,z);
time_recur = toc;

%solve using direct method
B = gamma_inv*eye(n);
for i=1:M
    B = B - a(:,i)*a(:,i)'+ b(:,i)*b(:,i)';
end
tic;
t = (B+sigma*eye(n))\z;
time_direct = toc;

%test recursion solution quality by computing the
%residual relative error ||(B + sigma*I)*x - z||/||z||
resid = norm( B*x+sigma*x - z)/norm(z);
fprintf('\n');
fprintf('Residual relative error: %8.2e\n\n',     resid);
fprintf('Recursion solve time:    %8.6f sec\n',   time_recur);
fprintf('Direct solve time:       %8.6f sec\n\n', time_direct);


