load matrix.dat

k = 3

X = matrix';
[d, N] = size(X)

%% part 1: plot the eigenvalues
% center the data, i.e. ensure that the mean in each dimension is zero
Xmean = mean(X, 2);
Xc = X - repmat(Xmean, 1, N);

% compute the sample covariance matrix
Sigmahat = (Xc * Xc')/N;

% compute the eigenvalues of the sample covariance matrix
[U, D] = eig(Sigmahat);

% extract the vector of eigenvalues from the returned matrix. we assume
% that these are sorted in ascending order
ev = diag(D);

% plot the eigenvalues
figure()
plot(ev)

%% part 2: compute Xhat and show images of X and Xhat
% truncate U to contain the eigenvectors corresponding to the k largest
% eigenvalues
Uk = U(:, end-k+1:end);

% compute the centered approximation
Xhatc = Uk * Uk' * Xc;

% restore the original mean
Xhat = Xhatc + repmat(Xmean, 1, N);

% compute the error
% directly using the frobenius norm
error_frobenius = (norm(Xhatc - Xc, 'fro')^2)/N
% and by summing the discarded eigenvalues
error_eigenvalues = sum(ev(1:end-k))

% draw the images
figure()
subplot(2,1,1)
image(X)
subplot(2,1,2)
image(Xhat)

% compression ratio is the size of the original data, divided by the size
% of the stored eigenvectors plus the projections of the documents
comp = d * N / (d*k + k*N)

%% part 3: find the angles between the projected documents
% it is better to compute the angles between the centered projections,
% however the notation used in the assignment implied computing the
% angles between the non-centered projections. both answers are computed
% here
for zzz = 1:2
    desc = {'Centered', 'Mean added'};
    disp(desc(zzz))

    for n = 1:N
    
        if zzz == 1
            a = Xhatc(:,1);
            b = Xhatc(:,n);
        else
            a = Xhat(:,1);
            b = Xhat(:,n);
        end
    
        dp = (a' * b) / (norm(a) * norm(b));
    
        rad = acos(dp);
        deg = rad * 180 / pi;
    
        disp(sprintf('%d: %g degrees, %g radians', n, deg, rad))
    
    end
end