klfda.m

这是一个关于Fisher线性判别分析的Matlab的m文件

KLFDA/demo_KLFDA.m                                                                                  100666  000000  000000         4647 11165617660  11411  0                                                                                                    ustar 00sugi                                                                                                                                                                                                                                                   % test_KLFDA.m
%
% (c) Masashi Sugiyama, Department of Compter Science, Tokyo Institute of Technology, Japan.
%     [e-mail],     http://sugiyama-www.cs.titech.ac.jp/~sugi/software/LFDA/

clear all;

rand('state',0);
randn('state',0);

%%%%%%%%%%%%%%%%%%%%%% Generating data
n1a=100;
n1b=100;
n2=100;
X1a=[randn(2,n1a).*repmat([1;2],[1 n1a])+repmat([-6;0],[1 n1a])];
X1b=[randn(2,n1b).*repmat([1;2],[1 n1b])+repmat([ 6;0],[1 n1b])];
X2= [randn(2,n2 ).*repmat([1;2],[1 n2 ])+repmat([ 0;0],[1 n2 ])];
X=[X1a X1b X2];
Y=[ones(n1a+n1b,1);2*ones(n2,1)];

%%%%%%%%%%%%%%%%%%%%%% Computing LFDA solution

%Gaussian kernel
K1=Kmatrix_Gauss(X,1);
[T1,Z1]=KLFDA(K1,Y,1);

% linear kernel
K2=X'*X;
[T2,Z2]=KLFDA(K2,Y,1);

%%%%%%%%%%%%%%%%%%%%%% Displaying original 2D data
figure(1)
clf
hold on

set(gca,'FontName','Helvetica')
set(gca,'FontSize',12)
h=plot(X(1,Y==1),X(2,Y==1),'bo');
h=plot(X(1,Y==2),X(2,Y==2),'rx');
axis equal
axis([-10 10 -10 10])
title('Original 2D data')

set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 12 12]);
print -depsc original_data

%%%%%%%%%%%%%%%%%%%%%% Displaying projected data (Gauss kernel)
figure(2)
clf

subplot(2,1,1)
hold on
set(gca,'FontName','Helvetica')
set(gca,'FontSize',12)
hist(Z1(Y==1),linspace(min(Z1),max(Z1),30));
h=get(gca,'Children');
set(h,'FaceColor','b')
axis([min(Z1) max(Z1) 0 inf])
title('1D projection by KLFDA (Gauss kernel)')

subplot(2,1,2)
hold on
set(gca,'FontName','Helvetica')
set(gca,'FontSize',12)
hist(Z1(Y==2),linspace(min(Z1),max(Z1),30));
h=get(gca,'Children');
set(h,'FaceColor','r')
axis([min(Z1) max(Z1) 0 inf])

set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 12 12]);
print -depsc projected_data_KLFDA_Gauss

%%%%%%%%%%%%%%%%%%%%%% Displaying projected data (linear kernel)
figure(3)
clf

subplot(2,1,1)
hold on
set(gca,'FontName','Helvetica')
set(gca,'FontSize',12)
hist(Z2(Y==1),linspace(min(Z2),max(Z2),30));
h=get(gca,'Children');
set(h,'FaceColor','b')
axis([min(Z2) max(Z2) 0 inf])
title('1D projection by KLFDA (linear kernel)')

subplot(2,1,2)
hold on
set(gca,'FontName','Helvetica')
set(gca,'FontSize',12)
hist(Z2(Y==2),linspace(min(Z2),max(Z2),30));
h=get(gca,'Children');
set(h,'FaceColor','r')
axis([min(Z2) max(Z2) 0 inf])

set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 12 12]);
print -depsc projected_data_KLFDA_linear

                                                                                         KLFDA/KLFDA.m                                                                                       100666  000000  000000         4613 11165656311  10372  0                                                                                                    ustar 00sugi                                                                                                                                                                                                                                                   function [T,Z]=KLFDA(K,Y,r,metric,kNN,reg)
