forked from cswhjiang/Transfer-Spectral-Clustering
-
Notifications
You must be signed in to change notification settings - Fork 0
/
kmeans_cl.m
146 lines (132 loc) · 4.44 KB
/
kmeans_cl.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
function [cluster_labels,qnew] = kmeans_cl(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
%
% Author : Dimitrios Zeimpekis, Efstratios Gallopoulos, 2006.
% http://scgroup.hpclab.ceid.upatras.gr/scgroup/Projects/TMG/
%
% Modified : Wen-Yen Chen (wychen@alumni.cs.ucsb.edu)
% Chih-Jen Lin (cjlin@csie.ntu.edu.tw)
% Frank Lin (frank@cs.cmu.edu)
%
% Parameter setting
%
iter = 0;
qold = inf;
threshold = 0.001;
%
% Check if with initial centers
%
if strcmp(centers, 'random')
cdisp('Random initialization...');
centers = random_init(data, num_clusters);
elseif isempty(centers)
cdisp('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
cdisp('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));
cdisp(mesg);
% Check if objective function value is less than/equal to threshold
if threshold >= abs((qnew-qold)/qold)
mesg = sprintf('\nkmeans converged!');
cdisp(mesg);
break;
end
qold = qnew;
end
cluster_labels = ind';
end
%-----------------------------------------------------------------------------
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)), :);
end
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;
end
function cdisp(msg)
if false
disp(msg);
end
end