简单易学的机器学习算法

一、复杂网络中的一些基本概念

1、复杂网络的表示

    在复杂网络的表示中，复杂网络可以建模成一个图，其中，表示网络中的节点的集合，表示的是连接的集合。在复杂网络中，复杂网络可以是无向图、有向图、加权图或者超图。

2、网络簇结构

    网络簇结构(network cluster structure)也称为网络社团结构(network community structure)，是复杂网络中最普遍和最重要的拓扑属性之一。网络簇是整个网络中的稠密连接分支，具有同簇内部节点之间相互连接密集，不同簇的节点之间相互连接稀疏的特征。

3、复杂网络的分类

    复杂网络主要分为：随机网络，小世界网络和无标度网络。

二、谱方法介绍

1、谱方法的思想

    在复杂网络的网络簇结构存在着同簇节点之间连接密集，不同簇节点之间连接稀疏的特征，是否可以根据这样的特征对网络中的节点进行聚类，使得同类节点之间的连接密集，不同类别节点之间的连接稀疏？

    在谱聚类中定义了“截”函数的概念，当一个网络被划分成为两个子网络时，“截”即指子网间的连接密度。谱聚类的目的就是要找到一种合理的分割，使得分割后形成若干子图，连接不同的子图的边的权重尽可能低，即“截”最小，同子图内的边的权重尽可能高。

2、“截”函数的具体表现形式

    “截”表示的是子网间的密度，即边比较少。以二分为例，将图聚类成两个类：类和类。假设用来表示图的划分，我们需要的结果为：

其中表示的是类别和之间的权重。对于个不同的类别，优化的目标为：

3、基本“截”函数的弊端

    对于上述的“截”函数，最终会导致不好的分割，如二分类问题:

上述的“截”函数通常会将图分割成一个点和其余个点。

4、其他的“截”函数的表现形式

    为了能够让每个类都有合理的大小，目标函数中应该使得足够大，则提出了或者：

其中表示类中包含的顶点的数目

三、Laplacian矩阵

1、Laplacian矩阵的定义

    拉普拉斯矩阵(Laplacian Matrix)，也称为基尔霍夫矩阵，是图的一种矩阵表示形式。

    对于一个有个顶点的图，其Laplacian矩阵定义为：

其中，为图的度矩阵，为图的邻接矩阵。

2、度矩阵的定义

   度矩阵是一个对角矩阵，主角线上的值由对应的顶点的度组成。

   对于一个有个顶点的图，其邻接矩阵为：

其度矩阵为：

其中。

3、Laplacian矩阵的性质

1. Laplacian矩阵是对称半正定矩阵；
2. Laplacian矩阵的最小特征值是，相应的特征向量是；
3. Laplacian矩阵有个非负实特征值：，且对于任何一个实向量，都有下面的式子成立：

性质3的证明：

4、不同的Laplacian矩阵

    除了上述的拉普拉斯矩阵，还有规范化的Laplacian矩阵形式：

四、Laplacian矩阵与谱聚类中的优化函数的关系

1、由Laplacian矩阵到“截”函数

    对于二个类别的聚类问题，优化的目标函数为：

定义向量，且

而已知：，则

而

而

其中，表示的是顶点的数目，对于确定的图来说是个常数。由上述的推导可知，由推导出了，由此可知：Laplacian矩阵与有优化的目标函数之间存在密切的联系。

2、新的目标函数

    由上式可得：

由于是个常数，故要求的最小值，即求的最小值。则新的目标函数为：

 

其中

3、转化到Laplacian矩阵的求解

    假设是Laplacian矩阵的特征值，是特征值对应的特征向量，则有：

在上式的两端同时左乘

已知，则，上式可以转化为：

要求，即只需求得最小特征值。由Laplacian矩阵的性质可知，Laplacian矩阵的最小特征值为。由Rayleigh-Ritz理论，可以取第2小特征值。

五、从二类别聚类到多类别聚类

