PCA_fixspace.m

function B = PCA_fixspace(X,XFix,A,k,method,epsilon)
%function B = PCA_fixspace(X,XFix [, [A], [k], [method],[epsilon] ])
%
% function computing a PCA basis of the set of column vectors X projected on
% the orthogonal complement of span(XFix) The computation is performed
% explicitly using the gram matrix.
%
% Parameters:
%   X:       the matrix of column vectors on which the PCA shall be applied
%   XFix:    matrix of column vectors spanning a fix space to which the vectors
%            in 'X' shall be orthonormalized before the PCA.
%   A:       inner product matrix defining the scalar product `\langle x, y
%            \rangle = y^t A x`. (default='1')
%   k:       maximum number of principal components to be computed. (default=
%            number of columns in 'X').
%   method:  string specifying the orthonormalization method, either 'qr' or
%            'gram-schmidt'. (default='gram-schmidt')
%   epsilon: tolerance for the absolute value of prinicipal compent's
%            eigenvalues. If the eigenvalue is below, the principal component
%            is discarded. (default='eps')
%
% Return values:
%   B:       matrix with computed principal component vectors as column vectors


% 5.10.2005 Bernard Haasdonk

%keyboard;

%epsilon = 1e-10;

% trivial weighting if matrix is not given:
if nargin < 3 || isempty(A)
  A = 1;
end;

if nargin < 4 || isempty(k)
  k = size(X,2);
end;

if nargin < 5 || isempty(method)
  method = 'qr';
end;

if nargin < 6
  epsilon = eps; % = 2e-16
end;

% the following is expensive, use of matlab-function instead ?
XFixON = orthonormalize(XFix,A,[],method);
%keyboard;
%XFixON = delzerocolumns(XFixON,[],A);

if ~isempty(XFix)
     Xo = X - XFixON * (XFixON' * A * X); 
else
    Xo = X;
end;

clear('XFix');
clear('XFixON');
clear('X');

%if size(Xo,1)<size(Xo,2) % ordinary correlation matrix
%  [e,v] = eig(Xo*Xo');
%  e = e(:,end:-1:1);
%  evalues = diag(v);
%  evalues = evalues(end:-1:1);
%  fi = find(abs(evalues)>epsilon);
%  e = e(:,fi);
%else  
% via gram matrix of orthogonalized trajectory
K = Xo'*A*Xo;
K = 0.5*(K+K'); % required for rounding problems

% generate descending list of eigenvectors/values:
if k<size(K,2)
  [ep,vp] = eigs(K,k);
  evalues = diag(vp);
else
  % use eigendecomposition:
  [ep,vp] = eig(K);
  ep = ep(:,end:-1:1);
  evalues = diag(vp);
  evalues = evalues(end:-1:1);
  % svd gives similar results:
  %[ep,vp,dummy] = svd(K);
  %evalues = diag(vp);
end;
fi = find(abs(evalues)>=epsilon);
%fi = 1:length(evalues);


% project Xo vectors on eigenvectors
B = Xo * ep(:,fi) * diag(evalues(fi).^(-0.5));
%end;

clear('Xo');

% ensure that only real valued vectors are returned
while (~isreal(B))
%  disp('complex eigenvector occured: please check!');
%  keyboard;
  B= B(:,1:end-1);
end;

% the following is necessary for correct scaling and improving the 
% orthogonality, i.e. e' A e => identity up to 1e-8
%                     B' A B => identity up to 1e-16
%(XFix,A,[],method);
B = orthonormalize(B,A,epsilon,method);
%max(max(abs(eye(size(B,2))-B' * A * B)))

%keyboard;