簡單易學的機器學習算法—譜聚類(Spectal Clustering)

一、復雜網絡中的一些基本概念

1、復雜網絡的表示

為了能夠讓每個類都有合理的大小，目標函數中應該使得足夠大，則提出了或者：

對于一個有n個頂點的圖，其Laplacian矩陣定義為：

1. Laplacian矩陣L的最小特征值是0，相應的特征向量是I；
2. Laplacian矩陣L有n個非負實特征值：，且對于任何一個實向量f，都有下面的式子成立：

定義向量，且

而已知：，則

由于是個常數，故要求的最小值，即求的最小值。則新的目標函數為：

 

假設是Laplacian矩陣L的特征值，F是特征值對應的特征向量，則有：

在上式的兩端同時左乘

已知，則，上式可以轉化為：

要求，即只需求得最小特征值。由Laplacian矩陣的性質可知，Laplacian矩陣的最小特征值為0。由Rayleigh-Ritz理論，可以取第2小特征值。

對于求解出來的特征向量中的每一個分量，根據每個分量的值來判斷對應的點所屬的類別：

[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
#coding:UTF-8  
'''''
Created on 2015年5月12日
 
@author: zhaozhiyong
'''  
from __future__ import division  
import scipy.io as scio  
from scipy import sparse  
from scipy.sparse.linalg.eigen import arpack#這里只能這么做，不然始終找不到函數eigs  
from numpy import *  
 
 
def spectalCluster(data, sigma, num_clusters):  
    print "將鄰接矩陣轉換成相似矩陣"  
    #先完成sigma ！= 0  
    print "Fixed-sigma譜聚類"  
    data = sparse.csc_matrix.multiply(data, data)  
 
    data = -data / (2 * sigma * sigma)  
      
    S = sparse.csc_matrix.expm1(data) + sparse.csc_matrix.multiply(sparse.csc_matrix.sign(data), sparse.csc_matrix.sign(data))     
      
    #轉換成Laplacian矩陣  
    print "將相似矩陣轉換成Laplacian矩陣"  
    D = S.sum(1)#相似矩陣是對稱矩陣  
    D = sqrt(1 / D)  
    n = len(D)  
    D = D.T  
    D = sparse.spdiags(D, 0, n, n)  
    L = D * S * D  
      
    #求特征值和特征向量  
    print "求特征值和特征向量"  
    vals, vecs = arpack.eigs(L, k=num_clusters,tol=0,which="LM")    
      
    # 利用k-Means  
    print "利用K-Means對特征向量聚類"  
    #對vecs做正規化  
    sq_sum = sqrt(multiply(vecs,vecs).sum(1))  
    m_1, m_2 = shape(vecs)  
    for i in xrange(m_1):  
        for j in xrange(m_2):  
            vecs[i,j] = vecs[i,j]/sq_sum[i]  
      
    myCentroids, clustAssing = kMeans(vecs, num_clusters)  
      
    for i in xrange(shape(clustAssing)[0]):  
        print clustAssing[i,0]  
      
 
def randCent(dataSet, k):  
    n = shape(dataSet)[1]  
    centroids = mat(zeros((k,n)))#create centroid mat  
    for j in range(n):#create random cluster centers, within bounds of each dimension  
        minJ = min(dataSet[:,j])   
        rangeJ = float(max(dataSet[:,j]) - minJ)  
        centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))  
    return centroids  
 
def distEclud(vecA, vecB):  
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)  
 
def kMeans(dataSet, k):  
    m = shape(dataSet)[0]  
    clusterAssment = mat(zeros((m,2)))#create mat to assign data points to a centroid, also holds SE of each point  
    centroids = randCent(dataSet, k)  
    clusterChanged = True  
    while clusterChanged:  
        clusterChanged = False  
        for i in range(m):#for each data point assign it to the closest centroid  
            minDist = inf; minIndex = -1  
            for j in range(k):  
                distJI = distEclud(centroids[j,:],dataSet[i,:])  
                if distJI < minDist:  
                    minDist = distJI; minIndex = j  
            if clusterAssment[i,0] != minIndex: clusterChanged = True  
            clusterAssment[i,:] = minIndex,minDist**2  
        #print centroids  
        for cent in range(k):#recalculate centroids  
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster  
            centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean   
    return centroids, clusterAssment  
 
 
