function [eigval, eigvec, scores, omega] = kpca(Xtrain, kernel_type, kernel_pars ,Xt,etype,n,centr)% Kernel Principal Component Analysis (KPCA)% % >> [eigval, eigvec] = kpca(X, kernel_fct, sig2)% >> [eigval, eigvec, scores] = kpca(X, kernel_fct, sig2, Xt)% % Compute the nb largest eigenvalues and the corresponding rescaled% eigenvectors corresponding with the principal components in the% feature space of the centered kernel matrix. To calculate the% eigenvalue decomposition of this N x N matrix, Matlab's% eig is called by default. The decomposition can also be% approximated by Matlab ('eigs') or by Nystr鰉's method ('eign')% using nb components. In some cases one wants to disable% ('original') the rescaling of the principal components in feature% space to unit length. % % The scores of a test set Xt on the principal components is computed by the call:% % >> [eigval, eigvec, scores] = kpca(X, kernel_fct, sig2, Xt)% % Full syntax% % >> [eigval, eigvec, empty, omega] = kpca(X, kernel_fct, sig2) % >> [eigval, eigvec, empty, omega] = kpca(X, kernel_fct, sig2, [],etype, nb) % >> [eigval, eigvec, empty, omega] = kpca(X, kernel_fct, sig2, [],etype, nb, rescaling) % >> [eigval, eigvec, scores, omega] = kpca(X, kernel_fct, sig2, Xt) % >> [eigval, eigvec, scores, omega] = kpca(X, kernel_fct, sig2, Xt,etype, nb) % >> [eigval, eigvec, scores, omega] = kpca(X, kernel_fct, sig2, Xt,etype, nb, rescaling)% % Outputs    %   eigval       : N (nb) x 1 vector with eigenvalues values%   eigvec       : N x N (N x nb) matrix with the principal directions%   Xt(*)        : Nt x nb matrix with the scores of the test data (or [])%   Omega(*)     : N x N centered kernel matrix% Inputs    %   X            : N x d matrix with the inputs of the training data%   kernel       : Kernel type (e.g. 'RBF_kernel')%   sig2         : Kernel parameter(s) (for linear kernel, use [])%   Xt(*)        : Nt x d matrix with the inputs of the test data (or [])%   etype(*)     : 'svd', 'eig'(*),'eigs','eign'%   nb(*)        : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation%   rescaling(*) : 'original size' ('o') or 'rescaling'(*) ('r')% % See also:%   bay_lssvm, bay_optimize, eign% Copyright (c) 2002,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.ac.be/sista/lssvmlab%% defaults %nb_data = size(Xtrain,1);eval('n=min(n,nb_data);','n=min(10,nb_data);')eval('centr;','centr=''rescaled'';');eval('etype;','etype=''eig'';');eval('Xt;','Xt=[];');%% tests%if ~isempty(Xt) & size(Xt,2)~=size(Xtrain,2),  error('Training points and test points need to have the same dimension');endif ~(strcmpi(etype,'svd') | strcmpi(etype,'eig') | strcmpi(etype,'eigs') | strcmpi(etype,'eign')),  error('Eigenvalue decomposition via ''svd'', ''eig'', ''eigs'' or ''eign''...');endif (strcmpi(etype,'svd') | strcmpi(etype,'eig') | strcmpi(etype,'eigs')),    omega = kernel_matrix(Xtrain,kernel_type, kernel_pars);    % centered kernel matrix : Z*omega*Z  %Zc = eye(nb_data) - ones(nb_data)./nb_data;  %omega = Zc*omega*Zc;  Meanvec = mean(omega);  MM = mean(Meanvec);  for i=1:nb_data,    omega(i,:) = omega(i,:)-Meanvec(i);  end  for i=1:nb_data,    omega(:,i) = omega(:,i)-Meanvec(i);  end  omega = omega+MM;      %  % eigenvalues are computed with more stable svd  %    % numerical stability issues  omega = (omega+omega')./2;    if strcmpi(etype,'svd'),    [eigvec, eigval,ff] = svd(omega); clear ff  elseif strcmpi(etype,'eig'),    [eigvec, eigval] = eig(omega);   elseif (strcmpi(etype,'eigs')),    [eigvec, eigval] = eigs(omega,n);  end  eigval = diag(eigval)./(nb_data-1);elseif strcmpi(etype,'eign'),  if nargout>1,    [eigvec,eigval] = eign(Xtrain,kernel_type,kernel_pars, n);   else    eigval = eign(Xtrain,kernel_type,kernel_pars, n);   end  omega = [];  eigval = (eigval)./(nb_data-1);  Meanvec = [];  MM = [];else  error('Unknown type for eigenvalue approximation');end%% only keep relevant eigvals & eigvec%peff = find(eigval>1000*eps);%eigval = eigval(peff);neff = length(peff);%if nargout>1, eigvec = eigvec(:,peff); end% rescaling the eigenvectorsif (centr(1) =='r' & nargout>1),  %disp('rescaling the eigvec');  for i=1:neff,    eigvec(:,i) = eigvec(:,i)./sqrt(eigvec(:,i)'*eigval(i)*eigvec(:,i));  endend%% compute scores%if ~isempty(Xt),  omega_t = kernel_matrix(Xtrain,kernel_type, kernel_pars,Xt);  if isempty(Meanvec),    if size(omega_t,2)>10,  Meanvec = mean(omega_t,2); MM = mean(Meanvec); end  end  for i=1:size(omega_t,1),    omega_t(i,:) = omega_t(i,:)-Meanvec(i);  end  for i=1:size(omega_t,2),    omega_t(:,i) = omega_t(:,i)-Meanvec(i);  end  omega_t = omega_t+MM;      scores = omega_t'*eigvec;else  scores = [];end