function hw7p2

rng(5)

%%% Set various parameters.
m = 200;
b = 70;
q = 2;

%%% Build a symmetric test matrix
A = fast_decay(m);

%%% Run power iteration to build an ON basis Q1 for the range of A.
Y = randn(m,b);
for k = 1:q
  Y = A*Y;
end
[Q1,~] = qr(Y,0);
    
%%% Compute the norm of the difference between A and the projection of A
%%% onto the first k columns of Q1.
E = zeros(1,b);
for k = 1:b
  E(k) = norm(A - Q1(:,1:k)*Q1(:,1:k)'*A);
end

%%% Plot the results
ss = svd(A);
semilogy(0:b,ss(1:(b+1)),'k',...
         1:b,E,'b','LineWidth',2)
     
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function builds a matrix A that is symmetric and whose eigenvalues
% decay rapidly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = fast_decay(m)

xxs = rand(2,m);
xxt = [1.8+rand(1,m);rand(1,m)];
DD1 = xxs(1,:)'*ones(1,m) - ones(m,1)*xxt(1,:);
DD2 = xxs(2,:)'*ones(1,m) - ones(m,1)*xxt(2,:);
A   = log(DD1.*DD1 + DD2.*DD2);
A   = A + A';
A   = (1/norm(A))*A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