if __name__ == '__main__':  
    # 導入數據集  
    matf = 'E://data_sc//corel_50_NN_sym_distance.mat'  
    dataDic = scio.loadmat(matf)  
    data = dataDic['A']  
    # 譜聚類的過程  
    spectalCluster(data, 20, 18)  

2、網上提供的一個Matlab代碼
[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
function [cluster_labels evd_time kmeans_time total_time] = sc(A, sigma, num_clusters)  
%SC Spectral clustering using a sparse similarity matrix (t-nearest-neighbor).  
%  
%   Input  : A              : N-by-N sparse distance matrix, where  
%                             N is the number of data  
%            sigma          : sigma value used in computing similarity,  
%                             if 0, apply self-tunning technique  
%            num_clusters   : number of clusters  
%  
%   Output : cluster_labels : N-by-1 vector containing cluster labels  
%            evd_time       : running time for eigendecomposition  
%            kmeans_time    : running time for k-means  
%            total_time     : total running time  
 
%  
% Convert the sparse distance matrix to a sparse similarity matrix,  
% where S = exp^(-(A^2 / 2*sigma^2)).  
% Note: This step can be ignored if A is sparse similarity matrix.  
%  
disp('Converting distance matrix to similarity matrix...');  
tic;  
n = size(A, 1);  
 
if (sigma == 0) % Selftuning spectral clustering  
  % Find the count of nonzero for each column  
  disp('Selftuning spectral clustering...');  
  col_count = sum(A~=0, 1)';  
  col_sum = sum(A, 1)';  
  col_mean = col_sum ./ col_count;  
  [x y val] = find(A);  
  A = sparse(x, y, -val.*val./col_mean(x)./col_mean(y)./2);  
  clear col_count col_sum col_mean x y val;  
else % Fixed-sigma spectral clustering  
  disp('Fixed-sigma spectral clustering...');  
  A = A.*A;  
  A = -A/(2*sigma*sigma);  
end  
 
% Do exp function sequentially because of memory limitation  
num = 2000;  
num_iter = ceil(n/num);  
S = sparse([]);  
for i = 1:num_iter  
  start_index = 1 + (i-1)*num;  
  end_index = min(i*num, n);  
  S1 = spfun(@exp, A(:,start_index:end_index)); % sparse exponential func  
  S = [S S1];  
  clear S1;  
end  
clear A;  
toc;  
 
%  
% Do laplacian, L = D^(-1/2) * S * D^(-1/2)  
%  
disp('Doing Laplacian...');  
D = sum(S, 2) + (1e-10);  
D = sqrt(1./D); % D^(-1/2)  
D = spdiags(D, 0, n, n);  
L = D * S * D;  
clear D S;  
time1 = toc;  
 
%  
% Do eigendecomposition, if L =  
%   D^(-1/2) * S * D(-1/2)    : set 'LM' (Largest Magnitude), or  
%   I - D^(-1/2) * S * D(-1/2): set 'SM' (Smallest Magnitude).  
%  
disp('Performing eigendecomposition...');  
OPTS.disp = 0;  
[V, val] = eigs(L, num_clusters, 'LM', OPTS);  
time2 = toc;  
 
%  
% Do k-means  
%  
disp('Performing kmeans...');  
% Normalize each row to be of unit length  
sq_sum = sqrt(sum(V.*V, 2)) + 1e-20;  
U = V ./ repmat(sq_sum, 1, num_clusters);  
clear sq_sum V;  
cluster_labels = k_means(U, [], num_clusters);  
total_time = toc;  
 
%  
% Calculate and show time statistics  
%  
evd_time = time2 - time1  
kmeans_time = total_time - time2  
total_time  
disp('Finished!');  

[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
function cluster_labels = k_means(data, centers, num_clusters)  
%K_MEANS Euclidean k-means clustering algorithm.  
%  
%   Input    : data           : N-by-D data matrix, where N is the number of data,  
%                               D is the number of dimensions  
%              centers        : K-by-D matrix, where K is num_clusters, or  
%                               'random', random initialization, or  
%                               [], empty matrix, orthogonal initialization  
%              num_clusters   : Number of clusters  
%  
%   Output   : cluster_labels : N-by-1 vector of cluster assignment  
%  
%   Reference: Dimitrios Zeimpekis, Efstratios Gallopoulos, 2006.  
%              http://scgroup.hpclab.ceid.upatras.gr/scgroup/Projects/TMG/  
 
%  
% Parameter setting  
%  
iter = 0;  
qold = inf;  
threshold = 0.001;  
 
%  
% Check if with initial centers  
%  
if strcmp(centers, 'random')  
  disp('Random initialization...');  
  centers = random_init(data, num_clusters);  
elseif isempty(centers)  
  disp('Orthogonal initialization...');  
  centers = orth_init(data, num_clusters);  
end  
 
%  
% Double type is required for sparse matrix multiply  
%  
data = double(data);  
centers = double(centers);  
 
%  
% Calculate the distance (square) between data and centers  
%  
n = size(data, 1);  
x = sum(data.*data, 2)';  
X = x(ones(num_clusters, 1), :);  
y = sum(centers.*centers, 2);  
Y = y(:, ones(n, 1));  
P = X + Y - 2*centers*data';  
 
%  
% Main program  
%  
while 1  
  iter = iter + 1;  
 
  % Find the closest cluster for each data point  
  [val, ind] = min(P, [], 1);  
  % Sum up data points within each cluster  
  P = sparse(ind, 1:n, 1, num_clusters, n);  
  centers = P*data;  
  % Size of each cluster, for cluster whose size is 0 we keep it empty  
  cluster_size = P*ones(n, 1);  
  % For empty clusters, initialize again  
  zero_cluster = find(cluster_size==0);  
  if length(zero_cluster) > 0  
    disp('Zero centroid. Initialize again...');  
    centers(zero_cluster, :)= random_init(data, length(zero_cluster));  
    cluster_size(zero_cluster) = 1;  
  end  
  % Update centers  
  centers = spdiags(1./cluster_size, 0, num_clusters, num_clusters)*centers;  
 
  % Update distance (square) to new centers  
  y = sum(centers.*centers, 2);  
  Y = y(:, ones(n, 1));  
  P = X + Y - 2*centers*data';  
 
  % Calculate objective function value  
  qnew = sum(sum(sparse(ind, 1:n, 1, size(P, 1), size(P, 2)).*P));  
  mesg = sprintf('Iteration %d:\n\tQold=%g\t\tQnew=%g', iter, full(qold), full(qnew));  
  disp(mesg);  
 
  % Check if objective function value is less than/equal to threshold  
  if threshold >= abs((qnew-qold)/qold)  
    mesg = sprintf('\nkmeans converged!');  
    disp(mesg);  
    break;  
  end  
  qold = qnew;  
end  
 
cluster_labels = ind';  
 
 
%-----------------------------------------------------------------------------  
function init_centers = random_init(data, num_clusters)  
%RANDOM_INIT Initialize centroids choosing num_clusters rows of data at random  
%  
%   Input : data         : N-by-D data matrix, where N is the number of data,  
%                          D is the number of dimensions  
%           num_clusters : Number of clusters  
%  
%   Output: init_centers : K-by-D matrix, where K is num_clusters  
rand('twister', sum(100*clock));  
init_centers = data(ceil(size(data, 1)*rand(1, num_clusters)), :);  
 
function init_centers = orth_init(data, num_clusters)  
%ORTH_INIT Initialize orthogonal centers for k-means clustering algorithm.  
%  
%   Input : data         : N-by-D data matrix, where N is the number of data,  
%                          D is the number of dimensions  
%           num_clusters : Number of clusters  
%  
%   Output: init_centers : K-by-D matrix, where K is num_clusters  
 
%  
% Find the num_clusters centers which are orthogonal to each other  
%  
Uniq = unique(data, 'rows'); % Avoid duplicate centers  
num = size(Uniq, 1);  
first = ceil(rand(1)*num); % Randomly select the first center  
init_centers = zeros(num_clusters, size(data, 2)); % Storage for centers  
init_centers(1, :) = Uniq(first, :);  
Uniq(first, :) = [];  
c = zeros(num-1, 1); % Accumalated orthogonal values to existing centers for non-centers  
% Find the rest num_clusters-1 centers  
for j = 2:num_clusters  
  c = c + abs(Uniq*init_centers(j-1, :)');  
  [minimum, i] = min(c); % Select the most orthogonal one as next center  
  init_centers(j, :) = Uniq(i, :);  
  Uniq(i, :) = [];  
  c(i) = [];  
end  
clear c Uniq;  

個人的一點認識：譜聚類的過程相當于先進行一個非線性的降維，然后在這樣的低維空間中再利用聚類的方法進行聚類。