Spectral Clustering

Spectral clustering groups items by building a weighted graph whose edges say how similar two items are, computing a few eigenvectors of a matrix derived from that graph (a graph Laplacian or a normalized affinity matrix), and running k-means on the rows of those eigenvectors. The eigenvectors give each item new coordinates in which items that are connected through chains of similar neighbours sit together, so the method can recover groups that are not convex in the original space: rings, spirals, curved manifolds. Normalized cuts is the image-segmentation form of the same idea and owns its derivation from a graph-cut criterion; this article covers the general data-clustering algorithm, its linear algebra, its analysis by Ng, Jordan and Weiss (2001), and what scikit-learn computes. Clustering places it among the other families of methods.

Problem

Given $n$ items and a number of clusters $K$, partition the items so that items in the same group are strongly connected and items in different groups weakly connected, without assuming any shape for a group. Ng, Jordan and Weiss (2001, §1) motivate the method against generative mixture models, which need simplifying assumptions such as Gaussian clusters and have many local optima, and note that k-means has similar problems. Their Figure 1 shows k-means failing on clusters that are not convex regions, where their spectral algorithm succeeds.

Inputs and Outputs

• Input: either feature vectors $\mathbf{s}_1, \dots, \mathbf{s}_n$ and a rule for turning them into similarities, or a precomputed symmetric affinity matrix $W$ with $w_{ij} = w_{ji} \ge 0$; and the number of clusters $K$.
• Output: a label for each item. Along the way the method produces an $n \times K$ embedding (one row per item) and the leading eigenvalues, which are themselves used to choose $K$ and to diagnose the graph.

Spectral clustering labels only the items it was run on; scikit-learn calls such methods transductive (see clustering).

Intuition

Think of the similarity graph as a network of springs, or of a random walker who moves from item to item with probability proportional to the affinity. If the graph falls into $K$ pieces with no edges between them, a walker never leaves its piece, and any function that is constant on each piece has zero “roughness” over the graph. The eigenvectors with eigenvalue 0 are exactly such piecewise-constant functions, so stacking $K$ of them gives every item of one piece the same row: the pieces collapse to $K$ points that k-means separates trivially, whatever their shape was in the original space.

Real graphs have weak edges between groups. Then the smallest eigenvalues are only close to 0, and the eigenvectors are close to piecewise constant. How close depends on how much larger the next eigenvalue is, the eigengap: a large gap means the groups are well separated relative to how tightly each is connected inside.

Algorithm

The algorithm of Ng, Jordan and Weiss (2001, §2), often called NJW, for points $S = \{\mathbf{s}_1, \dots, \mathbf{s}_n\}$ and $K$ clusters:

1. Form the affinity matrix $A_{ij} = \exp(-\lVert \mathbf{s}_i - \mathbf{s}_j \rVert^2 / 2\sigma^2)$ for $i \ne j$, and $A_{ii} = 0$.
2. Let $D$ be the diagonal matrix of row sums of $A$ and form $L = D^{-1/2} A D^{-1/2}$.
3. Find the $K$ largest eigenvectors $\mathbf{x}_1, \dots, \mathbf{x}_K$ of $L$ (chosen orthogonal when eigenvalues repeat) and stack them as the columns of $X \in \mathbb{R}^{n \times K}$.
4. Renormalize each row of $X$ to unit length: $Y_{ij} = X_{ij} / (\sum_j X_{ij}^2)^{1/2}$.
5. Cluster the rows of $Y$ into $K$ clusters with k-means or another algorithm that minimizes distortion.
6. Assign item $i$ to cluster $j$ if row $i$ of $Y$ was assigned to cluster $j$.

NJW call their $L$ what spectral graph theory would write as $I - L$; they note that the change only maps each eigenvalue $\lambda$ to $1 - \lambda$ and leaves the eigenvectors unchanged (§2, footnote 1). Murphy (§21.5.2) writes the same algorithm with $L_{\text{sym}} = I - D^{-1/2} W D^{-1/2}$ and its $K$ smallest eigenvectors.

Two other variants differ in steps 2–4:

• Unnormalized: take the $K$ eigenvectors of $L = D - W$ with the smallest eigenvalues and run k-means on their rows, without row normalization (Murphy, §21.5.2).
• Random walk (Shi–Malik): take the $K$ smallest eigenvectors of $L_{\text{rw}} = I - D^{-1} W$, equivalently the generalized eigenvectors of $(D - W)\mathbf{u} = \lambda D \mathbf{u}$, and run k-means on their rows, again without row normalization (Murphy, §21.5.4.2). This is the variant scikit-learn implements (see Implementation). As the Normalized Cuts article describes, these are the eigenvectors of Shi and Malik’s (2000) relaxation, applied there recursively two ways at a time.

Mathematical Formulation

Graph Laplacians. For a symmetric weight matrix $W$ with degrees $d_i = \sum_j w_{ij}$ and $D = \operatorname{diag}(d_i)$, the (unnormalized) graph Laplacian is $L = D - W$ (Murphy, eq. 20.125). For any vector $\mathbf{f}$ on the nodes,

$$\mathbf{f}^\top L \mathbf{f} = \frac{1}{2} \sum_{i,j} w_{ij} (f_i - f_j)^2 ,$$

the Dirichlet energy (Murphy, eq. 20.129). Since $w_{ij} \ge 0$, $L$ is symmetric and positive semidefinite, with eigenvalues $0 = \lambda_1 \le \lambda_2 \le \dots \le \lambda_n$; eigenvectors with small eigenvalues vary slowly over strongly connected nodes. The two normalized forms are

$$L_{\text{sym}} = D^{-1/2} L D^{-1/2} = I - D^{-1/2} W D^{-1/2}, \qquad L_{\text{rw}} = D^{-1} L = I - D^{-1} W$$

(Murphy, eqs. 21.30 and 21.31). They have the same eigenvalues, with eigenvectors related by $\mathbf{u} = D^{-1/2}\mathbf{x}$: if $L_{\text{sym}}\mathbf{x} = \lambda \mathbf{x}$ then $L_{\text{rw}} D^{-1/2}\mathbf{x} = \lambda D^{-1/2}\mathbf{x}$, which follows by multiplying on the left by $D^{-1/2}$. This is how scikit-learn obtains $L_{\text{rw}}$ eigenvectors from a symmetric solver. $P = D^{-1} W = I - L_{\text{rw}}$ is the transition matrix of the random walk that moves from $i$ to $j$ with probability $w_{ij}/d_i$ (Murphy, §21.5.4.2).

