谱聚类算法—Matlab代码

% =========================================================================
% 算 法 名 称: Spectral Clustering Algorithm
% 编 码 作 者: Lee Wen-Tsao
% 编 码 邮 箱: liwenchao36@163.com
% 输 入 参 数: 
%             W ---> 邻接矩阵
%             k ---> 簇数目
%             t ---> 拉普拉斯矩阵归一化处理类型
% =========================================================================
%% step1: 清理运行环境
clc;
clear;
close all;

%% step2: 读入数据
Iris = uiimport('iris.data');
Iris = cellfun(@(x) regexp(x,',','split'), Iris.iris,'UniformOutput',false);
data = cellfun(@(x) x(:,1:4),Iris,'UniformOutput',false);
data = str2double(reshape([data{:}],4,150)');

%% step3: 构造相似矩阵
H = pdist2(data, data, 'euclidean');
W = 1-exp(-(H.^2)./2);
triu_W = triu(W, 0)./(sum(triu(W, 0),2) + eps);
W = triu_W' + triu_W;

%% step4: 计算度矩阵
d = sum(W, 2);                           % 对W进行列求和
D = sparse(1:size(W,1), 1:size(W,2), d); % 然后将d中的元素放到对角线上

%% step5: 计算拉普拉斯矩阵
% 1.未标准化的拉普拉斯矩阵
L = D - W; 

% 2.正则拉普拉斯矩阵
t = 'Symmetric';
switch t
    case 'RandomWalk'
        % 避免除以0
        d(d==0) = eps;
        % 计算D的逆
        D = spdiags(1./d, 0, size(D, 1), size(D, 2));
        % 随机游走正则化拉普拉斯矩阵
        L = D*L;
    case 'Symmetric'
        % 避免除以0
        d(d==0) = eps;
        % 计算D^(1/2)
        D = spdiags(1./(d.^0.5), 0, size(D, 1), size(D, 2));
        % 对称正则化拉普拉斯矩阵
        L = D*L*D;
end

%% step5: 特征值和特征向量
% 1.V表示特征向量；lamda表示特征值
k = 3;
[U, lamda] = eigs(L, k, 'smallestabs');  % 不能这么求特征向量，特征向量有重数

if strcmp('Symmetric', t)
    % 对称拉普拉斯矩阵单位化
    U = bsxfun(@rdivide, U, sqrt(sum(U.^2, 2)));
end 
%% step6: 使用kmeans对函数分类
% 0. 问题定义
labels = zeros(size(U,1),1);
errors = zeros(k, 1);
expose = 1;

% 1. 初始化簇心
loc = randperm(size(U,1));
centroids = U(loc(1:k),:);
% 2. 迭代
N_iter = 1000;
for it=1:N_iter
    for i=1:size(U,1)
        dists = sqrt(sum((U(i,:) - centroids).^2, 2));     % 计算每个数据到k个簇心的距离
        [distMin, idx] = min(dists);                       % 寻找距离每个簇心的最小距离
        labels(i,:) = idx;                                 % 给每个数据标注
    end
    
    % 3. 计算误差率
    for j=1:k
        errors(j, :) = sum(sqrt(sum((U(j==labels, :)- centroids(j, :)).^2, 2)));
    end
    
    % 4. 可视化
    if expose
        disp(sum(errors));
    end
    
    % 5. 更新簇心
    for j=1:k
        centroids(j,:) = mean(U((j==labels),:),1);
    end
end

思考：

1. 为什么要使用拉普拉斯正则化？
   拉普拉斯正则化过程有两个：
   （1）随机游走拉普拉斯正则化
   （2）对称拉普拉斯正则化
2. 上述拉普拉斯正则化的理论基础是什么？
3. 这种降维方式的原理是什么呢？
4. 这种聚类算法效果为啥没有论文里说的那么好，问题出现在哪里？