%rank k approximation of symmetric image
%copyright D. Aristoff :)

A = zeros(100,100); %100x100 square image   

for i=1:100      %set image = symmetric circular band
    for j=1:100
        r = sqrt((i-50)^2+(j-50)^2);
        if r > 25 && r < 31
            A(i,j) = 1;
        end
    end
end

for i=1:100        %introduce noise into image
    for j=1:i
        if rand() < 0.5             %alter half the pixels randomly
            A(i,j) = mod(A(i,j)+rand,1);   %introduce noise into the image
            A(j,i) = A(i,j);        %make the noise symmetric...
        end
    end
end

[Q,D] = eig(A);                   %This diagonalizes A, so D = Q^T*A*Q
lambdas = transpose(diag(D));     %lambdas = vector of eigenvalues of A
[lam,perm] = sort(abs(lambdas),'descend');  %put eigenvalues in descending order
Q = Q(:,perm);                    %reorder eigenvectors accordingly
lambdas = lambdas(perm);          %reorder eigenvalues accordingly

k = 1;
while k > 0

k = input('Compute rank k approximation where k = ... ');
close all

Dk = diag([lambdas(1:k),zeros(1,100-k)]);  %Dk consists of k largest eigenvalues
Ak = Q*Dk*transpose(Q);                    %Ak is a rank k approximation of A

subplot(2,1,1)
imagesc(A)
axis square
title('original image')
colorbar

subplot(2,1,2)
imagesc(Ak)
axis square
title('rank k approximation')
colorbar

movegui(gcf,'northwest')

end