1、二类别聚类

    对于求解出来的特征向量中的每一个分量，根据每个分量的值来判断对应的点所属的类别：

2、多类别聚类

    对于求出来的前个特征向量，可以利用K-Means聚类方法对其进行聚类，若前个特征向量为，这样便由特征向量构成如下的特征向量矩阵：

将特征向量矩阵中的每一行最为一个样本，利用K-Means聚类方法对其进行聚类。

六、谱聚类的过程

1、基本的结构

    基于以上的分析，谱聚类的基本过程为：

1. 对于给定的图，求图的度矩阵和邻接矩阵；
2. 计算图的Laplacian矩阵；
3. 对Laplacian矩阵进行特征值分解，取其前个特征值对应的特征向量，构成的特征向量矩阵；
4. 利用K-Means聚类算法对上述的的特征向量矩阵进行聚类，每一行代表一个样本点。

2、利用相似度矩阵的构造方法

   上述的方法是通过图的度矩阵和邻接矩阵来构造Laplacian矩阵，也可以通过相似度矩阵的方法构造Laplacian矩阵，其方法如下：

   相似度矩阵是由权值矩阵得到：

其中

再利用相似度矩阵构造Laplacian矩阵：

其中为相似度矩阵的度矩阵。

注意：在第一种方法中，求解的是Laplacian矩阵的前个最小特征值对应的特征向量，在第二种方法中，求解的是Laplacian矩阵的前个最大特征值对应的特征向量

七、实验代码

1、自己实现的一个

[python] view plain copy

1. #coding:UTF-8  
2. ''''' 
3. Created on 2015年5月12日 
4.  
5. @author: zhaozhiyong 
6. '''  
7. from __future__ import division  
8. import scipy.io as scio  
9. from scipy import sparse  
10. from scipy.sparse.linalg.eigen import arpack#这里只能这么做，不然始终找不到函数eigs  
11. from numpy import *  
12.   
13.   
14. def spectalCluster(data, sigma, num_clusters):  
15.     print "将邻接矩阵转换成相似矩阵"  
16.     #先完成sigma ！= 0  
17.     print "Fixed-sigma谱聚类"  
18.     data = sparse.csc_matrix.multiply(data, data)  
19.   
20.     data = -data / (2 * sigma * sigma)  
21.       
22.     S = sparse.csc_matrix.expm1(data) + sparse.csc_matrix.multiply(sparse.csc_matrix.sign(data), sparse.csc_matrix.sign(data))     
23.       
24.     #转换成Laplacian矩阵  
25.     print "将相似矩阵转换成Laplacian矩阵"  
26.     D = S.sum(1)#相似矩阵是对称矩阵  
27.     D = sqrt(1 / D)  
28.     n = len(D)  
29.     D = D.T  
30.     D = sparse.spdiags(D, 0, n, n)  
31.     L = D * S * D  
32.       
33.     #求特征值和特征向量  
34.     print "求特征值和特征向量"  
35.     vals, vecs = arpack.eigs(L, k=num_clusters,tol=0,which="LM")    
36.       
37.     # 利用k-Means  
38.     print "利用K-Means对特征向量聚类"  
39.     #对vecs做正规化  
40.     sq_sum = sqrt(multiply(vecs,vecs).sum(1))  
41.     m_1, m_2 = shape(vecs)  
42.     for i in xrange(m_1):  
43.         for j in xrange(m_2):  
44.             vecs[i,j] = vecs[i,j]/sq_sum[i]  
45.       
46.     myCentroids, clustAssing = kMeans(vecs, num_clusters)  
47.       
48.     for i in xrange(shape(clustAssing)[0]):  
49.         print clustAssing[i,0]  
50.       
51.   
52. def randCent(dataSet, k):  
53.     n = shape(dataSet)[1]  
54.     centroids = mat(zeros((k,n)))#create centroid mat  
55.     for j in range(n):#create random cluster centers, within bounds of each dimension  
56.         minJ = min(dataSet[:,j])   
57.         rangeJ = float(max(dataSet[:,j]) - minJ)  
58.         centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))  
59.     return centroids  
60.   
61. def distEclud(vecA, vecB):  
62.     return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)  
63.   
64. def kMeans(dataSet, k):  
65.     m = shape(dataSet)[0]  
66.     clusterAssment = mat(zeros((m,2)))#create mat to assign data points to a centroid, also holds SE of each point  
67.     centroids = randCent(dataSet, k)  
68.     clusterChanged = True  
69.     while clusterChanged:  
70.         clusterChanged = False  
71.         for i in range(m):#for each data point assign it to the closest centroid  
72.             minDist = inf; minIndex = -1  
73.             for j in range(k):  
74.                 distJI = distEclud(centroids[j,:],dataSet[i,:])  
75.                 if distJI < minDist:  
76.                     minDist = distJI; minIndex = j  
77.             if clusterAssment[i,0] != minIndex: clusterChanged = True  
78.             clusterAssment[i,:] = minIndex,minDist**2  
79.         #print centroids  
80.         for cent in range(k):#recalculate centroids  
81.             ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster  
82.             centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean   
83.     return centroids, clusterAssment  
84.   
85.   
86. if __name__ == '__main__':  
87.     # 导入数据集  
88.     matf = 'E://data_sc//corel_50_NN_sym_distance.mat'  
89.     dataDic = scio.loadmat(matf)  
90.     data = dataDic['A']  
91.     # 谱聚类的过程  
92.     spectalCluster(data, 20, 18)  

