nt_pca_kmeans

[topcs,pwr]=nt_pca_kmeans(x,nkeep) - PCA preceded by kmeans for speed

  topcs: PCA matrix

  x: data (time*channels or time*channels*trials)
  nkeep: desired number of PCs

 The kmeans implementation of VLFeat appears to be very efficient, 
 so we use it to find a small set of clusters to which we apply PCA.

 This differs from the usual way of combining PCA and kmeans, which is
 to use PCA to initialize kmeans.

0001 function [topcs,pwr]=nt_pca_kmeans(x,shifts,nkeep)
0002 %[topcs,pwr]=nt_pca_kmeans(x,nkeep) - PCA preceded by kmeans for speed
0003 %
0004 %  topcs: PCA matrix
0005 %
0006 %  x: data (time*channels or time*channels*trials)
0007 %  nkeep: desired number of PCs
0008 %
0009 % The kmeans implementation of VLFeat appears to be very efficient,
0010 % so we use it to find a small set of clusters to which we apply PCA.
0011 %
0012 % This differs from the usual way of combining PCA and kmeans, which is
0013 % to use PCA to initialize kmeans.
0014 
0015 if nargin<3; error('!'); end
0016 if isempty(shifts); shifts=0; end
0017 if shifts
0018     x=nt_multishift(x,shifts);
0019 end
0020 if nkeep>size(x,2); error('!'); end
0021 
0022 x=nt_unfold(x);
0023 [nsamples,nchans]=size(x);
0024 
0025 if nkeep>nsamples; error('nkeep greater than number of samples'); end
0026 if nkeep>nchans; error('nkeep greater than number of channels'); end
0027 
0028 
0029 % normalize channels
0030 nrm=sqrt(mean(x.^2));
0031 tonrm=diag(1./nrm); % normalization matrix
0032 tonrm(find(isnan(x)))=0;
0033 x=x*tonrm;
0034 
0035 % cluster
0036 nclusters=round(sqrt(nchans)); 
0037 [C,A]=vl_kmeans(x,nclusters,'algorithm','elkan');
0038 
0039 % perform PCA on each cluster
0040 topcs1=zeros(nchans);
0041 pwr=zeros(1,nchans);
0042 idx=0;
0043 for k=1:nclusters
0044     %[nclusters k numel(find(A==k))]
0045     [m,p]=nt_pca0(x(:,find(A==k)));
0046     topcs1(find(A==k),idx+(1:size(m,2))) = m;
0047     pwr(idx+(1:size(m,2)))=p;
0048     idx=idx+size(m,2);
0049 end
0050 
0051 % sort by decreasing power
0052 [~,idx]=sort(pwr,'descend');
0053 
0054 topcs2=nt_pca0(nt_mmat(x,topcs1(:,idx(1:nkeep))));
0055 
0056 topcs=tonrm*(topcs1(:,idx(1:nkeep))*topcs2);
0057 
0058 
0059 
0060 
0061 
0062 
0063 
0064