上述的“截”函数通常会将图分割成一个点和其余
个点。
4、其他的“截”函数的表现形式
性质3的证明：
4、不同的Laplacian矩阵
除了上述的拉普拉斯矩阵，还有规范化的 Laplacian 矩阵形式：
四、Laplacian矩阵与谱聚类中的优化函数的关系
1、由Laplacian矩阵到“截”函数
1、自己实现的一个
#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代码
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!');
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;
个人的一点认识：谱聚类的过程相当于先进行一个非线性的降维，然后在这样的低维空间中再利用聚类的方法进行聚类。
欢迎大家一起讨论，如有问题欢迎留言，欢迎大家转载。
参考
1、从拉普拉斯矩阵说到谱聚类
2、谱聚类(spectral clustering)
3、谱聚类算法(Spectral Clustering