An example of Spectral clustering

494 days ago by Buyan

1.Plot the original data.

data(ruspini, package="cluster") ruspini <- ruspini[sample(1:nrow(ruspini)),] plot(ruspini) my.data <- scale(ruspini) plot(my.data) 
       

2. Find the similarity matrix S (if two data are very close, then the similarity will be closer to 1.)

s <- function(x1, x2, alpha=1) { exp(- alpha * norm(as.matrix(x1-x2), type="F")) } make.similarity <- function(my.data, similarity) { N <- nrow(my.data) S <- matrix(rep(NA,N^2), ncol=N) for(i in 1:N) { for(j in 1:N) { S[i,j] <- similarity(my.data[i,], my.data[j,]) } } S } S <- make.similarity(my.data, s) round(S[1:12,1:12],2) 
       
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,] 1.00 0.15 0.42 0.22 0.50 0.73 0.18 0.10 0.11  0.65  0.75  0.19
 [2,] 0.15 1.00 0.06 0.06 0.07 0.13 0.05 0.08 0.10  0.11  0.11  0.74
 [3,] 0.42 0.06 1.00 0.24 0.74 0.48 0.24 0.07 0.08  0.47  0.53  0.08
 [4,] 0.22 0.06 0.24 1.00 0.20 0.30 0.74 0.24 0.23  0.16  0.21  0.07
 [5,] 0.50 0.07 0.74 0.20 1.00 0.51 0.19 0.06 0.07  0.62  0.66  0.10
 [6,] 0.73 0.13 0.48 0.30 0.51 1.00 0.25 0.12 0.14  0.53  0.66  0.15
 [7,] 0.18 0.05 0.24 0.74 0.19 0.25 1.00 0.20 0.18  0.14  0.18  0.05
 [8,] 0.10 0.08 0.07 0.24 0.06 0.12 0.20 1.00 0.77  0.06  0.08  0.08
 [9,] 0.11 0.10 0.08 0.23 0.07 0.14 0.18 0.77 1.00  0.07  0.09  0.10
[10,] 0.65 0.11 0.47 0.16 0.62 0.53 0.14 0.06 0.07  1.00  0.79  0.15
[11,] 0.75 0.11 0.53 0.21 0.66 0.66 0.18 0.08 0.09  0.79  1.00  0.14
[12,] 0.19 0.74 0.08 0.07 0.10 0.15 0.05 0.08 0.10  0.15  0.14  1.00
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,] 1.00 0.15 0.42 0.22 0.50 0.73 0.18 0.10 0.11  0.65  0.75  0.19
 [2,] 0.15 1.00 0.06 0.06 0.07 0.13 0.05 0.08 0.10  0.11  0.11  0.74
 [3,] 0.42 0.06 1.00 0.24 0.74 0.48 0.24 0.07 0.08  0.47  0.53  0.08
 [4,] 0.22 0.06 0.24 1.00 0.20 0.30 0.74 0.24 0.23  0.16  0.21  0.07
 [5,] 0.50 0.07 0.74 0.20 1.00 0.51 0.19 0.06 0.07  0.62  0.66  0.10
 [6,] 0.73 0.13 0.48 0.30 0.51 1.00 0.25 0.12 0.14  0.53  0.66  0.15
 [7,] 0.18 0.05 0.24 0.74 0.19 0.25 1.00 0.20 0.18  0.14  0.18  0.05
 [8,] 0.10 0.08 0.07 0.24 0.06 0.12 0.20 1.00 0.77  0.06  0.08  0.08
 [9,] 0.11 0.10 0.08 0.23 0.07 0.14 0.18 0.77 1.00  0.07  0.09  0.10
[10,] 0.65 0.11 0.47 0.16 0.62 0.53 0.14 0.06 0.07  1.00  0.79  0.15
[11,] 0.75 0.11 0.53 0.21 0.66 0.66 0.18 0.08 0.09  0.79  1.00  0.14
[12,] 0.19 0.74 0.08 0.07 0.10 0.15 0.05 0.08 0.10  0.15  0.14  1.00

3. Compute an adjacency matrix A based on S.

This is usually done by applying a k-nearest neighboor filter to build a representation of a graph connecting just the closest dataset points.

