Spectral Clustering, in depth: graph Laplacians, the Fiedler vector, normalised cuts and k-NN graphs

The graph Laplacian

Start with a weighted, undirected graph: W[i][j] is the similarity between points i and j, and the degree d_i is the sum of row i. The graph Laplacian is L = D - W, where D is the diagonal degree matrix. Everything follows from one identity: for any vector f, f^T L f equals one half of the sum over all pairs of W[i][j] times (f_i - f_j)^2.

Three consequences follow. L is positive semidefinite, since that sum cannot be negative. The constant vector gives zero, so the smallest eigenvalue is 0. And a vector gives zero exactly when it is constant on every connected component. So the number of zero eigenvalues equals the number of connected components, and their eigenvectors are indicator vectors of the components (up to rotation). If your clusters were perfectly disconnected, the eigenvectors would simply list them.

Real clusters are joined by a few weak edges, so the eigenvalues are small but not zero, and the eigenvectors are close to piecewise constant. That is why spectral clustering takes the eigenvectors of the smallest eigenvalues of L. Some descriptions say largest, which is true only of the normalised adjacency matrix, the form the code below uses.

Worked example: two triangles and a bridge

Take two triangles, nodes 0-1-2 and 3-4-5, joined by a single bridge edge from 2 to 3, all weights 1. Nodes 2 and 3 have degree 3 and the rest have degree 2. NumPy gives the eigenvalues of L as 0, 0.438, 3, 3, 3 and 4.562. The second eigenvector, the Fiedler vector, is (-0.465, -0.465, -0.261, 0.261, 0.465, 0.465), up to sign. Its signs separate the triangles, and the bridge nodes sit closest to zero, which marks them as the least certain. Delete the bridge and the eigenvalues become 0, 0, 3, 3, 3, 3: two components, two zeros.

The second eigenvalue, 0.438, is the algebraic connectivity. It measures how close the graph is to falling apart, and the gap from it to the third eigenvalue, 3, says that two clusters is the natural answer.

Why it works: relaxing graph cuts

Clustering a graph means cutting few edges while keeping the parts balanced. The ratio cut divides the cut weight by the number of nodes on each side. The normalised cut (Shi and Malik, 2000) divides by each side's volume, the sum of its degrees. For the triangles the cut weight is 1 and each side has volume 7, so Ncut = 1/7 + 1/7 = 2/7. Minimising either is NP-hard.

The relaxation is the key step. Encode a two-way partition as a vector f that takes one value on each side. Then f^T L f is proportional to the ratio cut. Let f take any real values instead, keep it orthogonal to the constant vector, and the minimiser is the Fiedler vector of L. The same argument with the volume-weighted constraint gives the generalised problem L f = lambda D f, the eigenvectors of the random-walk Laplacian L_rw = I - D^-1 W. For k clusters, take the first k eigenvectors and round the relaxed solution back to discrete labels; k-means on the embedded rows is that rounding step.

L_rw and L_sym have the same eigenvalues, and their eigenvectors differ by a factor of D^1/2. L_sym is symmetric, which suits standard eigensolvers. Its eigenvectors are scaled by the square root of degree, so Ng, Jordan and Weiss normalise each row of the embedding to unit length before k-means. Skipping that step while using L_sym is a common silent bug.

Implementation with NumPy and SciPy

This implementation builds a symmetric k-nearest-neighbour graph with a k-d tree, computes the top eigenvectors of D^-1/2 W D^-1/2 with a sparse Lanczos solver, row-normalises them and runs k-means with restarts. The top eigenvalues of that matrix are one minus the bottom eigenvalues of L_sym, so asking for the largest algebraic eigenvalues avoids slow shift-invert mode.

import numpy as np
from scipy.sparse import csr_matrix, diags
from scipy.sparse.linalg import eigsh
from scipy.spatial import cKDTree
from scipy.cluster.vq import kmeans2

def knn_graph(X, k):
    _, idx = cKDTree(X).query(X, k + 1)       # first neighbour is the point
    n = len(X)
    rows = np.repeat(np.arange(n), k)
    A = csr_matrix((np.ones(n * k), (rows, idx[:, 1:].ravel())), shape=(n, n))
    return A.maximum(A.T)                     # edge if either side chose it

def spectral(X, n_clusters, k_nn=10, restarts=10, seed=0):
    A = knn_graph(X, k_nn)
    d = np.asarray(A.sum(1)).ravel()
    Dm = diags(1 / np.sqrt(d))
    M = Dm @ A @ Dm                           # L_sym = I - M
    v0 = np.random.default_rng(seed).standard_normal(len(X))   # reproducible
    mu, U = eigsh(M, k=n_clusters + 3, which="LA", v0=v0)
    order = np.argsort(-mu)
    lam = 1 - mu[order]                       # L_sym eigenvalues, ascending
    Uk = U[:, order[:n_clusters]]
    Uk /= np.linalg.norm(Uk, axis=1, keepdims=True)
    best = None
    for s in range(restarts):
        cent, lab = kmeans2(Uk, n_clusters, minit="++", seed=s)
        inertia = ((Uk - cent[lab]) ** 2).sum()
        if best is None or inertia < best[0]:
            best = (inertia, lab)
    return best[1], lam

