%% Chapter 7 Eigenvalue Problems
% _This is a MATLAB Live Editor version of the *Fundamentals of Numerical Computation*._ 
% Example
% Here is the adjacency matrix of a graph with 60 nodes from a built-in MATLAB 
% function.

clear
[A,v] = bucky;
size(A)
%% 
% The extra vector $\mathbf{v}$ gives particular coordinates for each node on 
% the unit sphere. Plotting the nodes and edges with these coordinates reveals 
% beautiful structure.

gplot(A,v), axis equal
%% 
% As you can see, the edges resemble those on a soccer ball. The same structure, 
% when made of carbon atoms at the nodes, is called a _buckyball_ in materials 
% science.
% Example
% MATLAB reads an image using the |imread| command. We'll use an image that 
% ships with MATLAB for demonstrations.

A = imread('peppers.png');
size(A)
%% 
% As stated above, the image has three ``layers'' for red, green, and blue. 
% We can deal with each layer as a matrix, or (as below) convert it to a single 
% matrix indicating shades of gray from black (0) to white (255). Then we have 
% to tell MATLAB to represent the entries as floating-point values rather than 
% integers. 

A = rgb2gray(A); % collapse from 3 dimensions to 2
A = double(A);   % convert to floating point
[m,n] = size(A)
%% 
% We now have a matrix. To show it as an image, we use the |imshow| command. 
% We can tell it explicitly which values correspond to black and white.

imshow(A,[0 255])
% Example
% The |eig| command will give just the eigenvalues if one output is given, or 
% both $\mathbf{V}$ and $\mathbf{D}$ if two are given.

A = pi*ones(2,2);
lambda = eig(A)
%%
[V,D] = eig(A)
%% 
% We can check the fact that this is an EVD.

norm( A - V*D/V )   % /V is like *inv(V)
%% 
% Even if the matrix is not diagonalizable, |eig| will run successfully, but 
% the matrix $\mathbf{V}$ will not be invertible.

[V,D] = eig([1 1;0 1])
rankV = rank(V)
%% 
% We will confirm the Bauer-Fike theorem on a triangular matrix. These tend 
% to be far from normal. 

