1 %% basics of linear equations
   2 
   3 % examples of linear equations
   4 
   5 % invertible
   6 M1 = [1 0; -1 1]
   7 y = [1 2]'
   8 % solution
   9 inv(M1)*y
  10 
  11 % non-invertible
  12 M2 = [1 -1; -1 1]
  13 inv(M2)*y  % doesn't work
  14 
  15 % invertible but possibly instable
  16 M3 = [1 -1; -1 1.1]
  17 inv(M3)*y
  18 
  19 % What does that mean for the effect of noise?
  20 n = randn(2,10); % matrix with random noise
  21 inv(M1)*n
  22 inv(M3)*n
  23 
  24 
  25 %% rank
  26 rank(M1)
  27 rank(M2)
  28 rank(M3)
  29 
  30 % What is the result of rank( zeros(1000) )? rank( eye(5) )?
  31 % Does the rank change when multiplying the matrix with a scalar?
  32 
  33 %% Basic SVD
  34 
  35 % get singular values for above matrices, compute "condition numbers"
  36 
  37 s(1)/s(2)
  38 
  39 s = svd(M2)
  40 s(1)/s(2)
  41 
  42 s = svd(M3)
  43 s(1)/s(2)
  44 
  45 %% SVD for square matrix
  46 
  47 M = eye(2)
  48 [U,S] = svd(M)
  49 
  50 [U,S] = svd(M1)
  51 
  52 [U,S] = svd(M2)
  53 
  54 [U,S] = svd(M3)
  55 
  56 % the vectors in U are orthonormal (orthogonal and unit length)
  57 U'*U
  58 
  59 %% SVD for rectangular matrix
  60 M = [ [2 2 0]', [0 0 1]' ]
  61 [U,S,V] = svd(M)
  62 
  63 % check whether decomposition works out
  64 U*S*V'
  65 
  66 
  67 %% Why is this useful - SVD as a filter
  68 % make some data
  69 x = [-2*pi:0.01:2*pi];
  70 data = [sin(x)' cos(x)'];
  71 
  72 plot(data)
  73 
  74 % replicate this matrix
  75 data = repmat(data, 1, 10);
  76 whos
  77 
  78 % add noise
  79 data = data + randn(size(data));
  80 plot(data)
  81 
  82 % 'economy size' svd
  83 [U,S,V] = svd( data, 0 );
  84 
  85 % plot eigenvalues
  86 sdiag = diag(S)
  87 figure
  88 plot( sdiag, 'x' )
  89 % if you want to plot variances
  90 svar = 100 * sdiag.^2 / sum( sdiag.^2 )
  91 plot(svar, 'x')
  92 
  93 % leave only first two SVD components
  94 sdiag2 = sdiag;
  95 sdiag2(3:end) = 0
  96 S2 = diag(sdiag2);
  97 
  98 % reassamble filtered data
  99 data2 = U*S2*V';
 100 figure
 101 plot(data2)
 102 
 103 %% Relationship PCA and SVD
 104 
 105 % The Matlab command princomp removes means from columns,
 106 % which differs from svd()
 107 % Therefore, let's start with a column-centred matrix
 108 
 109 X = [1 1 -2; 1 -2 1; -2 1 1]
 110 mean(X)
 111 % what's the rank of this matrix?
 112 
 113 % PCA
 114 [A,B,C] = princomp(X)
 115 
 116 % SVD
 117 [U,S,V] = svd(X)
 118 
 119 % note: a and v are the same, because for princomp rows are observations
 120 
 121 % s and c should contain the same values, but c is half of s-squared
 122 % reason: s-square is the "raw" variance, and princomp divides by the
 123 % degrees of freedom (here: 2)
 124 
 125 %% Why "eigen"vectors?
 126 % when you stick eigenvectors into a covariance matrix...
 127 (M*M')*U
 128 
 129 % you get the vectors back, multiplied by square of eigenvalues
 130 U*S.^2
 131 
 132 % same for V
 133 (M'*M)*V
 134 S.^2*V
 135 
 136 %% SVD to compute pseudoinverse
 137 
 138 % get M4 from above
 139 M4 = [-0.1 0.1; 1 1];
 140 [U,S,V] = svd( M4 )
 141 
 142 % invert the eigenvalues
 143 S(1,1) = 1/S(1,1)
 144 S(2,2) = 1/S(2,2)
 145 
 146 % compute the inverse
 147 myinv = V*S*U'
 148 
 149 % check whether it's correct
 150 inv(M4)
 151 
 152 % "dampen" the effect of small eigenvalues
 153 S(2,2) = 5
 154 
 155 % see how this affects the inverse
 156 myinv2 = V*S*U'
 157 
 158 % apply to some "data" vector y and note the difference
 159 y = [1 1]'
 160 
 161 myinv*y
 162 myinv2*y
 163 
 164 % the data that these solutions would produce
 165 M4*myinv*y      % exact
 166 M4*myinv2*y     % approximate