The ideal case. The eigenvectors of $L$ with eigenvalue 0 are spanned by the indicator vectors $\mathbf{1}_{S_1}, \dots, \mathbf{1}_{S_K}$ of the graph’s $K$ connected components (Murphy, Theorem 21.5.1): if $\mathbf{f}^\top L \mathbf{f} = 0$, every term $w_{ij}(f_i - f_j)^2$ vanishes, so $\mathbf{f}$ is constant along every path. For $L_{\text{rw}}$ the zero eigenspace is again spanned by the indicators, and for $L_{\text{sym}}$ by $D^{1/2}\mathbf{1}_{S_k}$ (Murphy, §21.5.2 and §21.5.4.2). The number of zero eigenvalues therefore counts the connected components. Murphy notes that the rows of the bottom eigenvectors of $L$ are the same vectors that Laplacian eigenmaps (Belkin and Niyogi; Murphy §20.4.9) use as an embedding.

Ng, Jordan and Weiss analyse the same case for their matrix. Zeroing all affinities between clusters gives a block-diagonal $\hat{A}$, and $\hat{L} = \hat{D}^{-1/2}\hat{A}\hat{D}^{-1/2}$ is block diagonal too; each block has a strictly positive principal eigenvector with eigenvalue 1, and its next eigenvalue is strictly smaller, since all affinities within a cluster are positive (§3.1). Eigenvalue 1 is repeated $K$ times, so $X$ is determined only up to a rotation $XR$, and individual eigenvectors are unstable: they “can easily change due to small perturbations”, while the subspace they span can be considerably more stable (§3.1). After row normalization the rotation no longer matters. Proposition 1 states that if the off-diagonal blocks are zero and each cluster is connected, the rows of $\hat{Y}$ for cluster $i$ all equal the same unit vector $\mathbf{r}_i$, and $\mathbf{r}_1, \dots, \mathbf{r}_K$ are mutually orthogonal: $K$ orthogonal points on the unit sphere, one per true cluster. Their Figure 1h shows the rows of $Y$ for two circles lying at 90° to each other.

Perturbation and the eigengap. With between-cluster affinities $E = A - \hat{A}$, matrix perturbation theory indicates that the subspace of the first $K$ eigenvectors of $\hat{L}$ is stable to small changes in $\hat{L}$ if and only if the eigengap $\delta = |\lambda_K - \lambda_{K+1}|$ is large (§3.2; NJW state it for $K = 3$ with their eigenvalue ordering). In the ideal case $\lambda_K = 1$ and $\lambda_{K+1}$ is the largest second eigenvalue of any block, so a large gap means that every cluster, taken alone, is hard to split. NJW make this precise with four assumptions (§3.2):

• A1: each cluster’s second eigenvalue is at most $1 - \delta$. A1.1, a lower bound $h(S_i)^2/2 \ge \delta$ on the Cheeger constant of each cluster, implies A1 by a standard result of spectral graph theory, and A1 is equivalent to the random walk on each cluster mixing quickly.
• A2: for every pair of distinct clusters, $\sum_{j \in S_{i_1}} \sum_{k \in S_{i_2}} A_{jk}^2 / (\hat{d}_j \hat{d}_k) \le \epsilon_1$, where $\hat{d}_j$ is point $j$‘s degree within its own cluster (eq. 6): the squared affinities between two clusters, each scaled by the two points’ within-cluster degrees, sum to a small number.
• A3: for every point, the ratio of its total affinity to other clusters to its within-cluster degree is at most $\epsilon_2$ times a factor that NJW note is $O(1)$ for densely connected clusters (eq. 7): every point is far more connected to its own cluster than to the others.
• A4: every point’s within-cluster degree is at least $1/C$ times the mean degree of its cluster: no point is much less connected than the others in its cluster.

Theorem 2 sets $\epsilon = \sqrt{K(K-1)\epsilon_1 + K\epsilon_2^2}$. Under A1–A4, if $\delta > (2 + \sqrt{2})\epsilon$, there are $K$ orthogonal unit vectors $\mathbf{r}_1, \dots, \mathbf{r}_K$ such that

$$\frac{1}{n} \sum_{i=1}^{K} \sum_{j=1}^{n_i} \lVert \mathbf{y}_j^{(i)} - \mathbf{r}_i \rVert_2^2 \le 4C\left(4 + 2\sqrt{K}\right)^2 \frac{\epsilon^2}{(\delta - \sqrt{2}\epsilon)^2}$$

(§3.2, eq. 8), where $\mathbf{y}_j^{(i)}$ is the row of $Y$ for the $j$-th point of cluster $S_i$: the rows lie close, on average, to $K$ orthogonal unit vectors, one per true cluster. The eigengap is thus the quantity that decides whether the embedding separates the clusters; it is also why choosing $K$ from the eigenvalues works when clusters are well pronounced, the eigengap heuristic that Clustering describes after von Luxburg (2007).

Graph cuts. Murphy (§21.5.1) defines the $K$-way normalized cut with a factor $\tfrac{1}{2}$,

$$\operatorname{Ncut}(S_1, \dots, S_K) = \frac{1}{2} \sum_{k=1}^{K} \frac{\operatorname{cut}(S_k, \bar{S}_k)}{\operatorname{vol}(S_k)}, \qquad \operatorname{vol}(A) = \sum_{i \in A} d_i ,$$

notes that minimizing it over binary assignments is NP-hard, and presents the eigenvectors as its continuous relaxation; for $K = 2$, Ncut equals $p(\bar{S} \mid S) + p(S \mid \bar{S})$, the probabilities that the random walk in its stationary distribution, starting in one set, moves to the other in one step (eq. 21.32). Shi and Malik’s two-way criterion, written without the $\tfrac{1}{2}$, its derivation and the relation of RatioCut to the unnormalized Laplacian are in Normalized Cuts. The scikit-learn user guide (§2.3.5) says that for two clusters SpectralClustering “solves a convex relaxation of the normalized cuts problem”; “convex” is the guide’s word, and the relaxed problem it solves is the eigenproblem above.

Parameters

• The similarity graph. The graph decides what “connected” means. Murphy (§21.5.1) notes that a node is typically connected only to its most similar neighbours, which keeps the graph sparse. NJW use a fully connected Gaussian graph. scikit-learn offers an RBF graph (affinity='rbf'), a $k$-nearest-neighbour connectivity graph ('nearest_neighbors'), and precomputed affinities or neighbour distances.
• Scale ($\sigma$ or $\gamma$). The Gaussian scale sets how fast affinity falls with distance. NJW (§4) choose $\sigma^2$ automatically: they search over it and keep the value whose rows of $Y$ form the tightest clusters (smallest k-means distortion), the property Theorem 2 predicts for the right scale. scikit-learn’s RBF graph is $\exp(-\gamma \lVert \mathbf{x} - \mathbf{y} \rVert^2)$ with gamma=1.0 by default, whatever the scale of the data, so $\gamma = 1/(2\sigma^2)$ is a fixed $\sigma^2 = 0.5$ in the data’s units. SpectralEmbedding instead defaults gamma to $1/n_{\text{features}}$.
• Neighbours ($k$). In scikit-learn’s 'nearest_neighbors' graph, n_neighbors=10 by default, and SpectralClustering.fit builds it with kneighbors_graph(..., include_self=True) (cluster/_spectral.py), which queries each point against the data including itself (neighbors/_graph.py), so each point’s 10 neighbours are itself and 9 others. The 0/1 matrix is symmetrized as $\tfrac{1}{2}(A + A^\top)$, so its weights are 1 for mutual neighbours and 0.5 for one-sided ones.
• Number of clusters $K$. Fixed in advance in scikit-learn; the user guide says the method works well for a small number of clusters and is not advised for many. The eigengap can suggest $K$.
• Label assignment. NJW initialize k-means with centroids chosen to be about 90° apart and run it once (§4, footnote 3). scikit-learn’s assign_labels is 'kmeans' (with n_init=10), 'discretize', which alternates between a discrete partition and a rotation of the embedding following Yu and Shi (2003) as its docstring describes, or 'cluster_qr', which picks one row per cluster by a pivoted QR factorization and assigns each item to its largest absolute entry after a rotation, following Damle, Minden and Ying (2019) as the docstring cites. The user guide describes k-means as able to match finer details but unstable without a fixed random_state, 'discretize' as fully reproducible but tending to make parcels of even, geometrical shape, and 'cluster_qr' as a deterministic alternative.