2、网上提供的一个Matlab代码

[plain] view plain copy

1. function [cluster_labels evd_time kmeans_time total_time] = sc(A, sigma, num_clusters)  
2. %SC Spectral clustering using a sparse similarity matrix (t-nearest-neighbor).  
3. %  
4. %   Input  : A              : N-by-N sparse distance matrix, where  
5. %                             N is the number of data  
6. %            sigma          : sigma value used in computing similarity,  
7. %                             if 0, apply self-tunning technique  
8. %            num_clusters   : number of clusters  
9. %  
10. %   Output : cluster_labels : N-by-1 vector containing cluster labels  
11. %            evd_time       : running time for eigendecomposition  
12. %            kmeans_time    : running time for k-means  
13. %            total_time     : total running time  
14.   
15. %  
16. % Convert the sparse distance matrix to a sparse similarity matrix,  
17. % where S = exp^(-(A^2 / 2*sigma^2)).  
18. % Note: This step can be ignored if A is sparse similarity matrix.  
19. %  
20. disp('Converting distance matrix to similarity matrix...');  
21. tic;  
22. n = size(A, 1);  
23.   
24. if (sigma == 0) % Selftuning spectral clustering  
25.   % Find the count of nonzero for each column  
26.   disp('Selftuning spectral clustering...');  
27.   col_count = sum(A~=0, 1)';  
28.   col_sum = sum(A, 1)';  
29.   col_mean = col_sum ./ col_count;  
30.   [x y val] = find(A);  
31.   A = sparse(x, y, -val.*val./col_mean(x)./col_mean(y)./2);  
32.   clear col_count col_sum col_mean x y val;  
33. else % Fixed-sigma spectral clustering  
34.   disp('Fixed-sigma spectral clustering...');  
35.   A = A.*A;  
36.   A = -A/(2*sigma*sigma);  
37. end  
38.   
39. % Do exp function sequentially because of memory limitation  
40. num = 2000;  
41. num_iter = ceil(n/num);  
42. S = sparse([]);  
43. for i = 1:num_iter  
44.   start_index = 1 + (i-1)*num;  
45.   end_index = min(i*num, n);  
46.   S1 = spfun(@exp, A(:,start_index:end_index)); % sparse exponential func  
47.   S = [S S1];  
48.   clear S1;  
49. end  
50. clear A;  
51. toc;  
52.   
53. %  
54. % Do laplacian, L = D^(-1/2) * S * D^(-1/2)  
55. %  
56. disp('Doing Laplacian...');  
57. D = sum(S, 2) + (1e-10);  
58. D = sqrt(1./D); % D^(-1/2)  
59. D = spdiags(D, 0, n, n);  
60. L = D * S * D;  
61. clear D S;  
62. time1 = toc;  
63.   
64. %  
65. % Do eigendecomposition, if L =  
66. %   D^(-1/2) * S * D(-1/2)    : set 'LM' (Largest Magnitude), or  
67. %   I - D^(-1/2) * S * D(-1/2): set 'SM' (Smallest Magnitude).  
68. %  
69. disp('Performing eigendecomposition...');  
70. OPTS.disp = 0;  
71. [V, val] = eigs(L, num_clusters, 'LM', OPTS);  
72. time2 = toc;  
73.   
74. %  
75. % Do k-means  
76. %  
77. disp('Performing kmeans...');  
78. % Normalize each row to be of unit length  
79. sq_sum = sqrt(sum(V.*V, 2)) + 1e-20;  
80. U = V ./ repmat(sq_sum, 1, num_clusters);  
81. clear sq_sum V;  
82. cluster_labels = k_means(U, [], num_clusters);  
83. total_time = toc;  
84.   
85. %  
86. % Calculate and show time statistics  
87. %  
88. evd_time = time2 - time1  
89. kmeans_time = total_time - time2  
90. total_time  
91. disp('Finished!');  

