mfbox_pkmeans.m

toolbox for spm 5 for data, model free analysis

function [centroids,clusters,err,qerrs]=mfbox_pkmeans(method,D,n,epochs,verbose)

% k-means clustering method in projective space 
%
% [centroids,clusters,err,qerrs] = mfbox_pkmeans(method,D,n[,epochs[,verbose]])
%
%  Input and output arguments ([]'s are optional):
%   method    (string) 'batch' or 'seq'
%   D         (matrix) size dlen x dim, the data
%             (struct) map or data struct
%   n         (scalar) number of centroids
%             (matrix) size n x dim, initial values for the centroids
%   [epochs]  (scalar) number of training epochs, 100 by default
%   [verbose] (scalar) verbose level, 0 by default
%
%   centroids (matrix) size n x dim, the k-means centroids
%   clusters  (vector) size dlen x 1, cluster number for each sample
%   err       (scalar) total quantization error for the data set
%   qerr      (vector) quantization error in the clusters
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% input arguments and initialization

data = pnorm(D); 

[dlen,dim] = size(data);

if (prod(size(n))==1)
  temp = randperm(dlen);
  centroids = data(temp(1:n),:);
else
  centroids = n;
  n = size(centroids,1);
end

if (nargin<4) epochs = 100; end
if (nargin<5) verbose = 0; end

lr = 0.5;                      % initial learning rate for sequential k-means
clusters  = zeros(1,dlen);       

switch method
 case 'seq' % sequential k-means
  len = epochs * dlen;
  l_rate = linspace(lr,0,len);
  order  = randperm(dlen);
  for iter=1:len
    x  = data(order(rem(iter,dlen)+1),:);
    dx = pdist(x,centroids);
    [dist,nearest] = min(dx);
    if (x*centroids(nearest,:)'<0) x = -x; end 
    o = pdist(x,centroids(nearest,:));
    if (o~=0) x = pnorm(x - (o*centroids(nearest,:))); end
    centroids(nearest,:) = pnorm(centroids(nearest,:) + l_rate(iter)*x);
  end 
 case 'batch' % batch k-means
  iter = 0;
  old_clusters = zeros(1,dlen);
  while (iter<epochs)
    [dummy,clusters] = min(pdist(centroids,data));
    for i=1:n,
      f = find(clusters==i);
      s = length(f);
      if (s) centroids(i,:) = psum(data(f,:)); end
    end    
    if ((iter>0)&&all(old_clusters==clusters))
      if (verbose) fprintf(1,'Convergence in %d iterations\n', iter); end
      break; 
    end    
    old_clusters = clusters;
    iter = iter+1;
  end  
end

% final partitioning and error
[qerrs,clusters] = min(pdist(centroids,data));
err = sum(qerrs);

if (size(clusters,1)==1) clusters = clusters'; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function d=pdist(x,y)
%
d = 1-(x*y').^2;
return

function x=pnorm(x)
%
n = sqrt(sum(x.^2,2));
x = x./(n*ones(1,size(x,2)));
return

function s=psum(x)
%
n = size(x,2);
[E,D] = eig(x'*x);

[dummy,q] = sort(diag(D));
s = E(:,q(n))';
return