% Iterative methods for isolated eigenpairs. Barnett 3/25/14-3/27/14

rng(1);         % initialize random number generator so repeatable
n = 100; A = randn(n,n); A = A+A';  % given symm real matrix
figure; imagesc(A); axis equal; colorbar; title('A');

% ---------------- Power iteration
x = randn(n,1);  % random starting vector
for k=1:1000
  x = A*x;
  x = x/norm(x);
  disp(x(1:3)') % look at first 3 components of eigenvector converging
end
% note code looks like the math; strive for this.

l = eig(A)         % solve for evs via our best algorithms!

[~,i] = sort(abs(l),'descend');  % get indices of ordering of size of evs
r = abs(l(i(2))/l(i(1)))   % predicted rate, can be terribly close to 1 (rerun)




% ------------ Inverse power iteration
x = randn(n,1);  % random starting vector
mu = -30;        % estimate for an eigenvalue
for k=1:20
  w = (A - mu*eye(n)) \ x;   % linear system solve each iter
  x = w/norm(w);
  disp(x(1:3)')             % much faster convergence
  %R(x)
end




% ---------- estimating eigenvalues
R = @(x) x'*A*x / norm(x)^2;  % Rayleigh quot - note use of inline func
% now insert R(x) into above



% ------------ Rayleigh quotient iteration
x = randn(n,1);  % random starting vector
mu = -29;  % estimate for an eigenvalue
for k=1:5
  w = (A - mu*eye(n)) \ x;   % solve
  x = w/norm(w);
  disp(x(1:3)') % look at first 3 components of eigenvector converging
  mu = R(x);  % update to best avail mu
end
ls(i) = mu;
% all digits in 5 iters! Note (A-mu.I) becomes singular, but not a problem...
end




% ------------ Damian suggested looping over starting mu's (takes 1 min):
mus = -50:0.01:50; ls = nan*mus; for i=1:numel(mus), mu = mus(i);
  x = randn(n,1);  % random starting vector
  for k=1:5
    w = (A - mu*eye(n)) \ x;   % solve
    x = w/norm(w);
    mu = R(x);  % update to best avail mu
  end
  ls(i) = mu;
end
figure; plot(mus,ls,'.'); xlabel('starting \mu'); ylabel('$\hat{\lambda_m}$','interpreter','latex');