Complexity

A dense affinity matrix has $n^2$ entries, and a full dense eigendecomposition costs $O(n^3)$; only the $K$ (or $K + 1$) extreme eigenvectors are needed, so sparse graphs and iterative eigensolvers are the practical route. scikit-learn’s eigen_solver is 'arpack' by default, with 'lobpcg' and 'amg' (which needs the pyamg package) as alternatives, and the user guide notes that the method is especially efficient when the affinity matrix is sparse and the AMG solver is used. For spectral embedding, the user guide gives the overall cost as $O[D \log(k) N \log(N)] + O[D N k^3] + O[d N^2]$ for $N$ points of input dimension $D$, $k$ neighbours and output dimension $d$; it lists the stages as graph construction, Laplacian construction and partial eigendecomposition.

Implementation

What scikit-learn 1.7.2 computes. Reading cluster/_spectral.py, manifold/_spectral_embedding.py and SciPy’s sparse/csgraph/_laplacian.py:

• SpectralClustering.fit builds affinity_matrix_ and calls _spectral_embedding(..., drop_first=False). A code comment there describes the result as the relaxed Ncut solution, obtained by solving $L_{\text{sym}}\mathbf{x} = \lambda\mathbf{x}$ and recovering $\mathbf{u} = D^{-1/2}\mathbf{x}$.
• _spectral_embedding calls SciPy’s csgraph.laplacian(adjacency, normed=True, return_diag=True), sets the diagonal of the result to 1, and asks ARPACK in shift-invert mode for the eigenvalues of $-L_{\text{sym}}$ nearest 1.0, with tolerance 0 when eigen_tol='auto'. It divides the eigenvectors by the returned square roots of the degrees (embedding / dd) and flips each vector’s sign so that its largest-magnitude entry is positive (_deterministic_vector_sign_flip). The first vector is kept (drop_first=False), because it is constant only for connected graphs, as the comment says.
• k-means then runs on these rows with no unit-length normalization. scikit-learn thus clusters $L_{\text{rw}}$ eigenvectors, the random-walk variant, although its references include NJW.
• The diagonal. The RBF affinity_matrix_ has ones on its diagonal ($\exp(0)$), and the kNN graph has the point as its own neighbour. Neither reaches the Laplacian: SciPy’s docstring says that diagonal entries of the adjacency matrix “are ignored and replaced with zeros” for normalization, _laplacian_dense zeroes them (np.fill_diagonal(m, 0)) before computing degrees, and _laplacian_sparse, used for the sparse kNN graph, subtracts them from the degrees. The Laplacian scikit-learn decomposes therefore has NJW’s $A_{ii} = 0$. Check C6 below confirms it numerically.
• _spectral_embedding warns “Graph is not fully connected, spectral embedding may not work as expected” when the graph has more than one component.

Example: grating patches. The example clusters image patches whose groups are rings. Each patch is a $15 \times 15$ grey sinusoidal grating, $I = 0.5 + 0.4 \sin(2\pi f (x\cos t + y \sin t) + \varphi)$ with $f = 0.2$ cycles per pixel, at one of three orientations $t \in \{0°, 60°, 120°\}$, with phase $\varphi$ drawn uniformly from $[0, 2\pi)$ and Gaussian pixel noise of standard deviation 0.05; 200 patches per orientation, 600 in all, and orientation is the label to recover. Writing $0.4\sin(u + \varphi) = 0.4(\cos\varphi \sin u + \sin\varphi \cos u)$ shows the construction: in the 225-dimensional pixel space, the patches of one orientation lie on a circle of squared radius $R^2 = 0.4^2 \cdot 225/2 = 18$ around the grey patch $0.5 \cdot \mathbf{1}$, in the plane spanned by the sine and cosine gratings of that orientation, and the three planes are nearly orthogonal. All three rings share one centre. A centroid method has nothing to work with: the mean of each orientation is close to the mean of all patches. Neighbouring phases, though, are close in pixel space, so a neighbourhood graph links each ring into a chain. Opposite phases on one ring are farther apart ($4R^2$ in squared distance, without noise) than patches on different rings ($2R^2$). Pixel coordinates are $x, y = 0, \dots, 14$, with $t$ measured from the $x$ axis.

The script compares k-means on raw pixels with three spectral variants implemented in NumPy (NJW; $L_{\text{rw}}$ as scikit-learn computes it; unnormalized $L$) and with scikit-learn’s SpectralClustering, on five graphs:

• scikit-learn’s 10-nearest-neighbour graph, re-implemented in NumPy (the point itself counted, symmetrized by averaging);
• the default RBF graph, $\gamma = 1$, on pixel values in $[0, 1]$, and on the same patches scaled to $[0, 255]$;
• a Gaussian graph with $\sigma$ set to the median pairwise distance, on both scales.

All spectral variants use $K = 3$ and the same KMeans(3, n_init=10, random_state=0); ARI (adjusted Rand index, 1 for a perfect match, about 0 for chance) is printed as the range over data seeds 0–4. For seed 0 it prints the first ten $L_{\text{rw}}$ eigenvalues, the number of connected components, the number of eigenvalues below $10^{-10}$, where the largest gap among the ten lies, and each eigenvector’s orientation share, the fraction of its variance explained by the orientation labels (1 for a vector constant on each orientation, 0 for one that averages to the same value on all three). A last block removes phase sectors: the phase is drawn from two opposite arcs, $[0, \pi\rho)$ and $[\pi, \pi + \pi\rho)$, for $\rho = 1$ (the full circle, the same data as above), 0.5 and 0.25, on the 10-NN graph.

The outputs came from NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 on one BLAS thread; the script runs in about 6 seconds. The script was planned, with its checks and hypotheses, before it first ran:

