rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
model_PCA_fixspace.m
1 function B = model_PCA_fixspace(model,model_data,X,XFix,k)
2 %function B = model_PCA_fixspace(model,model_data,X,XFix [, [A], [k] ])
3 %
4 % function computing a PCA basis of the set of column vectors X projected
5 % on the orthogonal complement of span(XFix)
6 % The computation is performed explicitly using the
7 % gram matrix.
8 % Additionally, the sparse matrix of the inner product A on the
9 % vector space can be provided, i.e.
10 % i.e. <x,x'> = x'*A*x
11 % Additionally, the number of principal components can be specified
12 
13 % 5.10.2005 Bernard Haasdonk
14 
15 %epsilon = 1e-10;
16 % epsilon = eps; % = 2e-16
17 
18 % trivial weighting if matrix is not given:
19 
20 A = model_data.W;
21 
22 if nargin < 4
23  k = size(X,2);
24 end;
25 
26 % the following is expensive, use of matlab-function instead ?
27 XFixON = model_orthonormalize(model,model_data,XFix);
28 %XFixON = delzerocolumns(XFixON,[],A);
29 
30 if ~isempty(XFix)
31  Xo = X - XFixON * (XFixON' * A * X);
32 else
33  Xo = X;
34 end;
35 
36 %if size(Xo,1)<size(Xo,2) % ordinary correlation matrix
37 % [e,v] = eig(Xo*Xo');
38 % e = e(:,end:-1:1);
39 % evalues = diag(v);
40 % evalues = evalues(end:-1:1);
41 % fi = find(abs(evalues)>epsilon);
42 % e = e(:,fi);
43 %else
44 % via gram matrix of orthogonalized trajectory
45 K = Xo'*A*Xo;
46 K = 0.5*(K+K'); % required for rounding problems
47 
48 % generate descending list of eigenvectors/values:
49 if k<size(K,2)
50  [ep,vp] = eigs(K,k);
51  evalues = diag(vp);
52 else
53  [ep,vp] = eig(K);
54  ep = ep(:,end:-1:1);
55  evalues = diag(vp);
56  evalues = evalues(end:-1:1);
57 end;
58 %fi = find(abs(evalues)>=epsilon);
59 fi = 1:length(evalues);
60 
61 % project Xo vectors on eigenvectors
62 B = Xo * ep(:,fi) * diag(evalues(fi).^(-0.5));
63 %end;
64 
65 % ensure that only real valued vectors are returned
66 while (~isreal(B))
67 % disp('complex eigenvector occured: please check!');
68 % keyboard;
69  B= B(:,1:end-1);
70 end;
71 
72 % the following is theoretically superfluous, but does improve the
73 % orthogonality, i.e. e' A e => identity up to 1e-8
74 % B' A B => identity up to 1e-16
75 %B = orthonormalize(B,A);
76 
77 %| \docupdate
function Xdel = delzerocolumns(X, epsilon, A)
function deleting zero-columns from matrix X.