On two noisy rings of 500 points each (radii 1 and 3, noise sigma 0.08), this returns 100% accuracy with k_nn = 10, where plain k-means on the coordinates scores 50.7%.

The similarity graph decides the answer

The eigenvector maths is the easy part. The graph decides the answer. Sweeping the neighbour count on the same rings:

The eigenvalues come from a dense solver on the full 1,000 by 1,000 L_sym, and the zero counts match the component counts exactly, as the theory says. Too few neighbours fragments the graph into more components than clusters. The zero eigenvalue then repeats, any basis of its eigenspace is equally valid, and the accuracy depends on which basis the solver happens to return: rerunning with five different Lanczos start vectors gave anything from coin-flip quality to 98%. Too many would eventually join the rings with short-cut edges. Here k_nn = 100 joins them into one component but still separates them, with a much smaller eigengap. Always print the component count before trusting a result.

Two other choices matter. A fully connected Gaussian kernel exp(-||x_i - x_j||^2 / 2 sigma^2) needs a well-chosen sigma and gives a dense n-by-n matrix. Self-tuning spectral clustering (Zelnik-Manor and Perona, 2004) sets a local scale per point, the distance to its seventh neighbour, which copes with clusters of different density. A mutual k-NN graph, which keeps an edge only if both points choose each other, cuts bridges through sparse noise but disconnects outliers.

To choose k, look for the eigengap: a jump after the k-th smallest eigenvalue. It works on well-separated data, as in the six-node example, and is weak on noisy data, where the jump may not exist. Treat it as a hint and check stability across resampled runs.

Scaling and out-of-sample points

A dense Laplacian on 100,000 points holds 10^10 float64 entries, about 80 GB, and a dense eigendecomposition costs O(n^3). Production pipelines avoid both. Build a sparse k-NN graph (approximate neighbour search beyond a few million points), and compute only k eigenvectors with Lanczos (eigsh) or LOBPCG with an algebraic multigrid preconditioner, which handles millions of nodes. The Nyström method samples m landmark points, decomposes the m-by-m block and extends to the rest. It is fast but sensitive to landmark choice. For very large graphs, a few steps of power iteration on the normalised adjacency (power iteration clustering) is a cheaper approximation.

Spectral clustering has no native predict step. To label new points, use the Nyström extension, or train a classifier on the clustered data and treat the clustering as labelling.

Operational guidance

Clustering has no ground truth in production, so validate structure instead of accuracy. Four checks catch most problems before anyone looks at the labels.

• Graph health. Log the number of connected components, the minimum and median degree, and the share of edges that are mutual. A jump in components after a data refresh usually means a feature changed scale.
• Spectrum. Log the first k + 3 eigenvalues on every run. A collapsed eigengap, or several eigenvalues at zero, explains most surprising clusterings.
• Stability. Rerun on bootstrap samples or different seeds and compare label agreement, for example with the adjusted Rand index. Clusters that change from run to run should not drive decisions.
• Per-cluster conductance. For each cluster, divide the weight of edges leaving it by its volume. A cluster with high conductance is barely separated from its neighbours, whatever its size.

Standardise features before building the graph, since distances decide every edge. Cache the k-NN graph, which usually dominates the run time, and keep its parameters in the run record with the eigenvalues, so a changed result can be traced to the graph or to the clustering step.

Failure modes and alternatives

Failure modes seen in practice:

• Isolated points. A zero-degree node makes D^-1/2 divide by zero. Symmetrised k-NN graphs avoid it; epsilon-ball and mutual graphs do not. Drop or attach such points first.
• More components than clusters. The k_nn = 3 and 5 rows above. Check the component count; raise k_nn or cluster per component.
• Unreliable iterative solvers on repeated eigenvalues. A Lanczos run from one start vector can miss copies of a highly repeated eigenvalue, which is exactly what a fragmented graph produces, and report a wrong spectrum. Check component counts with a graph routine, use a dense solver on small cases, or use LOBPCG with a block at least as large as the multiplicity.
• Degenerate eigenspaces. With repeated eigenvalues the solver may return any rotation of the basis, and signs are arbitrary. Never threshold individual eigenvectors for k above 2; k-means on the rows is rotation-tolerant.
• Missing row normalisation with L_sym, which pulls low-degree points towards the origin and mislabels them.
• Mixed scales. One dense and one sparse cluster defeat a single global sigma. Use local scaling or k-NN graphs.
• Noise bridges. A thin trail of points between clusters merges them. Density-based methods ignore such bridges more gracefully.

What to do next

1. Compute L and its eigenvalues for the six-node example by hand or in NumPy, and confirm 0 and 0.438.
2. Run the code on two rings, sweep k_nn, and print the component count each time.
3. Swap in L_rw by solving the generalised problem with scipy.linalg.eigh(L, D) on a small case, and compare labels with the NJW embedding.
4. Compare against k-means and DBSCAN on your own data.
5. Read eigenvector centrality for the other end of the spectrum, and Louvain for community detection on native graphs.
6. See the Stoer-Wagner minimum cut for why an unbalanced minimum cut alone is not a clustering.