• C1 (data): each orientation’s noise-free patches span exactly two dimensions, its two gratings have equal norm and are orthogonal (so each ring is a circle), gratings of different orientations have $|\cos| < 0.05$, every orientation’s mean lies within $0.15R$ of the overall mean, and $x = R^2/\sigma^2$ for the median $\sigma$ lies below the crossover $x^*$ defined in H3c.
• C2: the NumPy $L_{\text{rw}}$ embedding spans the same 9-dimensional subspace as scikit-learn’s spectral_embedding, on the kNN and $\gamma = 1$ graphs (largest sine of the principal angles below $10^{-6}$). Subspaces, not vectors, are compared because eigenvectors are defined only up to sign, and up to rotation where eigenvalues repeat.
• C3: the NumPy kNN graph equals scikit-learn’s affinity_matrix_ off the diagonal.
• C4: the number of eigenvalues below $10^{-10}$ equals the number of connected components, on every graph and seed (Theorem 21.5.1).
• C5: the median-$\sigma$ graph is the same on both scales (to $10^{-12}$), with identical labels.
• C6: scikit-learn’s RBF affinity_matrix_ has ones on its diagonal, yet its embedding matches NumPy with $A_{ii} = 0$ (below $10^{-6}$) and not NumPy with the ones kept in the degrees (above $10^{-3}$).
• C7: NJW’s Proposition 1 on their ideal case, the $\gamma = 1$ graph with all cross-orientation affinities set to zero: the rows of $Y$ for each orientation coincide, and the three points are orthogonal (both to $10^{-10}$).
• H1: on the kNN graph, k-means has ARI below 0.1 on every seed, and NJW, $L_{\text{rw}}$, $L$ and scikit-learn each above 0.95.
• H2a: on the kNN graph (seed 0), $L_{\text{rw}}$ has exactly three eigenvalues below $10^{-10}$ and $\lambda_4 > 10^{-3}$. H2b, from a planning argument: the largest gap among the first ten eigenvalues is not after $k = 3$ but after $k = 9$, because each ring adds a pair of nearly equal eigenvalues for each Fourier mode of phase, growing with the mode number, so three zeros are followed by six eigenvalues of the first modes.
• H3a: with $\gamma = 1$ on $[0, 1]$ data, all four spectral methods exceed ARI 0.95. H3b: with $\gamma = 1$ on $[0, 255]$ data every off-diagonal affinity underflows to 0, the NumPy Laplacians are undefined, and scikit-learn’s ARI is below 0.1. H3c, from a planning argument: with the median $\sigma$, all four methods have ARI below 0.1 on both scales. The argument models the rings as exact circles in orthogonal planes. Then within-ring affinities are $e^{-\nu - x} e^{x\cos\Delta}$ for a phase difference $\Delta$, cross-ring ones $e^{-\nu - x}$, with $x = R^2/\sigma^2$ and $\nu$ a common noise term. Averaging over phase gives the vectors that contrast orientations an eigenvalue of $W$ proportional to $I_0(x) - 1$, and each ring’s first Fourier modes ($\cos\varphi$, $\sin\varphi$) one proportional to $I_1(x)$, where $I_0$ and $I_1$ are modified Bessel functions. Below the crossover $x^*$ where $I_1(x) = I_0(x) - 1$, the six phase modes therefore have smaller $L_{\text{rw}}$ eigenvalues than the two orientation modes, and eigenvectors 2 and 3 encode phase. Its measure: on the seed-0 median-$\sigma$ graph, eigenvectors 2–7 have orientation share below 0.1 and eigenvectors 8 and 9 above 0.9. For $\gamma = 1$ on $[0, 1]$ data, $\sigma^2 = 0.5$ and $x = 36$.
• H4: with phase sectors removed, the number of connected components and of eigenvalues below $10^{-10}$ is 3 at $\rho = 1$ and 6 at $\rho = 0.5$ and 0.25, on every seed, so the eigengap moves from 3 to 6. ARI at $K = 3$ is printed without a hypothesis.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                     # one BLAS thread: reproducible sums
import warnings
import numpy as np
from scipy.linalg import eigh, subspace_angles
from scipy.sparse.csgraph import connected_components
from scipy.special import i0, i1
from scipy.optimize import brentq
from sklearn.cluster import KMeans, SpectralClustering
from sklearn.manifold import spectral_embedding
from sklearn.metrics import adjusted_rand_score as ari

warnings.filterwarnings("ignore", message="Graph is not fully connected")
P, F, K = 15, 0.2, 3                        # patch size, cycles per pixel, clusters
yy, xx = np.mgrid[0:P, 0:P].astype(float)
THETAS = np.deg2rad([0, 60, 120])

def basis(t):
    """The two unit-amplitude gratings (sine and cosine phase) at orientation t."""
    u = 2 * np.pi * F * (xx * np.cos(t) + yy * np.sin(t))
    return np.sin(u).ravel(), np.cos(u).ravel()

def make(seed, n_per=200, rho=1.0, noise=0.05):
    """n_per noisy 15x15 gratings per orientation. rho = 1: phase uniform on [0, 2 pi);
    rho < 1: phase on two opposite arcs [0, pi rho) and [pi, pi + pi rho)."""
    rng = np.random.default_rng(seed)
    X, y = [], []
    for k, t in enumerate(THETAS):
        b = rng.integers(0, 2, n_per)
        phi = np.pi * b + rng.uniform(0, np.pi * rho, n_per) if rho < 1 else rng.uniform(0, 2 * np.pi, n_per)
        s, c = basis(t)                     # 0.4 sin(u + phi) = 0.4 (cos phi sin u + sin phi cos u)
        clean = 0.5 + 0.4 * (np.cos(phi)[:, None] * s + np.sin(phi)[:, None] * c)
        X.append(clean + rng.normal(0, noise, clean.shape)); y.append(np.full(n_per, k))
    return np.vstack(X), np.concatenate(y)

def sqdist(X):
    g = (X ** 2).sum(1)
    return np.maximum(g[:, None] + g[None, :] - 2 * X @ X.T, 0)

def knn_graph(X, k=10):
    """scikit-learn's construction: each point's k nearest, itself included, 0/1,
    symmetrized as (A + A^T) / 2. The diagonal is dropped (the Laplacian ignores it)."""
    n = len(X)
    idx = np.argsort(sqdist(X), axis=1, kind="stable")[:, :k]
    A = np.zeros((n, n)); A[np.repeat(np.arange(n), k), idx.ravel()] = 1
    W = 0.5 * (A + A.T); np.fill_diagonal(W, 0)
    return W

def gauss_graph(D2, sigma2):
    W = np.exp(-D2 / (2 * sigma2)); np.fill_diagonal(W, 0)   # NJW: A_ii = 0
    return W

def median_sigma2(D2):
    return np.median(D2[np.triu_indices(len(D2), 1)])         # sigma = median pairwise distance

