了凡春秋USTC
谱聚类
http://chunqiu.blog.ustc.edu.cn/?p=505
最近忙着写文章，好久不写博客了。最近看到一个聚类方法--谱聚类，号称现代聚类方法，看到它简洁的公式推导、实现代码，不禁要尝试一把。关于它的理论，google一搜有很多博客讲，这里就不赘述了，反正最后还是归结为一个SVD分解问题，参考网址如下
看到下面的图，乐了
上网搜索，看到了几个相关的工具箱，下面一一介绍：
1.这个是加利福尼亚大学一个实验室为他们的论文《Parallel Spectral Clustering in Disributed Systems》而写的测试代码，这篇文章发在IEEE Transactions on Pattern Analysis and Machine Intelligence上，机器学习方面的顶级期刊。不过看了看作者，怎么感觉都是中国人啊，下载地址见http://alumni.cs.ucsb.edu/~wychen/sc.html#FAQ
压缩包内有Readme，有工具箱的使用介绍，不过由于是为文章效果而写的，所以感觉并不能称得上工具箱吧，我自己编了个测试例子试了试，效果、速度确实不错：
原始待聚类数据
聚类后
代码如下：
% 谱聚类算法测试
clc
clear
close all
% 构造原始数据
ang = 0:pi/200:pi/2*3;
x1 = sin(ang) + randn(size(ang))*0.1;
y1 = cos(ang) + randn(size(ang))*0.1;
x2 = 2.3*sin(ang+pi) + randn(size(ang))*0.1;
y2 = 2.3*cos(ang) + randn(size(ang))*0.1;
x3 = 6*sin(ang+pi/2*3+pi/6) + randn(size(ang))*0.1;
y3 = 5*cos(ang+pi/2*3) + randn(size(ang))*0.1;
X = [x1, x2, x3; y1, y2, y3]';
figure
plot(X(:, 1), X(:, 2), '.')
axis equal
% 谱聚类
num_neighbors = 20;
block_size = 5;
gen_nn_distance(X, num_neighbors, block_size, 0); % 计算相似矩阵，保存到NN_sym_distance.mat文件中
filename = [num2str(num_neighbors), '_NN_sym_distance'];
load(filename)
num_clusters = 3;
sigma = 10;
cluster_labels = sc(A, sigma, num_clusters); % 谱聚类
figure
plot(X(cluster_labels == 1, 1), X(cluster_labels == 1, 2), 'r.', ...
X(cluster_labels == 2, 1), X(cluster_labels == 2, 2), 'b.', ...
X(cluster_labels == 3, 1), X(cluster_labels == 3, 2), 'g.')
axis equal
2.有个matlab exchange上的工具箱，是德国斯图加特大学一个博士整齐其博士论文的研究而写的GUI，下载链接如下（压缩包内有其博士论文）
同样，Readme里有详细使用说明，只要运行MAIN.m文件即可弹出GUI。运行可能会出个警告
影响不大，Readme里是这么说的：
recommended: export_fig configured to save to pdf files (Mathworks File Exchange ID #23629)
如果你不用保存图片，不用理会。
GUI如下
对其给的twomoons的例子做测试，结果如下
额，是不是有点自己打脸的赶脚！自己给的例子，效果也不够完美。。。估计调整一些参数能得到理想的效果吧！
Matlab exchange上还有个谱聚类的文件，提供了谱聚类各种算法的简单实现，对于学习谱聚类具体实现挺有帮助！
3.很多介绍谱聚类的博客都有提到算法的matlab实现，因为实现代码太简单了，如下
% 谱聚类简单实现
function idx = spectral_clustering(W, k)
D = diag(sum(W));
L = D-W;
opt = struct('issym', true, 'isreal', true);
[V, ~] = eigs(L, D, k, 'SM', opt);
idx = kmeans(V, k);
end
谱聚类复杂的理论最终归结为代码就这几句，当然这是考虑最简化的情况，如果考虑细节，就得靠上面提到的工具箱了。对于这个代码，我做了测试，效果还行（使用第一个代码里的计算相似矩阵的函数），不过多次测试的效果说明，不同数据对效果是有影响的，有时就得不到好的分离效果。
% 谱聚类算法测试
clc
clear
close all
% 构造原始数据
ang = 0:pi/200:pi/2*3;
x1 = sin(ang) + randn(size(ang))*0.1;
y1 = cos(ang) + randn(size(ang))*0.1;
x2 = 2.3*sin(ang+pi) + randn(size(ang))*0.1;
y2 = 2.3*cos(ang) + randn(size(ang))*0.1;
x3 = 6*sin(ang+pi/2*3+pi/6) + randn(size(ang))*0.1;
y3 = 5*cos(ang+pi/2*3) + randn(size(ang))*0.1;
X = [x1, x2, x3; y1, y2, y3]';
figure
plot(X(:, 1), X(:, 2), '.')
axis equal
% 谱聚类
num_neighbors = 20;
block_size = 5;
gen_nn_distance(X, num_neighbors, block_size, 0); % 计算相似矩阵，保存到NN_sym_distance.mat文件中
filename = [num2str(num_neighbors), '_NN_sym_distance'];
load(filename)
idx = spectral_clustering(A, 3);
figure
plot(X(idx == 1, 1), X(idx == 1, 2), 'r.', ...
X(idx == 2, 1), X(idx == 2, 2), 'b.', ...
X(idx == 3, 1), X(idx == 3, 2), 'g.')
axis equal
如果对那个博士论文的例子做测试，效果如下
貌似这个图确实很难分，这个效果不如那个德国博士的效果。
% 谱聚类算法测试
clc
clear
close all
% 构造原始数据
M = csvread('twomoon-2d-50s.test.csv');
X = M(:, 2:3);
figure
plot(X(:, 1), X(:, 2), '.')
axis equal
% 谱聚类
num_neighbors = 20;
block_size = 5;
gen_nn_distance(X, num_neighbors, block_size, 0); % 计算相似矩阵，保存到NN_sym_distance.mat文件中
filename = [num2str(num_neighbors), '_NN_sym_distance'];
load(filename)
idx = spectral_clustering(A, 2);
figure
plot(X(idx == 1, 1), X(idx == 1, 2), 'r.', ...
X(idx == 2, 1), X(idx == 2, 2), 'b.')
axis equal
如果用第一个工具箱的方法再测试这个例子
效果也就那样！
% 谱聚类算法测试
clc
clear
close all
% 构造原始数据
M = csvread('twomoon-2d-50s.test.csv');
X = M(:, 2:3);
figure
plot(X(:, 1), X(:, 2), '.')
axis equal
% 谱聚类
num_neighbors = 20;
block_size = 5;
gen_nn_distance(X, num_neighbors, block_size, 0); % 计算相似矩阵，保存到NN_sym_distance.mat文件中
filename = [num2str(num_neighbors), '_NN_sym_distance'];
load(filename)
num_clusters = 2;
sigma = 10;
cluster_labels = sc(A, sigma, num_clusters); % 谱聚类
figure
plot(X(cluster_labels == 1, 1), X(cluster_labels == 1, 2), 'r.', ...
X(cluster_labels == 2, 1), X(cluster_labels == 2, 2), 'b.')
axis equal