kpca.m

KPCA的程序,比较好用,有详细的注释,希望对大家有用!

function [s,UU,salid]=kpca(x,ker,par,x2)
%
%  This function solves the Kernal Principal Component Analysis (KPCA).
%
%     [s,UU,Ncom,salid]=kpca(x,ker,par,x2)
%
% This function solves the SVM using the IRWLS procedure. 
% Parameters:
%
%              x: matrix of vectors. Each row of x1 represents a d-dimensional vector. 
%                 The number of rows is the number of training samples.
%              ker: the kernel to be employed. It has to be a string out of: 'linear', 'poly_h', 'poly_i' and 'rbf'.
%                 Type help kernel for further information.
%              par: is the parameter of the used kernel. type help kernel for further information.
%              x2: matrix of vectors, whihc one wants to extract the components from using the components of x. 
%                 Each row of x1 represents a d-dimensional vector. 
%                 The number of rows is the number of training samples.
%
% Outputs:
%
%              s: The eigenvalues.
%              UU: The eigenvectors.
%              salid: The projection on the eigenvector os the input x2
%              vectors in the feature space. If x2 is measing it gives the
%              projection of x1.
%              rho: the margin.
%              T: running time.
%              Lp: A vector that show the value of the primal \|w\|^2+C\sum_{i=0}^n \xi_i in each iteration
%               of the algorithm.
%
% Example:
%
%       with linear kernel:
%               x1=.3*randn(2,50);
%               x2=[.3*randn(1,50)+2;.3*randn(1,50)];
%               x3=[.3*randn(1,50)+1;.3*randn(1,50)+2];
%               x=[x1';x2';x3'];
% 
%               v=-1:.1:3;
%               for i=1:length(v)
%                   X(i:length(v):length(v)^2,1)=v';
%                   X(1+length(v)*(i-1):length(v)*i,2)=v';
%               end
% 
%               [s,UU,salid]=kpca(x,'rbf',.5,X);
% 
%               M1=zeros(length(v));
%               for i=1:length(v)
%                   M1(i,:)=salid(1+(i-1)*length(v):length(v)*i,1)';
%               end
%               figure(1)              
%               pcolor(v,v,M1');
%               shading interp
%               hold on
%               plot(x(:,1),x(:,2),'k.','markersize',20)
%               contour(v,v,M1',[-1.1:.2:1.1],'k');
%               hold off
%               

% Author:   Fernando Perez-Cruz ([e-mail])
% Version:  1.0
% Date:     20th January 2004.

[n,d]=size(x);
K=kernel(ker,x,x,par);

uno=ones(n)/n;

K2=(K-uno*K-K*uno+uno*K*uno);

[U,S,V]=svd(K2);
s=diag(S);

Ncom=length(find(s>10^-10));
UU=zeros(n,Ncom);
for i=1:Ncom
    UU(:,i)=U(:,i)/sqrt(s(i));
end
s=s(1:Ncom);

if(nargin==4)
    salid=kernel(ker,x2,x,par)*UU;
else
    salid=K*UU;
end