Toy visualization of spectral clustering

clear all 
close all
set(0,'DefaultFigureWindowStyle','docked')

rng('default') % For reproducibility

% Parameters for data generation
N = 300;  % Size of each cluster
radii = [2 5 8]; % the radii of each disk, the number of radii defines the
                 % number of disks (i.e. number of classes)
K = length(radii);
theta = linspace(0,2*pi,N)'; %variable for parametric definition of disks data'

X = []; % Noisy 2-D circular data set
for i = 1 : length(radii)
    Xi = radii(i)*[cos(theta),sin(theta)]+ rand(N,2);
    X = [X;Xi];
end

figure
title('ground truth')
hold on
for i = 1:K
    I = (i-1)*N+1 : i*N;
    plot(X(I,1),X(I,2),'o');
end
hold off
pbaspect([1 1 1])
xticks('')
yticks('')


[IND, C] = kmeans(X,K);
figure
title('kmeans')
hold on
for i = 1:K
    I = (IND==i);
    plot(X(I,1),X(I,2),'o');
    plot(C(i,1),C(i,2),'rx','markersize',15);
end
hold off
pbaspect([1 1 1])
xticks('')
yticks('')


%% build the knn similarity matrix
A = zeros(length(X));
k = 15; %number of nearest neighbors
for i = 1 : length(X)
    for j = i+1 : length(X)
        A(i,j) = exp(-norm(X(i,:)-X(j,:),2));
    end
    [val,idx] = sort(A(i,:),'descend');
    A(i,idx(k+1:end)) = 0;
end
A = spones(A+A');


%% build the Laplacian matrix'
D = diag(sum(A));
L = D-A;
[U,~] = eigs(L,K,'SM');


%% Define the new embedding, the first colum is redoundant
NewX = U;

%% if K = 2 or K=3 we can draw the spectral embedding
if K<=3
    figure
    title('new embedding')
    hold on
    for i = 1:K
        I = (i-1)*N+1 : i*N;
        if K == 3
            plot(NewX(I,2),NewX(I,3),'o');
        elseif K==2
            plot(NewX(I,1),NewX(I,2),'o');
        end
    end
    hold off
    pbaspect([1 1 1])
    xticks('')
    yticks('')
end

%% Apply kmeans to the new embedding
[IND, C] = kmeans(NewX,K);

%% if K = 2 or K=3 we can draw the effect of kmeans on the spectral embedding
if K <= 3
    figure
    title('kmeans on spectral embedding')
    hold on
    for i = 1:K
        I = (IND==i);
        if K == 3
            plot(NewX(I,2),NewX(I,3),'o');
        elseif K==2
            plot(NewX(I,1),NewX(I,2),'o');
            plot(C(i,1),C(i,2),'rx','markersize',15);
        end 
    end
    hold off
    pbaspect([1 1 1])
    xticks('')
    yticks('')
end


%% Plot the effect of spectral clustering
figure
title('spectral clustering')
hold on
for i = 1:K
    I = (IND==i);
    plot(X(I,1),X(I,2),'o');
end
hold off
pbaspect([1 1 1])
xticks('')
yticks('')