def lsym_eig(W):
    d = W.sum(1)
    Lsym = np.eye(len(W)) - W / np.sqrt(d)[:, None] / np.sqrt(d)[None, :]
    lam, V = eigh(Lsym)                     # ascending; same eigenvalues as L_rw
    return lam, V, d

def km(Z):
    return KMeans(K, n_init=10, random_state=0).fit_predict(Z)

def spectral_all(W):
    """Labels from NJW, L_rw and unnormalized L on one graph, plus L_rw eigenpairs."""
    lam, V, d = lsym_eig(W)
    Xn = V[:, :K]                                           # = K largest of D^-1/2 W D^-1/2
    Y = Xn / np.maximum(np.linalg.norm(Xn, axis=1, keepdims=True), 1e-300)
    U = V[:, :K] / np.sqrt(d)[:, None]                      # L_rw: u = D^-1/2 x
    _, Vu = eigh(np.diag(d) - W)
    return {"NJW": km(Y), "L_rw": km(U), "L": km(Vu[:, :K])}, lam, V / np.sqrt(d)[:, None]

def sin_max(A, B):
    return np.sin(subspace_angles(A, B)).max()

def share(u, y):
    """Fraction of u's variance explained by the orientation labels."""
    u = u - u.mean(); tot = u @ u
    return sum(len(u[y == k]) * u[y == k].mean() ** 2 for k in range(K)) / tot if tot > 1e-30 else np.nan

# ---- Data checks -----------------------------------------------------------------
print("data checks")
for k, t in enumerate(THETAS):
    s, c = basis(t)
    phi = np.random.default_rng(99).uniform(0, 2 * np.pi, 200)
    Z = 0.4 * (np.cos(phi)[:, None] * s + np.sin(phi)[:, None] * c)
    sv = np.linalg.svd(Z, compute_uv=False)
    print(f"  {np.rad2deg(t):3.0f} deg: rank-2 energy fraction {np.sum(sv[:2]**2) / np.sum(sv**2):.12f}, "
          f"|s|/|c| {np.linalg.norm(s) / np.linalg.norm(c):.4f}, cos(s, c) {abs(s @ c) / np.linalg.norm(s) / np.linalg.norm(c):.4f}")
B = [v / np.linalg.norm(v) for t in THETAS for v in basis(t)]
print(f"  max |cos| between gratings of different orientations {max(abs(B[i] @ B[j]) for i in range(6) for j in range(6) if i // 2 != j // 2):.4f}")
R2 = 0.4 ** 2 * P * P / 2
xstar = brentq(lambda x: i1(x) - (i0(x) - 1), 0.5, 10)
print(f"  R^2 = {R2:.1f}; crossover x* (I1 = I0 - 1) = {xstar:.4f}")

# NJW's ideal case: zero the cross-orientation affinities of the gamma = 1 graph (seed 0).
X, y = make(0)
Ahat = gauss_graph(sqdist(X), 0.5) * (y[:, None] == y[None, :])
lam, V, d = lsym_eig(Ahat)
Yi = V[:, :K] / np.linalg.norm(V[:, :K], axis=1, keepdims=True)
Ci = np.array([Yi[y == k].mean(0) / np.linalg.norm(Yi[y == k].mean(0)) for k in range(K)])
print(f"  ideal block graph: lambda_1-3 max {lam[:K].max():.1e}, lambda_4 {lam[K]:.3e}; NJW rows: max deviation "
      f"from their orientation's point {max(np.abs(Yi[y == k] - Yi[y == k].mean(0)).max() for k in range(K)):.1e}, "
      f"max |cos| between the three points {np.abs(Ci @ Ci.T - np.eye(K)).max():.1e}")

# ---- Five data seeds: ARI per graph and method --------------------------------------
rows = ["kNN (10, self incl.)", "gamma=1, 0-1 data", "gamma=1, 0-255 data", "median sigma, 0-1", "median sigma, 0-255"]
cols = ["NJW", "L_rw", "L", "sklearn"]
res = {r: {c: [] for c in cols} for r in rows}
kmeans_raw = {"0-1": [], "0-255": []}
checks = {"offset": [], "x": [], "C2": [], "C3": [], "C5": [], "C5lab": [], "C6d": [], "C6a0": [], "C6a1": [],
          "C4": [], "zero_255": []}
for seed in range(5):
    X, y = make(seed)
    X255 = 255 * X
    gm = X.mean(0)
    checks["offset"].append(max(np.linalg.norm(X[y == k].mean(0) - gm) for k in range(K)) / np.sqrt(R2))
    kmeans_raw["0-1"].append(ari(y, km(X))); kmeans_raw["0-255"].append(ari(y, km(X255)))
    D2, D2b = sqdist(X), sqdist(X255)
    checks["x"].append(R2 / median_sigma2(D2))
    graphs = {rows[0]: knn_graph(X), rows[1]: gauss_graph(D2, 0.5), rows[2]: gauss_graph(D2b, 0.5),
              rows[3]: gauss_graph(D2, median_sigma2(D2)), rows[4]: gauss_graph(D2b, median_sigma2(D2b))}
    checks["C5"].append(np.abs(graphs[rows[3]] - graphs[rows[4]]).max())
    sk = {rows[0]: SpectralClustering(K, affinity="nearest_neighbors", n_neighbors=10, random_state=0).fit(X),
          rows[1]: SpectralClustering(K, gamma=1.0, random_state=0).fit(X)}
    try:
        sk[rows[2]] = SpectralClustering(K, gamma=1.0, random_state=0).fit(X255)
    except Exception as e:
        sk[rows[2]] = type(e).__name__
    for r in rows[3:]:
        sk[r] = SpectralClustering(K, affinity="precomputed", random_state=0).fit(graphs[r])
    labs = {}
    for r, W in graphs.items():
        if (W.sum(1) == 0).any():           # isolated nodes: D^-1/2 undefined
            checks["zero_255"].append((int((W > 0).sum()), int((W.sum(1) == 0).sum())))
            for c in cols[:3]:
                res[r][c].append(np.nan)
        else:
            lab, lam, U = spectral_all(W)
            labs[r] = lab
            for c in cols[:3]:
                res[r][c].append(ari(y, lab[c]))
            cc = connected_components(W > 0, directed=False)[0]
            checks["C4"].append(int((lam < 1e-10).sum()) == cc)
        res[r]["sklearn"].append(ari(y, sk[r].labels_) if not isinstance(sk[r], str) else np.nan)
    checks["C5lab"].append(all((labs[rows[3]][c] == labs[rows[4]][c]).all() for c in cols[:3]))
    Ask = sk[rows[0]].affinity_matrix_.toarray(); np.fill_diagonal(Ask, 0)
    checks["C3"].append(np.abs(Ask - graphs[rows[0]]).max())
    for r in rows[:2]:
        _, V, d = lsym_eig(graphs[r])
        U = V / np.sqrt(d)[:, None]
        E = spectral_embedding(graphs[r], n_components=9, drop_first=False, random_state=0)
        checks["C2"].append(sin_max(U[:, :9], E))
    Ar = sk[rows[1]].affinity_matrix_
    checks["C6d"].append(bool((Ar.diagonal() == 1).all()))
    E = spectral_embedding(Ar, n_components=9, drop_first=False, random_state=0)
    _, V, d = lsym_eig(graphs[rows[1]]); checks["C6a0"].append(sin_max(V[:, :9] / np.sqrt(d)[:, None], E))
    d1 = Ar.sum(1); _, V1 = eigh(np.eye(len(Ar)) - Ar / np.sqrt(d1)[:, None] / np.sqrt(d1)[None, :])
    checks["C6a1"].append(sin_max(V1[:, :9] / np.sqrt(d1)[:, None], E))