n = 15;
lambda = (1:n)';
A = triu( ones(n,1)*lambda' );
% Example
% The Bauer-Fike theorem provides an upper bound on the condition number of 
% these eigenvalues.

[V,D] = eig(A); 
kappa = cond(V)
%% 
% The theorem suggests that eigenvalue changes may be up to 7 orders of magnitude 
% larger than a perturbation to the matrix. A few random experiments show that 
% effects of nearly that size are not hard to observe.

for k = 1:3
    E = randn(n);  E = 1e-7*E/norm(E);
    mu = eig(A+E);
    max_change = norm( sort(mu)-lambda, inf )
end
% Example
% We verify some of the fundamental SVD properties using the built-in MATLAB 
% command |svd|.

A = vander(1:5);
A = A(:,1:4)
%%
[U,S,V] = svd(A);
norm(U'*U - eye(5))
%%
norm(V'*V - eye(4))
%%
sigma = diag(S)
%%
[ norm(A) sigma(1) ]
%%
[ cond(A) sigma(1)/sigma(4) ]
% Example
% The following matrix is not hermitian.

A = [0 2; -2 0]
%% 
% It has an eigenvalue decomposition with a unitary matrix of eigenvectors, 
% though, so it is normal. 

[V,D] = eig(A);
norm( V'*V - eye(2) )
%% 
% The eigenvalues are pure imaginary.

lambda = diag(D)
%% 
% The singular values are the complex magnitudes of the eigenvalues.

svd(A)
% Example
% We construct a real symmetric matrix with known eigenvalues by using the QR 
% factorization to produce a random orthogonal set of eigenvectors. 

n = 30;
lambda = (1:n)';  
D = diag(lambda);
[V,R] = qr(randn(n));   % get a random orthogonal V
A = V*D*V';
%% 
% The condition number of these eigenvalues is one. Thus the effect on them 
% is bounded by the norm of the perturbation to $\mathbf{A}$. 

for k = 1:3
    E = randn(n); E = 1e-4*E/norm(E);
    mu = sort(eig(A+E));
    max_change = norm(mu-lambda,inf)
end
% Example
% Let's start with a known set of eigenvalues and an orthogonal eigenvector 
% basis.

D = diag([-6 -1 2 4 5]);
[V,R] = qr(randn(5));
A = V*D*V';    % note that V' = inv(V)
%% 
% Now we will take the QR factorization and just reverse the factors.

[Q,R] = qr(A);
A = R*Q;
%% 
% It turns out that this is a similarity transformation, so the eigenvalues 
% are unchanged.

eig(A)
%% 
% What's remarkable is that if we repeat the transformation many times, the 
% process converges to $\mathbf{D}$. 

for k = 1:15
    [Q,R] = qr(A);
    A = R*Q;
end
A
% Example
% We construct a symmetric matrix with a known EVD.

n = 20;
lambda = (1:n)';  D = diag(lambda);
[V,~] = qr(randn(n));   % get a random orthogonal V
A = V*D*V';
%% 
% The Rayleigh quotient of an eigenvector is its eigenvalue.

R = @(x) (x'*A*x)/(x'*x);
format long, R(V(:,7))
%% 
% The Rayleigh quotient's value is much closer to an eigenvalue than its input 
% is to an eigenvector. In this experiment, each additional digit of accuracy 
% in the eigenvector estimate gives two more digits to the eigenvalue estimate.

delta = 1./10.^(1:4)';
quotient = 0*delta;
for k = 1:4
    e = randn(n,1);  e = delta(k)*e/norm(e);
    x = V(:,7)+e;
    quotient(k) = R(x);
end
table(delta,quotient)
% Example
% We make an image from some text, then reload it as a matrix.

figure
tobj = text(0,0,'Hello world','fontsize',44);
ex = get(tobj,'extent');      % ignore this line 
axis([ex(1) ex(1)+ex(3) ex(2) ex(2)+ex(4)]), axis off   % ignore this line 
saveas(gcf,'hello.png')
A = imread('hello.png');
A = double(rgb2gray(A));
imagesc(A),  colormap gray
[m,n] = size(A)
%% 
% Next we show that the singular values decrease exponentially, until they reach 
% zero (more precisely, are about $\sigma_1 \varepsilon_{\rm mach}$). For all 
% numerical purposes, this determines the rank of the matrix.

[U,S,V] = svd(A);
sigma = diag(S);
semilogy(sigma,'.')
title('singular values'), axis tight        % ignore this line 
xlabel('i'), ylabel('\sigma_i')  % ignore this line 
r = find(sigma/sigma(1) > 10*eps,1,'last')
%% 
% The rapid decrease suggests that we can get fairly good low-rank approximations. 

figure
for i = 1:4
    subplot(2,2,i)
    k = 2*i;
    Ak = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';
    imshow(Ak,[0 255])
    title(sprintf('rank = %d',k))
end
%% 
% Consider how little data is needed to reconstruct these images. For rank 8, 
% for instance, we have 8 left and right singular vectors plus 8 singular values, 
% for a compression ratio of better than 25:1. 

compression = 8*(m+n+1) / (m*n)
% Example
% This matrix describes the votes on bills in the 111th session of the United 
% States Senate. (The data set was obtained from |voteview.com|.) Each row is 
% one senator and each column is a vote item.

load voting    % get from the book's website
[m,n] = size(A);
%% 
% If we visualize the votes (white is ``yea,'' black is ``nay,'' and gray is 
% anything else), we can see great similarity between many rows, reflecting party 
% unity.

figure
imagesc(A)
colormap gray
title('Votes in 111th U.S. Senate')   % ignore this line
ylabel('senator'), xlabel('bill')    % ignore this line
%% 
% We use singular value ``energy'' to quantify the decay rate of the values. 

[U,S,V] = svd(A);
sigma = diag(S);
tau = cumsum(sigma.^2) / sum(sigma.^2);
plot(tau(1:16),'.')
xlabel('k'), ylabel('\tau_k')  % ignore this line
title('Fraction of singular value energy')     % ignore this line
%% 
% The first and second singular triples contain about 58% and 17% respectively 
% of the energy of the matrix. All others have far less effect, suggesting that 
% the information is primarily two-dimensional. The first left and right singular 
% vectors also contain interesting structure.

subplot(211), plot(U(:,1),'.')
axis tight  % ignore this line
xlabel('senator number'), title('left singular vector')  % ignore this line
subplot(212), plot(V(:,1),'.')
axis tight  % ignore this line
xlabel('bill number'), title('right singular vector')  % ignore this line
%% 
% Both vectors have values greatly clustered near $\pm C$ for a constant $C$. 
% These can be roughly interpreted as how partisan a particular senator or bill 
% was, and for which political party. 
% 
% Projecting the senators' vectors into the first two $\mathbf{V}$-coordinates 
% gives a particularly nice way to reduce them to two dimensions. Political scientists 
% label these dimensions ``partisanship'' and ``bipartisanship.'' Here we color 
% them by actual party affiliation (also given in the data file): red for Republican, 
% blue for Democrat, and black for independent. 

clf
x1 = V(:,1)'*A';   x2 = V(:,2)'*A'; 
scatter(x1(Dem),x2(Dem),20,'b'),  hold on
scatter(x1(Rep),x2(Rep),20,'r')
scatter(x1(Ind),x2(Ind),20,'k')
xlabel('partisanship'),  ylabel('bipartisanship')  % ignore this line
title('111th US Senate in 2D')