kpca_calc.m

KPCA是一种非线性的盲源分离方法

function basis = kpca_calc(xs,kernel,d,kmataxis);
%KPCA_CALC calculates a kernel PCA basis.
%
%   usage
%      basis = kpca_calc(xs,kernel,d);
%
%   input
%      xs       matrix of column vectors
%      kernel   a chosen kernel, default = {'gaussian',1}
%      d        number of eigenvectors (give for efficiency),
%               default = size(xs,2)
%      kmataxis is a figure handle where the kernel matrix will be
%               plotted (default = 0 no plot)
%
%   output
%      basis    struct containing the following entries
%                   basis.V        eigenvectors
%                   basis.Lambda   eigenvalues
%                   basis.xs       used vectors
%                   basis.kernel   used kernel
%
%   see also
%      kpca_plot, kpca_map
%
%   STH * 12MAR2002

if ~exist('kernel')|isempty(kernel), kernel = {'gaussian',1}; end
if ~exist('d')|isempty(d), d = size(xs,2); end
if ~exist('kmataxis')|isempty(kmataxis), kmataxis = 0; end

% d can't be larger than the number of samples
if d>size(xs,2)
  warning('d is larger than the number of samples, resetting d')
  d = size(xs,2);
end

xsc = size(xs,2);   % column of xs

% calculate the kernel matrix
K = kpca_matrix(xs,xs,kernel);
if kmataxis>0
  cf = gcf;
  figure(kmataxis)
  imagesc(K)
  figure(cf)
end

% center the kernel matrix
sk = size(K,1);             % note, K is square matrix
rowK = sum(K)/sk;           % the sums of the columns
allK = sum(K(:))/(sk*sk);   % the sum of all entries
K = K - repmat(rowK,[sk 1]) - repmat(rowK',[1 sk]) + repmat(allK,[sk sk]);

% find the eigenvectors and eigenvalues
switch 2
 case 1
  [V,Lambda] = jdqr(K/sk,d);
 case 2
  opts.disp = 0;
  [V,Lambda,flag] = eigs(K/sk,d,'LM',opts);
  if flag
    warning([mfilename ': not all eigenvalues converged'])
  end
end

% we can not assume that the eigenvalues are sorted
[dummy, ind] = sort(-diag(Lambda));
Lambda = Lambda(ind,ind);
V = V(:,ind);

% due to numerical instabilities some eigenvalues might be negative
% or smaller than eps, we want to ignore those
valid = find(diag(Lambda)<2*eps);
if length(valid)<1
  % all eigenvalues are valid, keep d unchanged
else
  % some are not valid
  d = valid(1)-1;
  warning([mfilename ': some eigenvalues of kernel matrix are less than eps'])
end
clear valid

% cut off those eigenvalues and eigenvectors
V = V(:,1:d);
Lambda = Lambda(1:d,1:d);

% normalize the eigenvectors in feature space
V = V*inv(sqrtm(sk*Lambda));

% assign struct
basis.V = V;
basis.Lambda = Lambda;
basis.xs = xs;
basis.kernel = kernel;