make.affinity <- function(S, n.neighboors=2) { N <- length(S[,1]) if (n.neighboors >= N) { # fully connected A <- S } else { A <- matrix(rep(0,N^2), ncol=N) for(i in 1:N) { # for each line # only connect to those points with larger similarity best.similarities <- sort(S[i,], decreasing=TRUE)[1:n.neighboors] for (s in best.similarities) { j <- which(S[i,] == s) A[i,j] <- S[i,j] A[j,i] <- S[i,j] # to make an undirected graph, ie, the matrix becomes symmetric } } } A } A <- make.affinity(S, 15) # use 15 neighboors (includes self) round(A[1:11,1:11],3) 
       
       [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 1.000    0 0.000 0.000 0.504 0.730 0.000 0.000 0.000 0.655 0.754
 [2,] 0.000    1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [3,] 0.000    0 1.000 0.000 0.745 0.000 0.000 0.000 0.000 0.000 0.000
 [4,] 0.000    0 0.000 1.000 0.000 0.000 0.744 0.000 0.000 0.000 0.000
 [5,] 0.504    0 0.745 0.000 1.000 0.505 0.000 0.000 0.000 0.625 0.662
 [6,] 0.730    0 0.000 0.000 0.505 1.000 0.000 0.000 0.000 0.529 0.662
 [7,] 0.000    0 0.000 0.744 0.000 0.000 1.000 0.000 0.000 0.000 0.000
 [8,] 0.000    0 0.000 0.000 0.000 0.000 0.000 1.000 0.774 0.000 0.000
 [9,] 0.000    0 0.000 0.000 0.000 0.000 0.000 0.774 1.000 0.000 0.000
[10,] 0.655    0 0.000 0.000 0.625 0.529 0.000 0.000 0.000 1.000 0.790
[11,] 0.754    0 0.000 0.000 0.662 0.662 0.000 0.000 0.000 0.790 1.000
       [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 1.000    0 0.000 0.000 0.504 0.730 0.000 0.000 0.000 0.655 0.754
 [2,] 0.000    1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [3,] 0.000    0 1.000 0.000 0.745 0.000 0.000 0.000 0.000 0.000 0.000
 [4,] 0.000    0 0.000 1.000 0.000 0.000 0.744 0.000 0.000 0.000 0.000
 [5,] 0.504    0 0.745 0.000 1.000 0.505 0.000 0.000 0.000 0.625 0.662
 [6,] 0.730    0 0.000 0.000 0.505 1.000 0.000 0.000 0.000 0.529 0.662
 [7,] 0.000    0 0.000 0.744 0.000 0.000 1.000 0.000 0.000 0.000 0.000
 [8,] 0.000    0 0.000 0.000 0.000 0.000 0.000 1.000 0.774 0.000 0.000
 [9,] 0.000    0 0.000 0.000 0.000 0.000 0.000 0.774 1.000 0.000 0.000
[10,] 0.655    0 0.000 0.000 0.625 0.529 0.000 0.000 0.000 1.000 0.790
[11,] 0.754    0 0.000 0.000 0.662 0.662 0.000 0.000 0.000 0.790 1.000

4. The degree matrix.

D <- diag(apply(A, 1, sum)) # sum rows round(D[1:11,1:11],2) 
       
       [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 10.61 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [2,]  0.00 9.18  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [3,]  0.00 0.00 11.71 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [4,]  0.00 0.00  0.00 9.72  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [5,]  0.00 0.00  0.00 0.00 17.2  0.00  0.00  0.00  0.00  0.00  0.00
 [6,]  0.00 0.00  0.00 0.00  0.0 10.03  0.00  0.00  0.00  0.00  0.00
 [7,]  0.00 0.00  0.00 0.00  0.0  0.00 13.17  0.00  0.00  0.00  0.00
 [8,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00 11.78  0.00  0.00  0.00
 [9,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00 10.76  0.00  0.00
[10,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00 11.26  0.00
[11,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00 12.11
       [,1] [,2]  [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 10.61 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [2,]  0.00 9.18  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [3,]  0.00 0.00 11.71 0.00  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [4,]  0.00 0.00  0.00 9.72  0.0  0.00  0.00  0.00  0.00  0.00  0.00
 [5,]  0.00 0.00  0.00 0.00 17.2  0.00  0.00  0.00  0.00  0.00  0.00
 [6,]  0.00 0.00  0.00 0.00  0.0 10.03  0.00  0.00  0.00  0.00  0.00
 [7,]  0.00 0.00  0.00 0.00  0.0  0.00 13.17  0.00  0.00  0.00  0.00
 [8,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00 11.78  0.00  0.00  0.00
 [9,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00 10.76  0.00  0.00
[10,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00 11.26  0.00
[11,]  0.00 0.00  0.00 0.00  0.0  0.00  0.00  0.00  0.00  0.00 12.11

5. Find the Laplacian matrix

L <- D - A round(L[1:12,1:12],1) 
       
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]  9.6  0.0  0.0  0.0 -0.5 -0.7  0.0  0.0  0.0  -0.7  -0.8   0.0
 [2,]  0.0  8.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0  -0.7
 [3,]  0.0  0.0 10.7  0.0 -0.7  0.0  0.0  0.0  0.0   0.0   0.0   0.0
 [4,]  0.0  0.0  0.0  8.7  0.0  0.0 -0.7  0.0  0.0   0.0   0.0   0.0
 [5,] -0.5  0.0 -0.7  0.0 16.2 -0.5  0.0  0.0  0.0  -0.6  -0.7   0.0
 [6,] -0.7  0.0  0.0  0.0 -0.5  9.0  0.0  0.0  0.0  -0.5  -0.7   0.0
 [7,]  0.0  0.0  0.0 -0.7  0.0  0.0 12.2  0.0  0.0   0.0   0.0   0.0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 10.8 -0.8   0.0   0.0   0.0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -0.8  9.8   0.0   0.0   0.0
[10,] -0.7  0.0  0.0  0.0 -0.6 -0.5  0.0  0.0  0.0  10.3  -0.8   0.0
[11,] -0.8  0.0  0.0  0.0 -0.7 -0.7  0.0  0.0  0.0  -0.8  11.1   0.0
[12,]  0.0 -0.7  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0   8.9
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]  9.6  0.0  0.0  0.0 -0.5 -0.7  0.0  0.0  0.0  -0.7  -0.8   0.0
 [2,]  0.0  8.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0  -0.7
 [3,]  0.0  0.0 10.7  0.0 -0.7  0.0  0.0  0.0  0.0   0.0   0.0   0.0
 [4,]  0.0  0.0  0.0  8.7  0.0  0.0 -0.7  0.0  0.0   0.0   0.0   0.0
 [5,] -0.5  0.0 -0.7  0.0 16.2 -0.5  0.0  0.0  0.0  -0.6  -0.7   0.0
 [6,] -0.7  0.0  0.0  0.0 -0.5  9.0  0.0  0.0  0.0  -0.5  -0.7   0.0
 [7,]  0.0  0.0  0.0 -0.7  0.0  0.0 12.2  0.0  0.0   0.0   0.0   0.0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 10.8 -0.8   0.0   0.0   0.0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -0.8  9.8   0.0   0.0   0.0
[10,] -0.7  0.0  0.0  0.0 -0.6 -0.5  0.0  0.0  0.0  10.3  -0.8   0.0
[11,] -0.8  0.0  0.0  0.0 -0.7 -0.7  0.0  0.0  0.0  -0.8  11.1   0.0
[12,]  0.0 -0.7  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0   8.9

6. Find the eigenvalues of Laplacian matrix to determine the value of “k”

evL <- eigen(L, symmetric=TRUE) plot(1:10, rev(evL$values)[1:10], log="y") plot(1:10, rev(evL$values)[1:10], log="y") 
       
경고 메시지가 손실되었습니다
In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot
경고 메시지가 손실되었습니다
In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot
경고 메시지가 손실되었습니다
In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot
경고 메시지가 손실되었습니다
In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot

7. Find the corresponding eigenvectors of k smallest eigenvalues and use k-means clustering to find the appropriate clusters.

k <- 4 Z <- evL$vectors[,(ncol(evL$vectors)-k+1):ncol(evL$vectors)] km <- kmeans(Z, centers=k) plot(my.data, col=km$cluster) plot(my.data, col=km$cluster) 
       

k <- 3 Z <- evL$vectors[,(ncol(evL$vectors)-k+1):ncol(evL$vectors)] km <- kmeans(Z, centers=k) plot(my.data, col=km$cluster) plot(my.data, col=km$cluster)