print("checks over seeds 0-4")
print(f"  max ||orientation mean - overall mean|| / R: {max(checks['offset']):.4f}")
print(f"  x = R^2 / sigma_median^2: {min(checks['x']):.4f} to {max(checks['x']):.4f}")
print(f"  kNN affinity, NumPy vs scikit-learn, max |diff|: {max(checks['C3'])}")
print(f"  NumPy L_rw vs spectral_embedding, 9-dim subspaces, max sin angle: {max(checks['C2']):.1e}")
print(f"  rbf affinity_matrix_ diagonal all 1: {all(checks['C6d'])}; embedding vs NumPy with A_ii = 0: "
      f"{max(checks['C6a0']):.1e}, with A_ii = 1 in the degrees: {min(checks['C6a1']):.1e}")
print(f"  median sigma, 0-1 vs 0-255 graphs, max |W diff|: {max(checks['C5']):.1e}; identical labels: {all(checks['C5lab'])}")
print(f"  gamma=1 on 0-255 data, (nonzero affinities, isolated nodes): {sorted(set(checks['zero_255']))}")
print(f"  eigenvalues < 1e-10 equal connected components, every graph and seed: {all(checks['C4'])}")
print()
print(f"ARI over seeds 0-4 (min-max); k-means on raw pixels: 0-1 data {min(kmeans_raw['0-1']):.3f}-{max(kmeans_raw['0-1']):.3f}, "
      f"0-255 data {min(kmeans_raw['0-255']):.3f}-{max(kmeans_raw['0-255']):.3f}")
print(f"{'graph':22s}" + "".join(f"{c:>15s}" for c in cols))
for r in rows:
    cells = ["n/a" if np.isnan(res[r][c]).all() else f"{np.nanmin(res[r][c]):.3f}-{np.nanmax(res[r][c]):.3f}" for c in cols]
    print(f"{r:22s}" + "".join(f"{v:>15s}" for v in cells))

# ---- Seed 0: spectrum and what each eigenvector encodes -----------------------------
X, y = make(0)
D2 = sqdist(X)
print()
print("seed 0, first ten L_rw eigenvalues; largest gap after k; orientation share of eigenvectors 1-10")
for name, W in (("kNN", knn_graph(X)), ("gamma=1", gauss_graph(D2, 0.5)), ("median sigma", gauss_graph(D2, median_sigma2(D2)))):
    lam, V, d = lsym_eig(W)
    U = V / np.sqrt(d)[:, None]
    cc = connected_components(W > 0, directed=False)[0]
    print(f"  {name}: {cc} connected components, {(lam < 1e-10).sum()} eigenvalues < 1e-10, "
          f"largest gap after k = {np.argmax(np.diff(lam[:10])) + 1}")
    print("    lambda " + " ".join(f"{v:.3e}" for v in lam[:10]))
    print("    share  " + " ".join(f"{share(U[:, j], y):9.3f}" for j in range(10)))

# ---- Phase thinning on the kNN graph ---------------------------------------------
print()
print("phase on two opposite arcs of width pi*rho, kNN graph, seeds 0-4")
for rho in (1.0, 0.5, 0.25):
    comp, zero, gap, a = [], [], [], []
    for seed in range(5):
        X, y = make(seed, rho=rho)
        W = knn_graph(X)
        lab, lam, _ = spectral_all(W)
        comp.append(connected_components(W > 0, directed=False)[0]); zero.append(int((lam < 1e-10).sum()))
        gap.append(int(np.argmax(np.diff(lam[:10])) + 1)); a.append(ari(y, lab["L_rw"]))
    print(f"  rho={rho:<4}: components {comp}, eigenvalues < 1e-10 {zero}, largest gap after k {gap}, "
          f"ARI L_rw (K=3) {' '.join(f'{v:.3f}' for v in a)}")

Output:

data checks
    0 deg: rank-2 energy fraction 1.000000000000, |s|/|c| 1.0000, cos(s, c) 0.0000
   60 deg: rank-2 energy fraction 1.000000000000, |s|/|c| 1.0000, cos(s, c) 0.0000
  120 deg: rank-2 energy fraction 1.000000000000, |s|/|c| 1.0000, cos(s, c) 0.0000
  max |cos| between gratings of different orientations 0.0387
  R^2 = 18.0; crossover x* (I1 = I0 - 1) = 3.1719
  ideal block graph: lambda_1-3 max 6.7e-16, lambda_4 8.546e-03; NJW rows: max deviation from their orientation's point 0.0e+00, max |cos| between the three points 0.0e+00
checks over seeds 0-4
  max ||orientation mean - overall mean|| / R: 0.1151
  x = R^2 / sigma_median^2: 0.4844 to 0.4851
  kNN affinity, NumPy vs scikit-learn, max |diff|: 0.0
  NumPy L_rw vs spectral_embedding, 9-dim subspaces, max sin angle: 1.0e+00
  rbf affinity_matrix_ diagonal all 1: True; embedding vs NumPy with A_ii = 0: 8.2e-14, with A_ii = 1 in the degrees: 9.2e-03
  median sigma, 0-1 vs 0-255 graphs, max |W diff|: 5.4e-15; identical labels: True
  gamma=1 on 0-255 data, (nonzero affinities, isolated nodes): [(0, 600)]
  eigenvalues < 1e-10 equal connected components, every graph and seed: False

ARI over seeds 0-4 (min-max); k-means on raw pixels: 0-1 data 0.011-0.113, 0-255 data 0.011-0.113
graph                             NJW           L_rw              L        sklearn
kNN (10, self incl.)      1.000-1.000    1.000-1.000    1.000-1.000    0.619-1.000
gamma=1, 0-1 data         1.000-1.000    1.000-1.000    1.000-1.000    1.000-1.000
gamma=1, 0-255 data               n/a            n/a            n/a   -0.002-0.001
median sigma, 0-1         0.129-0.204    0.126-0.201    0.070-0.174    0.126-0.195
median sigma, 0-255       0.129-0.204    0.126-0.201    0.070-0.174    0.126-0.195

