% Faces
% Download faces from http://www.cs.nyu.edu/~roweis/data/olivettifaces.mat
load olivettifaces.mat;
imy = 64;
imx = 64;
n = size(faces,2);
X = faces/255;

montage(reshape(X,[imy imx 1 n]),'DisplayRange',[]);

% SVD
[U,S,V] = svd(X);

figure(1); clf;
plot(diag(S));

figure(2); clf;
k = 10;
montage(reshape(U(:,1:k),[imy imx 1 k]),'DisplayRange',[]);

% Show best rank-k approximation of image set
Xr = U(:,1:k) * S(1:k,1:k) * V(:,1:k)';
montage(reshape(Xr,[imy imx 1 n]));

% Show reconstructions of the ith image
i  = 1;
Xr = U(:,1)*S(1,1)*V(i,1);
for k = 2:n,
    Xr(:,k) = Xr(:,k-1) + U(:,k)*S(k,k)*V(i,k);
end
k = 50;
figure(1); clf;
montage(reshape(Xr(:,1:k),[imy imx 1 k]));
figure(2); clf;
montage(reshape(U(:,1:k),[imy imx 1 k]),'DisplayRange',[]);
