譜聚類算法—Matlab代碼

% =========================================================================
% 算 法 名 稱: Spectral Clustering Algorithm
% 編 碼 作 者: Lee Wen-Tsao
% 編 碼 郵 箱: [e-mail]
% 輸 入 參 數: 
%             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. 這種聚類算法效果為啥沒有論文里說的那么好，問題出現在哪里？