2 How to use spectral clustering
This chapter shows how to actually run spectral clustering, in R and in Python, on the kind of inputs you are likely to have in practice: ordinary feature data, a similarity matrix you already have, and a fused network from Similarity Network Fusion.
We keep using the same 300-observation crescent dataset from the intuition chapter, stored once in data/two-moons.csv.
The files used in the examples are available here:
- two-moons.csv — the crescent data used throughout the chapter
- snf-fused.csv — the fused patient similarity matrix in Example C
- snf-simulated-patients.csv — the matching patient subtype labels
Save the files in a data folder beside your R script or Python notebook so that the paths in the examples work unchanged.
2.1 Which packages, and why
For Python, there is really one obvious choice: sklearn.cluster.SpectralClustering. It is actively maintained, it accepts either raw features or a precomputed affinity matrix through the same interface, and it is what most people mean when they say “spectral clustering in Python.” All Python examples in this chapter use it, from scikit-learn 1.9.
R does not have one single dominant implementation, so this chapter uses two, each suited to a different situation:
kernlab::speccbuilds its own similarity matrix from an RBF kernel with an automatically chosen bandwidth, and is the natural choice when you are starting from raw feature data, similar to sklearn’s default behavior.SNFtool::spectralClusteringexpects a similarity matrix you already built, and offers a choice of three Laplacian conventions. Because it already appears in the companion Similarity Network Fusion course, it is the natural choice whenever the input is a precomputed similarity matrix, including a fused SNF network.
Other maintained options exist and are worth knowing about even though we do not use them here. RSpectra provides a fast, memory-efficient way to compute only the eigenvectors you actually need, which matters once a similarity matrix gets large. igraph implements a variety of graph-clustering algorithms and is a convenient way to inspect connected components. We use igraph briefly in this chapter for exactly that, checking whether a similarity graph is fully connected.
2.3 Example A: ordinary feature data
Let’s first confirm what the intuition chapter showed us visually: k-means struggles on this data, and spectral clustering does not.
set.seed(42)
km <- kmeans(X, centers = 2, nstart = 10)
ari(km$cluster, y_true)[1] 0.261from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
kmeans = KMeans(n_clusters=2, n_init=10, random_state=42)
kmeans_labels = kmeans.fit_predict(X)
adjusted_rand_score(y_true, kmeans_labels)0.261Cluster labels are arbitrary names: “cluster 1” in one run can be “cluster 2”, or in Python “cluster 0”, in another. The adjusted Rand index (ARI) compares two partitions without caring which name each group was given: 1 means perfect agreement, values near 0 mean roughly chance-level agreement. 0.261 confirms what Figure 1.2 already showed: k-means is not finding the two crescents. Now the spectral version.
library(kernlab)
set.seed(42)
spectral_fit <- specc(X, centers = 2)
ari(spectral_fit, y_true)
table(as.integer(spectral_fit))[1] 1
150 150from sklearn.cluster import SpectralClustering
spectral = SpectralClustering(
n_clusters=2,
affinity="nearest_neighbors",
n_neighbors=10,
random_state=42,
)
spectral_labels = spectral.fit_predict(X)
adjusted_rand_score(y_true, spectral_labels)1.0Both give a perfect ARI of 1, and both split the 300 observations into two groups of 150, exactly matching the two crescents. It is worth being honest, though, that specc and SpectralClustering got there in different ways. specc built a dense similarity between every pair of points using an RBF kernel, and automatically searched over the kernel bandwidth for one that produces well-separated clusters. SpectralClustering with affinity="nearest_neighbors" instead connects each point only to its 10 nearest neighbours, the same kind of network we built by hand in the intuition chapter. These are two different affinity constructions arriving at the same answer on this dataset. They will not always agree this closely.
2.3.1 What the inputs and outputs actually are
X: an n × p matrix of raw features, one row per observation. This is what bothspeccandSpectralClustering(affinity="nearest_neighbors")expect when you give them feature data directly.centers(R) /n_clusters(Python): how many clusters to look for. Spectral clustering does not discover this number on its own any more than k-means does. You supply it. We look at a data-driven way to get a sense of a reasonable value later in this chapter.- The output of both functions is a plain vector of integer cluster labels, one per observation, in the same order as the rows of
X. Nothing more.
2.3.2 Parameters worth understanding before you tune them
affinity (Python) controls how the similarity graph is built from raw features. "nearest_neighbors" builds the kind of network used throughout this chapter. "rbf" (the default if you do not set affinity) instead computes a similarity between every pair of points using a Gaussian kernel, controlled by gamma. There is no single correct choice. It depends on whether your clusters are better described by local neighbourhoods or by a smooth global similarity.
n_neighbors (Python, only used when affinity="nearest_neighbors") is exactly the neighbourhood size from Figure 1.10. Too small and the graph fragments. Too large and unrelated groups blend together. The next section shows this numerically.
kpar (R, kernlab::specc) controls the kernel used to build the similarity matrix. The default, kpar = "automatic", is what we used above: it searches over a range of RBF bandwidths and keeps the one that produces the tightest k-means fit in the resulting spectral coordinates. You can also pass kpar = "local" for a bandwidth that adapts to each point’s local neighbourhood, similar in spirit to affinity="nearest_neighbors".
random_state (Python) / a preceding set.seed() (R) affects two things: which eigensolver starting vector is used, and how k-means is initialized on the resulting coordinates. Spectral clustering is not fully deterministic without fixing this, because the final k-means step depends on where its centroids start.
2.4 Example E: the graph you build is the graph you get
The previous chapter showed this visually with three neighbourhood sizes. Let’s now see the effect numerically.
library(igraph)
# Same k-NN convention scikit-learn uses internally: each point is its own
# neighbour, and the directed graph is symmetrized by averaging.
knn_affinity <- function(X, k) {
n <- nrow(X)
Dm <- as.matrix(dist(X))
C <- matrix(0, n, n)
for (i in 1:n) {
nearest <- order(Dm[i, ])[1:(k + 1)]
C[i, nearest] <- 1
}
0.5 * (C + t(C))
}
for (k in c(2, 5, 10, 15, 30, 100)) {
Wk <- knn_affinity(X, k = k)
n_comp <- components(graph_from_adjacency_matrix(
Wk, mode = "undirected", weighted = TRUE))$no
labels_k <- SNFtool::spectralClustering(Wk, K = 2, type = 3)
cat(sprintf("n_neighbors=%4d connected_components=%3d ARI=%.3f\n",
k, n_comp, ari(labels_k, y_true)))
}n_neighbors= 2 connected_components= 20 ARI=0.013
n_neighbors= 5 connected_components= 2 ARI=1.000
n_neighbors= 10 connected_components= 1 ARI=1.000
n_neighbors= 15 connected_components= 1 ARI=1.000
n_neighbors= 30 connected_components= 1 ARI=0.716
n_neighbors= 100 connected_components= 1 ARI=0.443from scipy.sparse.csgraph import connected_components
for k in [2, 5, 10, 15, 30, 100]:
sc_k = SpectralClustering(
n_clusters=2, affinity="nearest_neighbors", n_neighbors=k, random_state=42
)
labels_k = sc_k.fit_predict(X)
n_comp = connected_components(sc_k.affinity_matrix_)[0]
ari_k = adjusted_rand_score(y_true, labels_k)
print(f"n_neighbors={k:4d} connected_components={n_comp:3d} ARI={ari_k:.3f}")n_neighbors= 2 connected_components= 97 ARI=0.001
n_neighbors= 5 connected_components= 3 ARI=0.716
n_neighbors= 10 connected_components= 2 ARI=1.000
n_neighbors= 15 connected_components= 1 ARI=1.000
n_neighbors= 30 connected_components= 1 ARI=0.716
n_neighbors= 100 connected_components= 1 ARI=0.443The overall pattern is exactly what Figure 1.10 showed: too few neighbours fragments the graph into many pieces and clustering fails, a moderate neighbourhood works well, and too many neighbours blends the two crescents together and clustering fails again, just for the opposite reason.
Look closely at n_neighbors = 2 and n_neighbors = 5. The connected-component counts differ between R and Python, and at n_neighbors = 5 so does the ARI (1.000 in R, 0.716 in Python). Both languages implement the same rule, “connect each point to its 10 nearest neighbours, including itself, then symmetrize,” but R’s order()-based neighbour search and scikit-learn’s tree-based search do not necessarily break ties the same way, so the exact set of edges near the fragmentation boundary can differ. Once the graph is safely past that boundary (n_neighbors = 10 and above here), both languages agree closely. Near a fragmentation or over-connection boundary, do not expect bit-identical behaviour across implementations, and do not trust a single neighbourhood value without checking a few nearby ones.
2.5 Example B: a similarity matrix you already have
Sometimes you do not start from feature data at all. A similarity matrix might come from a domain-specific score, a kernel computed elsewhere, or, as in the next section, another pipeline entirely. Both packages accept this directly.
W10 <- knn_affinity(X, k = 10)
labels_precomputed <- SNFtool::spectralClustering(W10, K = 2, type = 3)
ari(labels_precomputed, y_true)[1] 1import numpy as np
from sklearn.neighbors import kneighbors_graph
C10 = kneighbors_graph(X, n_neighbors=10, include_self=True)
W10 = np.asarray((0.5 * (C10 + C10.T)).todense())
sc_precomputed = SpectralClustering(n_clusters=2, affinity="precomputed", random_state=42)
labels_precomputed = sc_precomputed.fit_predict(W10)
adjusted_rand_score(y_true, labels_precomputed)1.0SNFtool::spectralClustering(affinity, K, type) takes the similarity matrix as its first argument, no feature data at all. type = 3 selects the symmetric normalized Laplacian, the convention this course uses throughout. We return to the three types in the normalization section of the mathematics chapter. SpectralClustering(affinity="precomputed") does the same thing in Python: whatever matrix you hand it is used directly as the similarity graph, with no further processing.
Both W10 matrices above have entries of 1 on the diagonal, because include_self=True (Python) and the matching R construction count each point as one of its own nearest neighbours. That self-similarity feeds directly into each observation’s degree, and both packages compute the degree simply as a row sum of whatever matrix you give them. If your own similarity matrix has a different convention on the diagonal, zero, or some other constant, check that it means what you think it means, since it is not discarded anywhere in the pipeline.
2.6 Example C: an SNF fused similarity matrix
Similarity Network Fusion produces exactly the kind of object spectral clustering is built to consume: a similarity matrix, already integrating several data types into one patient network. Here we reuse the 60-patient fused network from the companion SNF course, data/snf-fused.csv, built from simulated RNA and methylation data with three known patient subtypes.
fused <- as.matrix(read.csv("data/snf-fused.csv", row.names = 1, check.names = FALSE))
patients <- read.csv("data/snf-simulated-patients.csv")
snf_labels <- SNFtool::spectralClustering(fused, K = 3, type = 3)
table(Inferred = snf_labels, Truth = patients$TrueSubtype)
ari(snf_labels, patients$TrueSubtype) Truth
Inferred 1 2 3
1 0 1 19
2 20 0 0
3 0 19 1
[1] 0.9017fused = pd.read_csv("data/snf-fused.csv", index_col=0)
patients = pd.read_csv("data/snf-simulated-patients.csv")
sc_snf = SpectralClustering(n_clusters=3, affinity="precomputed", random_state=1)
snf_labels = sc_snf.fit_predict(fused.to_numpy())
print(pd.crosstab(snf_labels, patients["TrueSubtype"],
rownames=["Inferred"], colnames=["Truth"]))
adjusted_rand_score(patients["TrueSubtype"], snf_labels)Truth 1 2 3
Inferred
0 20 0 0
1 0 0 19
2 0 20 1
0.9496Both recover the three known subtypes almost perfectly, with a small, expected difference in ARI (0.90 in R, 0.95 in Python) that traces back to the same source documented in the SNF course: SNFtool and the packages used to build the fused matrix do not normalize identically at every step, so the two fused matrices being clustered are close but not numerically identical. The clustering itself, given a fused matrix, is otherwise exactly the same operation as Example B, just with K = 3 instead of K = 2.
The companion SNF course also clusters this fused network with PAM (partitioning around medoids). PAM expects a dissimilarity: small values mean similar, large values mean different. A fused SNF network is the opposite convention, a similarity, where large values mean similar. Using it with PAM requires an explicit conversion, for instance \(d_{ij} = 1 - W_{ij}\), chosen so that it makes sense for the similarity at hand. Spectral clustering needs no such conversion, because it was built from the beginning to work directly from a similarity, or affinity, matrix.
2.7 Looking inside the algorithm: the spectral embedding
Every example so far only asked for final cluster labels. But the labels are the last, easy step. The interesting object spectral clustering actually produces along the way is the embedding itself, the new coordinates from Figure 1.7. Neither specc nor spectralClustering hands this back to you directly, so here we build it by hand, using exactly the construction from the mathematics chapter: a normalized Laplacian, its smallest eigenvectors, rescaled by \(D^{-1/2}\).
W <- knn_affinity(X, k = 10)
degree <- rowSums(W)
Dinv_sqrt <- diag(1 / sqrt(degree))
n <- nrow(W)
L_sym <- diag(n) - Dinv_sqrt %*% W %*% Dinv_sqrt
eig <- eigen(L_sym, symmetric = TRUE)
ord <- order(eig$values)
eigenvalues <- eig$values[ord]
eigenvectors <- eig$vectors[, ord]
round(eigenvalues[1:5], 5)
# After D^{-1/2} rescaling, the first direction is constant; skip it.
embedding <- Dinv_sqrt %*% eigenvectors[, 2:3]
dim(embedding)
round(embedding[1:3, ], 4)[1] 0.00000 0.00038 0.00255 0.00371 0.01105
[1] 300 2
[,1] [,2]
[1,] 0.0172 -0.0002
[2,] 0.0207 0.0053
[3,] 0.0100 -0.0067degree = W.sum(axis=1)
D_inv_sqrt = np.diag(1 / np.sqrt(degree))
n = W.shape[0]
L_sym = np.eye(n) - D_inv_sqrt @ W @ D_inv_sqrt
eigenvalues, eigenvectors = np.linalg.eigh(L_sym)
print(np.round(eigenvalues[:5], 5))
# After D^{-1/2} rescaling, the first direction is constant; skip it.
embedding = D_inv_sqrt @ eigenvectors[:, 1:3]
print(embedding.shape)
print(np.round(embedding[:3], 4))[-0. 0.0006 0.00303 0.00494 0.01363]
(300, 2)
[[ 0.0164 -0.0001]
[ 0.0209 -0.0054]
[ 0.0075 0.0064]]embedding has exactly one row per observation and one column per eigenvector kept, here two. Row 1 of the R output, (0.0172, -0.0002), is observation 1’s new position: not a gene, not a measurement, just a coordinate that reflects where observation 1 sits in the similarity network relative to everyone else. The smallest eigenvalue is (numerically) zero. After multiplication by \(D^{-1/2}\), its eigenvector becomes the constant, uninformative direction described in the mathematics chapter, which is why it is skipped here and the next two directions are used instead. The second eigenvalue itself differs slightly between R (0.00038) and Python (0.0006), for the same reason the two languages’ connected-component counts differed in the neighbourhood-size table above: knn_affinity()’s order()-based tie-breaking and scikit-learn’s tree-based search do not necessarily pick identical neighbour sets.
Feeding this straight into k-means recovers the same result as before, because this is, up to numerical detail, exactly what specc and SpectralClustering do internally:
set.seed(42)
km_embed <- kmeans(embedding, centers = 2, nstart = 10)
ari(km_embed$cluster, y_true)[1] 1km_on_embedding = KMeans(n_clusters=2, n_init=10, random_state=42).fit(embedding)
adjusted_rand_score(y_true, km_on_embedding.labels_)1.02.8 Diagnostics: how many clusters?
Spectral clustering does not choose the number of clusters for you, but the eigenvalues offer a useful hint, the eigengap already previewed conceptually in the intuition chapter and derived properly in the mathematics chapter. Let’s check it on the SNF fused network, where we happen to know the true answer is three.
d_fused <- rowSums(fused)
Dinv_fused <- diag(1 / sqrt(d_fused))
n_fused <- nrow(fused)
L_fused <- diag(n_fused) - Dinv_fused %*% fused %*% Dinv_fused
vals_fused <- sort(eigen(L_fused, symmetric = TRUE, only.values = TRUE)$values)
round(vals_fused[1:8], 4)
round(diff(vals_fused[1:8]), 4)[1] 0.0000 0.1792 0.3826 0.4168 0.4269 0.4387 0.4432 0.4516
[1] 0.1792 0.2034 0.0343 0.0101 0.0118 0.0045 0.0084d_fused = fused.to_numpy().sum(axis=1)
Dinv_fused = np.diag(1 / np.sqrt(d_fused))
n_fused = fused.shape[0]
L_fused = np.eye(n_fused) - Dinv_fused @ fused.to_numpy() @ Dinv_fused
vals_fused = np.sort(np.linalg.eigvalsh(L_fused))
print(np.round(vals_fused[:8], 4))
print(np.round(np.diff(vals_fused[:8]), 4))[-0. 0.1792 0.3826 0.4168 0.4269 0.4387 0.4432 0.4516]
[0.1792 0.2034 0.0343 0.0101 0.0118 0.0045 0.0084]The largest gap sits after the second eigenvalue (0.1792, then 0.2034), which on its own would suggest two clusters, not three. This matches what the companion SNF course found using SNFtool::estimateNumberOfClustersGivenGraph(): the eigengap heuristic favours two groups here, because two of the three true subtypes are closer to each other in the fused network than either is to the third. The eigengap is a genuinely useful diagnostic, and it is not a substitute for knowing your data. Treat it as one signal among several, alongside domain knowledge and stability under resampling, not as an automatic answer.
2.9 Summary
So the shape of the workflow is basically build or obtain a similarity matrix, hand it to a spectral clustering function (directly, from raw features, or through a precomputed affinity), and read off cluster labels. The next chapter explains how this workflow works.