%% basics of linear equations
% examples of linear equations
% invertible
M1 = [1 0; -1 1]
y = [1 2]'
% solution
inv(M1)*y
% non-invertible
M2 = [1 -1; -1 1]
inv(M2)*y % doesn't work
% invertible but possibly instable
M3 = [1 -1; -1 1.1]
inv(M3)*y
%% rank
rank(M1)
rank(M2)
rank(M3)
% What is the result of rank( zeros(1000) )?
% Does the rank change when multiplying the matrix with a scalar?
%% SVD for square matrix
M4 = [-0.5 0.5; 1 1];
[U,S] = svd( M4 )
% condition number
S(1,1)/S(2,2)
% change the similarity of columns and see what happens
M4 = [-0.1 0.1; 1 1];
[U,S] = svd( M4 )
% condition number
S(1,1)/S(2,2)
M4 = [-1 1; 1 1];
[U,S] = svd( M4 )
% condition number
S(1,1)/S(2,2)
% the vectors in U are orthonormal (orthogonal and unit length)
U'*U
%% SVD for rectangular matrix
M = [ [2 2 0]', [0 0 1]' ]
[U,S,V] = svd(M)
% check whether decomposition works out
U*S*V'
%% Why "eigen"vectors?
% when you stick eigenvectors into a covariance matrix...
(M*M')*U
% you get the vectors back, multiplied by square of eigenvalues
U*S.^2
% same for V
(M'*M)*V
S.^2*V
%% Why is this useful - SVD as a filter
% make some data
x = [-2*pi:0.01:2*pi];
data = [sin(x)' cos(x)'];
plot(data)
% replicate this matrix
data = repmat(data, 1, 10);
whos
% add noise
data = data + randn(size(data));
plot(data)
% 'economy size' svd
[U,S,V] = svd( data, 0 );
% plot eigenvalues
sdiag = diag(S)
figure
plot( sdiag, 'x' )
% if you want to plot variances
svar = 100 * sdiag.^2 / sum( sdiag.^2 )
plot(svar, 'x')
% leave only first two SVD components
sdiag2 = sdiag;
sdiag2(3:end) = 0
S2 = diag(sdiag2);
% reassamble filtered data
data2 = U*S2*V';
figure
plot(data2)
%% Why is this useful - SVD to compute pseudoinverse
% get M4 from above
M4 = [-0.1 0.1; 1 1];
[U,S,V] = svd( M4 )
% invert the eigenvalues
S(1,1) = 1/S(1,1)
S(2,2) = 1/S(2,2)
% compute the inverse
myinv = V*S*U'
% check whether it's correct
inv(M4)
% "dampen" the effect of small eigenvalues
S(2,2) = 5
% see how this affects the inverse
myinv2 = V*S*U'
% apply to some "data" vector y and note the difference
y = [1 1]'
myinv*y
myinv2*y
% the data that these solutions would produce
M4*myinv*y % exact
M4*myinv2*y % approximate
%% Relationship PCA and SVD
% The Matlab command princomp removes means from columns,
% which differs from svd()
% Therefore, let's start with a column-centred matrix
X = [1 1 -2; 1 -2 1; -2 1 1]
mean(X)
% what's the rank of this matrix?
% PCA
[A,B,C] = princomp(X)
% SVD
[U,S,V] = svd(X)
% note: a and v are the same, because for princomp rows are observations
% s and c should contain the same values, but c is half of s-squared
% reason: s-square is the "raw" variance, and princomp divideds by the
% degrees of freedom (here: 2)