%
% Kernel Local Fisher Discriminant Analysis for Supervised Dimensionality Reduction
%
% Usage:
%       [T,Z]=KLFDA(K,Y,r,metric,kNN,reg)
%
% Input:
%    K:      n x n kernel matrix
%            n --- the number of samples 
%    Y:      n dimensional vector of class labels
%    r:      dimensionality of reduced space (default: d)
%    metric: type of metric in the embedding space (default: 'weighted')
%            'weighted'        --- weighted eigenvectors 
%            'orthonormalized' --- orthonormalized
%            'plain'           --- raw eigenvectors
%    kNN:    parameter used in local scaling method (default: 7)
%    reg:    regularization parameter (default: 0.001)
%
% Output:
%    T: d x r transformation matrix (Z=T'*X)
%    Z: r x n matrix of dimensionality reduced samples 
%
% (c) Masashi Sugiyama, Department of Compter Science, Tokyo Institute of Technology, Japan.
%     [e-mail],     http://sugiyama-www.cs.titech.ac.jp/~sugi/software/LFDA/

if nargin<2
  error('Not enough input arguments.')
end

if nargin<3
  r=1;
end

if nargin<4
  metric='weighted';
end

if nargin<5
  kNN=7;
end

if nargin<6
  reg=0.0001;
end

[n ndum]=size(K);
opts.disp = 0; 
tSb=zeros(n,n);
tSw=zeros(n,n);

for c=unique(Y')
  Kcc=K(Y==c,Y==c);
  Kc=K(:,Y==c);
  nc=size(Kcc,1);

  % Define classwise affinity matrix
  Kccdiag=diag(Kcc);
  distance2=repmat(Kccdiag,1,nc)+repmat(Kccdiag',nc,1)-2*Kcc;
  [sorted,index]=sort(distance2);
  kNNdist2=sorted(kNN+1,:);
  sigma=sqrt(kNNdist2);
  localscale=sigma'*sigma;
  flag=(localscale~=0);
  A=zeros(nc,nc);
  A(flag)=exp(-distance2(flag)./localscale(flag));

  Kc1=sum(Kc,2);
  Z=Kc*(repmat(sum(A,2),[1 n]).*Kc')-Kc*A*Kc';
  tSb=tSb+Z/n+Kc*Kc'*(1-nc/n)+Kc1*Kc1'/n;
  tSw=tSw+Z/nc;
end

K1=sum(K,2);
tSb=tSb-K1*K1'/n-tSw;

tSb=(tSb+tSb')/2;
tSw=(tSw+tSw')/2;

%[eigvec,eigval_matrix]=eig(tSb,tSw+reg*eye(size(tSw)));
opts.disp = 0; 
[eigvec,eigval_matrix]=eigs(tSb,tSw+reg*eye(size(tSw)),r,'la',opts);

eigval=diag(eigval_matrix);
[sort_eigval,sort_eigval_index]=sort(eigval);
T0=eigvec(:,sort_eigval_index(end:-1:1));

switch metric %determine the metric in the embedding space
  case 'weighted'
   T=T0.*repmat(sqrt(sort_eigval(end:-1:1))',[n,1]);
  case 'orthonormalized'
   [T,dummy]=qr(T0,0);
  case 'plain'
   T=T0;
end

Z=T'*K;

                                                                                                                     KLFDA/Kmatrix_Gauss.m                                                                               100666  000000  000000          307 10542644713  12307  0                                                                                                    ustar 00sugi                                                                                                                                                                                                                                                   function K=Kmatrix_Gauss(X,sigma)

  if nargin<=1
    sigma=1;
  end

  [d,n]=size(X);
  X2=sum(X.^2,1);
  distance2=repmat(X2,n,1)+repmat(X2',1,n)-2*X'*X;
  K=exp(-distance2/(2*sigma^2));