seed 0, first ten L_rw eigenvalues; largest gap after k; orientation share of eigenvectors 1-10
  kNN: 3 connected components, 3 eigenvalues < 1e-10, largest gap after k = 9
    lambda 2.220e-16 6.661e-16 6.661e-16 2.768e-03 3.336e-03 3.952e-03 4.054e-03 4.147e-03 4.200e-03 1.220e-02
    share      1.000     1.000     1.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000
  gamma=1: 1 connected components, 3 eigenvalues < 1e-10, largest gap after k = 9
    lambda 3.026e-17 1.823e-14 2.093e-14 8.546e-03 8.767e-03 9.027e-03 1.243e-02 1.387e-02 1.752e-02 3.380e-02
    share      1.000     1.000     1.000     0.007     0.003     0.019     0.002     0.000     0.004     0.001
  median sigma: 1 connected components, 1 eigenvalues < 1e-10, largest gap after k = 1
    lambda -1.332e-15 9.117e-01 9.132e-01 9.191e-01 9.220e-01 9.289e-01 9.298e-01 9.831e-01 9.832e-01 9.920e-01
    share        nan     0.011     0.005     0.011     0.001     0.000     0.002     0.969     0.980     0.000

phase on two opposite arcs of width pi*rho, kNN graph, seeds 0-4
  rho=1.0 : components [3, 3, 3, 3, 3], eigenvalues < 1e-10 [3, 3, 3, 3, 3], largest gap after k [9, 9, 9, 9, 9], ARI L_rw (K=3) 1.000 1.000 1.000 1.000 1.000
  rho=0.5 : components [6, 6, 6, 6, 7], eigenvalues < 1e-10 [6, 6, 6, 6, 7], largest gap after k [6, 6, 6, 6, 7], ARI L_rw (K=3) 0.181 0.201 0.195 0.231 0.124
  rho=0.25: components [6, 6, 6, 6, 6], eigenvalues < 1e-10 [6, 6, 6, 6, 6], largest gap after k [6, 6, 6, 6, 6], ARI L_rw (K=3) 0.181 0.201 0.195 0.206 0.211

C1 held. Each ring is an exact circle in its own plane: rank-2 energy fraction 1.000000000000, equal norms and orthogonal gratings. The planes of different orientations are close to orthogonal (largest $|\cos|$ 0.0387), every orientation’s mean lies within $0.1151R$ of the overall mean, and the median-distance scale gives $x$ = 0.4844 to 0.4851, far below $x^* = 3.1719$.

C3, C5, C6 and C7 held. The NumPy kNN graph is scikit-learn’s exactly. The median-$\sigma$ graph is the same on both scales to $5.4 \times 10^{-15}$ with identical labels, as it must be: scaling the data by 255 scales every distance and the median by 255. scikit-learn’s RBF affinity has ones on its diagonal, yet its embedding matches NumPy with $A_{ii} = 0$ to $8.2 \times 10^{-14}$, against $9.2 \times 10^{-3}$ when the ones are kept in the degrees, so the self-affinities do not reach the Laplacian, as SciPy’s source says. In NJW’s ideal case the rows of $Y$ for each orientation coincide exactly and the three points are exactly orthogonal, Proposition 1 to the printed precision; the ideal graph’s $\lambda_1$ to $\lambda_3$ are below $10^{-15}$, with $\lambda_4 = 8.546 \times 10^{-3}$.

C2 failed. The largest sine of the principal angles between the NumPy and scikit-learn 9-dimensional embeddings is 1.0, against a tolerance of $10^{-6}$: on at least one of the two graphs and five seeds, the two embeddings do not span the same subspace. The printed value is the maximum over both graphs and all seeds, so the output does not show which graph or seed failed.

C4 failed. Equality of near-zero eigenvalues and connected components did not hold on every graph and seed. The seed-0 lines show one failing case: the $\gamma = 1$ graph has 1 connected component and 3 eigenvalues below $10^{-10}$ ($3.0 \times 10^{-17}$, $1.8 \times 10^{-14}$, $2.1 \times 10^{-14}$). In our reading, this follows from the construction: two patches on different rings are at a squared distance of about $2R^2$ plus the noise term $2 \cdot 225 \cdot 0.05^2 = 1.125$, about 37, so $\gamma = 1$ gives them an affinity near $e^{-37}$, positive in floating point and so enough to connect the graph, but so small that the graph is numerically three blocks. The theorem counts exact zeros, and a threshold on computed eigenvalues cannot tell a missing edge from a negligible one.

H1 failed in its k-means part and its scikit-learn part, and held for the three NumPy spectral methods. k-means reaches ARI 0.011 to 0.113, above the 0.1 bar on at least one seed, though still far from the true grouping; in our reading, because the three orientation means nearly coincide (C1), the true partition has nearly the largest within-cluster scatter of any partition, so k-means prefers groupings that cut across the rings. NJW, $L_{\text{rw}}$ and $L$ score 1.000 on every seed, which held. On every seed the kNN graph has exactly three connected components (the $\rho = 1$ row of the last block, which uses the same data), one per orientation by the seed-0 orientation shares of 1.000, so the zero eigenspace is spanned by the orientation indicators and every variant recovers them. scikit-learn’s SpectralClustering on the same graph (C3) scores 0.619 to 1.000, which failed: on at least one seed it does not recover the orientations. Since its graph equals the NumPy one, the difference arises after the graph is built; the output does not show where.

H2a held: three eigenvalues below $10^{-10}$, then $\lambda_4 = 2.768 \times 10^{-3}$. H2b held: after the three zeros come six eigenvalues between $2.768 \times 10^{-3}$ and $4.200 \times 10^{-3}$, then $\lambda_{10} = 1.220 \times 10^{-2}$, so the largest gap among the first ten lies after $k = 9$, on every seed (the $\rho = 1$ row). The six eigenvectors after the zeros have orientation share 0.000: they do not separate the orientations (in our reading, they are the phase modes of the planning argument). Read naively, the eigengap heuristic would return 9, although the graph has three components and the zero eigenvalues say so exactly. Its answer depends on how many eigenvalues are inspected and on whether one looks for exact zeros or for the largest jump.

H3a and H3b held. With $\gamma = 1$ on $[0, 1]$ data every method scores 1.000. The spectrum is that of three nearly disconnected rings: three eigenvalues below $10^{-13}$, then eigenvectors with orientation share near 0. On $[0, 255]$ data, the same $\gamma$ gives 0 nonzero off-diagonal affinities and 600 isolated nodes. In our reading from the construction, squared distances between patches there are of the order of $255^2$ times the noise term of 1.125, some 73,000, and $e^{-x}$ is zero in double precision for $x$ above about 745. The NumPy Laplacians are then undefined, scikit-learn returns labels without an error, and its ARI is $-0.002$ to $0.001$, chance. scikit-learn’s default $\gamma$ is a length scale in the units of the data, not a property of the method.