[plain] view plain copy

1. function cluster_labels = k_means(data, centers, num_clusters)  
2. %K_MEANS Euclidean k-means clustering algorithm.  
3. %  
4. %   Input    : data           : N-by-D data matrix, where N is the number of data,  
5. %                               D is the number of dimensions  
6. %              centers        : K-by-D matrix, where K is num_clusters, or  
7. %                               'random', random initialization, or  
8. %                               [], empty matrix, orthogonal initialization  
9. %              num_clusters   : Number of clusters  
10. %  
11. %   Output   : cluster_labels : N-by-1 vector of cluster assignment  
12. %  
13. %   Reference: Dimitrios Zeimpekis, Efstratios Gallopoulos, 2006.  
14. %              http://scgroup.hpclab.ceid./scgroup/Projects/TMG/  
15.   
16. %  
17. % Parameter setting  
18. %  
19. iter = 0;  
20. qold = inf;  
21. threshold = 0.001;  
22.   
23. %  
24. % Check if with initial centers  
25. %  
26. if strcmp(centers, 'random')  
27.   disp('Random initialization...');  
28.   centers = random_init(data, num_clusters);  
29. elseif isempty(centers)  
30.   disp('Orthogonal initialization...');  
31.   centers = orth_init(data, num_clusters);  
32. end  
33.   
34. %  
35. % Double type is required for sparse matrix multiply  
36. %  
37. data = double(data);  
38. centers = double(centers);  
39.   
40. %  
41. % Calculate the distance (square) between data and centers  
42. %  
43. n = size(data, 1);  
44. x = sum(data.*data, 2)';  
45. X = x(ones(num_clusters, 1), :);  
46. y = sum(centers.*centers, 2);  
47. Y = y(:, ones(n, 1));  
48. P = X + Y - 2*centers*data';  
49.   
50. %  
51. % Main program  
52. %  
53. while 1  
54.   iter = iter + 1;  
55.   
56.   % Find the closest cluster for each data point  
57.   [val, ind] = min(P, [], 1);  
58.   % Sum up data points within each cluster  
59.   P = sparse(ind, 1:n, 1, num_clusters, n);  
60.   centers = P*data;  
61.   % Size of each cluster, for cluster whose size is 0 we keep it empty  
62.   cluster_size = P*ones(n, 1);  
63.   % For empty clusters, initialize again  
64.   zero_cluster = find(cluster_size==0);  
65.   if length(zero_cluster) > 0  
66.     disp('Zero centroid. Initialize again...');  
67.     centers(zero_cluster, :)= random_init(data, length(zero_cluster));  
68.     cluster_size(zero_cluster) = 1;  
69.   end  
70.   % Update centers  
71.   centers = spdiags(1./cluster_size, 0, num_clusters, num_clusters)*centers;  
72.   
73.   % Update distance (square) to new centers  
74.   y = sum(centers.*centers, 2);  
75.   Y = y(:, ones(n, 1));  
76.   P = X + Y - 2*centers*data';  
77.   
78.   % Calculate objective function value  
79.   qnew = sum(sum(sparse(ind, 1:n, 1, size(P, 1), size(P, 2)).*P));  
80.   mesg = sprintf('Iteration %d:\n\tQold=%g\t\tQnew=%g', iter, full(qold), full(qnew));  
81.   disp(mesg);  
82.   
83.   % Check if objective function value is less than/equal to threshold  
84.   if threshold >= abs((qnew-qold)/qold)  
85.     mesg = sprintf('\nkmeans converged!');  
86.     disp(mesg);  
87.     break;  
88.   end  
89.   qold = qnew;  
90. end  
91.   
92. cluster_labels = ind';  
93.   
94.   
95. %-----------------------------------------------------------------------------  
96. function init_centers = random_init(data, num_clusters)  
97. %RANDOM_INIT Initialize centroids choosing num_clusters rows of data at random  
98. %  
99. %   Input : data         : N-by-D data matrix, where N is the number of data,  
100. %                          D is the number of dimensions  
101. %           num_clusters : Number of clusters  
102. %  
103. %   Output: init_centers : K-by-D matrix, where K is num_clusters  
104. rand('twister', sum(100*clock));  
105. init_centers = data(ceil(size(data, 1)*rand(1, num_clusters)), :);  
106.   
107. function init_centers = orth_init(data, num_clusters)  
108. %ORTH_INIT Initialize orthogonal centers for k-means clustering algorithm.  
109. %  
110. %   Input : data         : N-by-D data matrix, where N is the number of data,  
111. %                          D is the number of dimensions  
112. %           num_clusters : Number of clusters  
113. %  
114. %   Output: init_centers : K-by-D matrix, where K is num_clusters  
115.   
116. %  
117. % Find the num_clusters centers which are orthogonal to each other  
118. %  
119. Uniq = unique(data, 'rows'); % Avoid duplicate centers  
120. num = size(Uniq, 1);  
121. first = ceil(rand(1)*num); % Randomly select the first center  
122. init_centers = zeros(num_clusters, size(data, 2)); % Storage for centers  
123. init_centers(1, :) = Uniq(first, :);  
124. Uniq(first, :) = [];  
125. c = zeros(num-1, 1); % Accumalated orthogonal values to existing centers for non-centers  
126. % Find the rest num_clusters-1 centers  
127. for j = 2:num_clusters  
128.   c = c + abs(Uniq*init_centers(j-1, :)');  
129.   [minimum, i] = min(c); % Select the most orthogonal one as next center  
130.   init_centers(j, :) = Uniq(i, :);  
131.   Uniq(i, :) = [];  
132.   c(i) = [];  
133. end  
134. clear c Uniq;  

个人的一点认识：谱聚类的过程相当于先进行一个非线性的降维，然后在这样的低维空间中再利用聚类的方法进行聚类。

欢迎大家一起讨论，如有问题欢迎留言，欢迎大家转载。

参考

1、从拉普拉斯矩阵说到谱聚类(http://blog.csdn.NET/v_july_v/article/details/40738211)

2、谱聚类(spectral clustering)(http://www.cnblogs.com/FengYan/archive/2012/06/21/2553999.html)

3、谱聚类算法(Spectral Clustering)(http://www.cnblogs.com/sparkwen/p/3155850.html)