1 function [B,evalues] = PCA_fixspace_evalues(X,XFix,A,k,method,epsilon)
2 %
function [B,evalues] = PCA_fixspace(X,XFix [, [A], [k], [method],[epsilon] ])
4 %
function computing a PCA basis of the set of column vectors X projected on
5 % the orthogonal complement of span(XFix) The computation is performed
6 % explicitly using the gram matrix.
9 % X: the matrix of column vectors on which the PCA shall be applied
10 % XFix: matrix of column vectors spanning a fix space to which the vectors
11 % in 'X' shall be orthonormalized before the PCA.
12 % A: inner product matrix defining the scalar product `\langle x, y
13 % \rangle = y^t A x`. (default='1')
14 % k: maximum number of principal components to be computed. (default=
15 % number of columns in 'X').
16 % method:
string specifying the orthonormalization method, either 'qr' or
17 % 'gram-schmidt'. (default='gram-schmidt')
18 % epsilon: tolerance for the absolute value of prinicipal compent's
19 % eigenvalues. If the eigenvalue is below, the principal component
20 % is discarded. (default='eps')
23 % B: matrix with computed principal component vectors as column vectors
24 % evalues: eigenvalues of the gram matrix, necessary for calculating
28 % 5.10.2005 Bernard Haasdonk
34 % trivial weighting if matrix is not given:
35 if nargin < 3 || isempty(A)
39 if nargin < 4 || isempty(k)
43 if nargin < 5 || isempty(method)
48 epsilon = eps; % = 2e-16
51 % the following is expensive, use of matlab-function instead ?
52 XFixON = orthonormalize(XFix,A,[],method);
57 Xo = X - XFixON * (XFixON' * A * X);
66 %if size(Xo,1)<size(Xo,2) % ordinary correlation matrix
67 % [e,v] = eig(Xo*Xo');
70 % evalues = evalues(end:-1:1);
71 % fi = find(abs(evalues)>epsilon);
74 % via gram matrix of orthogonalized trajectory
76 K = 0.5*(K+K'); % required for rounding problems
78 % generate descending list of eigenvectors/values:
83 % use eigendecomposition:
87 evalues = evalues(end:-1:1);
88 % svd gives similar results:
89 %[ep,vp,dummy] = svd(K);
92 fi = find(abs(evalues)>=epsilon);
93 %fi = 1:length(evalues);
96 % project Xo vectors on eigenvectors
97 B = Xo * ep(:,fi) * diag(evalues(fi).^(-0.5));
102 % ensure that only real valued vectors are returned
104 % disp('complex eigenvector occured: please check!');
109 % the following is necessary for correct scaling and improving the
110 % orthogonality, i.e. e' A e => identity up to 1e-8
111 % B' A B => identity up to 1e-16
113 B = orthonormalize(B,A,epsilon,method);
114 %max(max(abs(eye(size(B,2))-B' * A * B)))
function Xdel = delzerocolumns(X, epsilon, A)
function deleting zero-columns from matrix X.