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

A = zeros(150,100);  %150x100 pixel image   

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

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

[U,S,V] = svd(A);                    %This is the SVD of A, so A = U*S*V^T
sigs = transpose(diag(S));           %lambdas = vector of eigenvalues of A
[sigs,perm] = sort(sigs,'descend');  %put eigenvalues in descending order
V = V(:,perm);                       %reorder eigenvectors of A^T*A accordingly
perm = [perm,linspace(max(perm)+1,150,150-max(perm))];
U = U(:,perm);                       %reorder eigenvectors of A*A^T accordingly

k = 1;
while k > 0

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

Dk = diag([sigs(1:k),zeros(1,100-k)]);  %Dk consists of k largest singular values
Dk = [Dk;zeros(50,100)];                %Augment Dk so that it has the same size as A
Ak = U*Dk*transpose(V);                 %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