H3c failed on its ARI bar, while the measure of its mechanism held. With the median $\sigma$, the four methods reach ARI 0.070 to 0.204, identical on both scales (C5): far from the true grouping, but above the 0.1 bar. The orientation shares bear out the planning argument: eigenvectors 2 to 7 have shares of 0.000 to 0.011, so they carry almost no orientation information (the phase modes of the argument, in our reading), and eigenvectors 8 and 9 have 0.969 and 0.980, the orientation modes, which come too late for $K = 3$. In our reading from the construction, two thirds of all pairs lie on different rings, at a squared distance of about $2R^2$ plus noise, so the median sets $\sigma^2$ near that value ($R^2/x \approx 37$ from the printed $x$); a patch’s links to the other rings (about $e^{-1/2}$) are then stronger than its links to the far side of its own ring, at about $4R^2$ (about $e^{-1}$).

H4 failed. At $\rho = 1$ every seed has 3 components and 3 zero eigenvalues, and at $\rho = 0.25$ every seed has 6, which held; at $\rho = 0.5$ seeds 0–3 have 6 and seed 4 has 7, which failed the “6 on every seed” part. In every case the count of zero eigenvalues equals the count of components. At $\rho = 0.5$ and 0.25 the largest gap among the first ten lies right after them (6, or 7 on seed 4), so the eigengap moves with the components; at $\rho = 1$ it lies after $k = 9$ (H2b). With $K = 3$ fixed, $L_{\text{rw}}$’s ARI falls from 1.000 to between 0.124 and 0.231. In our reading, once the zero eigenspace has more dimensions than $K$, its first three eigenvectors are an arbitrary basis of part of it, so k-means groups arcs without regard to which ring they belong to.

The two runs gave identical output (C8 held).

Process. The set-up and hypotheses H1–H4 were saved before any script ran. A planning script then ran the same data, graphs and seeds and printed the checked quantities, to set the tolerances of C1–C7; its output also carried scikit-learn’s warning that a graph was not fully connected. Before running, that script was edited once, to move C7 from the kNN graph to NJW’s ideal block graph. It printed a sine of 1.00 between the NumPy and scikit-learn 3-dimensional kNN embeddings for seed 0, and the plan records that “C2 is therefore planned on 9-dimensional subspaces only”; H2a and H4 were kept as written. The orientation-share measure of H3c’s mechanism was added after that planning script had printed $x$ and $x^*$, which the saved hypotheses had estimated as about 0.5 and 3.1. The example departs from the batch plan’s description in seven ways:

1. Phase thinning is done with two opposite arcs, so that the gaps have a size known by construction.
2. The kNN graph is scikit-learn’s own construction.
3. $\gamma = 1$ and the median $\sigma$ are each run on both scales.
4. Every method runs on every graph.
5. The orientation share and checks C1, C6 and C7 were added.
6. The scikit-learn embedding is compared on 9-dimensional subspaces.
7. ARI is reported over five seeds.

The script shown was run unchanged after its pilot.

Properties and Behavior

• Shape-free clusters. A cluster only has to be connected in the graph. In the example, three rings with one common centre are separated exactly while k-means is near chance.
• Exact in the ideal case, approximate otherwise. With $K$ components, the zero eigenspace is spanned by indicators (Theorem 21.5.1), and NJW’s rows collapse to $K$ orthogonal points (Proposition 1). With weak links between groups, NJW’s Theorem 2 bounds how far the rows move in terms of the eigengap.
• Subspaces, not vectors. Individual eigenvectors of a repeated or nearly repeated eigenvalue are arbitrary up to rotation. The subspace can be considerably more stable (NJW, §3.1), which is why implementations that differ only in their eigensolver can return different eigenvectors.
• Connected components count as clusters. Every component adds a zero eigenvalue. A graph with more components than $K$, as in the example’s thinned rings, gives an embedding that cannot separate the intended groups. In Clustering’s stippled-shapes example, scikit-learn’s 10-NN graph had five components for six shapes, with two bars in one component.
• Relation to kernel PCA. NJW (§5) note similarities with kernel PCA (see kernel methods) with a Gaussian kernel. The main difference in the first steps is the normalization of $A$ into $L$ and of the rows of $X$, which they say improves performance. Murphy (§21.5.4.1) adds that kernel PCA uses the largest eigenvectors of $W$, while spectral clustering uses the smallest eigenvectors of $D - W$, and that spectral clustering tends to give better results in practice.

Limitations

• The graph is the model. The result depends on the affinity, its scale and its connectivity parameter, and the example shows each failure: a $\gamma$ wrong for the data’s units isolates every node, a scale set by the median distance ranks phase above orientation, and a neighbour graph with too many components hands k-means an arbitrary basis.
• $K$ must be chosen. scikit-learn needs $K$ in advance. The eigengap heuristic is ambiguous when clusters are not well pronounced (von Luxburg, as Clustering reports), and in the example the largest gap among the first ten eigenvalues pointed to 9, not 3.
• Cost. Dense affinities take $O(n^2)$ memory, and the eigenproblem limits the method to moderate $n$ unless the graph is sparse. The scikit-learn user guide advises against it for many clusters.
• Transductive. New items need the graph and eigenvectors recomputed, or an extension the method itself does not provide.
• Label assignment adds randomness. k-means on the embedding depends on its initialization; the user guide describes 'discretize' as fully reproducible and 'cluster_qr' as deterministic.

Variants

• Meilă and Shi’s random-walk algorithm. NJW (§4) compare their method with it. As they describe it, it normalizes $A$‘s rows to sum to 1, uses that matrix’s eigenvectors and does not renormalize the rows of $X$. NJW suggest that it may cluster badly when the clusters’ degrees vary greatly, and their implementation of it used 2,000 k-means restarts. scikit-learn’s $L_{\text{rw}}$ embedding has the same eigenvectors.
• Kannan, Vempala and Vetta’s algorithm, which identifies clusters with individual singular vectors, very frequently gave poor results in NJW’s experiments (§5, Figure 1l).
• Recursive two-way cuts. Shi and Malik’s normalized cuts split the graph with one eigenvector at a time and recurse; see Normalized Cuts, which also covers the image versions with pixel affinities, multiscale graphs and spectral boundary cues.
• Laplacian eigenmaps. The same eigenvectors used as a low-dimensional embedding rather than for clustering (Belkin and Niyogi; Murphy, §20.4.9). scikit-learn’s SpectralEmbedding implements it, with a kNN graph by default and the first, constant vector dropped.
• Other label assignments. Yu and Shi’s discretization and Damle, Minden and Ying’s QR-based assignment replace k-means in scikit-learn (Parameters).
• Learned affinities. In motion segmentation, spectral clustering runs on affinities between point trajectories, including subspace-clustering affinities, and that article reports that spectral methods, with their preference for balanced partitions, tend to miss very small segments that a multicut formulation can find. Image segmentation and superpixels use its normalized-cuts form on pixel graphs.