%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a demo driver for gen_recursionSolve_qN.m.
%
% For details, see
%
% "Solving limited-memory BFGS systems with
% generalized diagonal upates"
% by Jennifer Erway and Roummel 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
epsil = 1e-6; %lower bound on gamma*sigma
%generate a positive-definite diagonal matrix
D = rand(n,1);
%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
%ensure min(D)>=sigma
%set 1/gamma and ensure that gamma*sigma > epsil
minD = min(D);
gamma_inv = norm(Y(:,M))^2/(S(:,M)'*Y(:,M));
gamma_inv = min(gamma_inv, minD/epsil);
%right-hand side
z = randn(n,1);
%unroll the L-BFGS formula
[a,b] = unrolling(S,Y,gamma_inv);
tic;
%solve (B + D)x=z
[x] = gen_recursionSolve_qN(a,b,gamma_inv,D,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+diag(D))\z;
time_direct = toc;
%test recursion solution quality by computing the
%residual relative error ||(B + D)*x - z||/||z||
resid = norm( B*x+D.*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);