What it actually does
Stop thinking about distance to a centre and start thinking about connectivity. Build a graph where nearby points are joined, then look for a way to cut it into pieces that severs as few connections as possible. Two crescents that interlock have no separating hyperplane and no useful centroids, but they are two obviously distinct connected regions, and cutting the graph finds them. The trick that makes this computable is that the best cut, once relaxed from a hard yes-or-no assignment to a continuous one, is given by the smallest eigenvectors of the graph Laplacian — so a combinatorial problem becomes a linear-algebra one. Those eigenvectors are new coordinates in which the crescents have become blobs, and then k-means finishes the job.
- Clusters are connected regions in the affinity graph — which is a much weaker assumption than being blobs, and is why the method exists
- The affinity function reflects the similarity that matters; it is the model, and a poor affinity cannot be repaired downstream
- The graph is connected, or nearly so — disconnected components take over the bottom of the spectrum
- The number of clusters is known or visible as an eigengap, which on real data is frequently ambiguous
- The dataset is small enough for an n-by-n affinity and its decomposition, which is the binding constraint in practice
Loss function
Cut the graph into K pieces so that few connections are severed, and normalize by each piece’s size so the answer is not simply to snip off one isolated vertex. That combinatorial problem is NP-hard. The second form is the relaxation actually solved: drop the requirement that H be a hard assignment, keep only the orthogonality constraint, and the optimum becomes the bottom eigenvectors of the Laplacian. The relaxation is where the honesty is owed — there is no approximation guarantee tying the relaxed solution back to the true minimum cut, and the k-means run afterwards is a heuristic on top of a heuristic. It works remarkably well and it is not solving the stated objective.
- total weight of the edges severed by separating cluster k from everything else
- total degree inside cluster k — the balance term, without which the optimum is to cut off a single point
- the graph Laplacian D - W, whose quadratic form measures exactly how much a labelling disagrees across edges
- the relaxed cluster indicator: continuous values instead of a hard assignment, which is what makes the problem solvable
How it is minimized
Eigendecomposition of the graph Laplacian, followed by k-means in the embedded space
- the affinity matrix — the actual model, and the thing that decides whether any of this works
- diagonal degree matrix; the normalization is what stops high-degree vertices dominating the embedding
- the k-th smallest eigenvector of the Laplacian — a coordinate in the new space
- the n-by-K embedding; its rows are the points, re-expressed so that connected regions become compact
There is no iteration over the objective at all: the relaxation turns the partition problem into an eigenproblem, and the eigenproblem has a closed-form answer. Recognizing that is the transferable idea — the same move appears in PCA, in kernel methods, and throughout numerical linear algebra. Three implementation decisions carry real consequences. The Laplacian normalization matters: the random-walk and symmetric forms behave far better than the unnormalized one when degrees vary, which on any realistic graph they do. Only the bottom K eigenvectors are needed, so a Lanczos method is the right tool rather than a full decomposition. And the embedding rows should be normalized to unit length before k-means, which is the step that turns elongated eigenvector geometry into the spherical shape k-means requires.
affinity construction
The model. A k-nearest-neighbour graph, a full RBF kernel, or a domain-supplied similarity — this choice matters far more than anything else here (kNN graph with 10 to 20 neighbours, symmetrized)
sigma (RBF bandwidth)
How quickly affinity decays with distance. No good default exists: too small disconnects the graph, too large connects everything and the structure vanishes (median pairwise distance, or local scaling per point)
K
Number of clusters. The eigengap heuristic — look for a jump in the sorted eigenvalues — is the one genuinely principled selection rule in this section, and it is often ambiguous in practice
Laplacian type
Unnormalized, symmetric, or random-walk. The normalized forms are almost always the right choice when vertex degrees vary
row normalization
Scale each embedding row to unit length before k-means. Small step, large effect on whether the final clustering is sensible
The eigendecomposition is exact and deterministic; the k-means afterwards is not, so the method inherits that entry’s dependence on initialization and restarts. The characteristic failures are all about the graph. If sigma is too small the graph fragments into many components, every component contributes a near-zero eigenvalue, and the bottom eigenvectors describe connectivity rather than structure — the clustering then simply recovers the fragments. If sigma is too large everything connects to everything, the eigengap closes, and the embedding is noise. A single bridging point between two genuine clusters can merge them entirely, because connectivity is transitive in a way distance is not. And the relaxation itself carries no bound: nothing guarantees the eigenvector solution is close to the true minimum cut, which is why this is judged empirically rather than argued from theory.
Affinity construction: O(n^2·d) for a dense graph, or O(n·log n·d) with an approximate neighbour index. Memory O(n^2) dense, or O(n·k) for a kNN graph — which is the difference between a few thousand points and a few hundred thousand. Eigendecomposition: O(n^3) dense, reduced to roughly O(n·k·K) per Lanczos iteration on a sparse graph. That n-squared to n-cubed profile is the ceiling, and Nyström or landmark approximations are what any large-scale spectral method is actually doing.
How, where, and why it is used
Time-Series Forecasting
AdaptedGroup series before modelling them, using a similarity graph rather than Euclidean distance. That distinction is the whole reason to prefer this over k-means here: the affinity can be built from correlation, from dynamic time warping, or from any domain similarity, none of which admits a mean and therefore none of which k-means can use.
- Grouping many series by a correlation or DTW-based affinity, so that one model can be pooled across each group
- Regime discovery over a similarity graph of windowed segments, where regimes are not separable by distance to a centre
- Sensor or asset grouping where a domain-specific similarity exists and a metric does not
It is the option that accepts a similarity you can compute but not average, which is the common situation with time series — DTW alignment being the standard example. It also finds groups that are not blobs, and series clusters frequently are not. The costs are the usual ones and they bind quickly: an n-by-n affinity over tens of thousands of series is not affordable, and there is no way to assign a new series without recomputing the embedding, so this is an offline analysis rather than a component of a serving path.
- Z-normalize each series before computing affinity, or the graph recovers magnitude rather than shape
- Build a k-nearest-neighbour graph rather than a full affinity matrix; it is both cheaper and usually better, since distant pairs contribute noise
- Fit on the training window only — an affinity computed over the full history leaks future co-movement into the grouping
- Check for disconnected components before clustering; they will otherwise consume the bottom eigenvectors and the clustering will simply recover them
Rolling-origin backtest of the pooled forecasts, not a clustering score. Inspect the eigenvalue spectrum as a diagnostic: no clear gap means the data does not support the number of groups being requested.
- A bandwidth that fragments the graph, after which the bottom eigenvectors describe connected components rather than structure
- No out-of-sample assignment — a new series cannot be placed without redoing the whole computation
- Quadratic memory reached quietly as the series count grows between refits
Anomaly Detection
Not applicableThe method partitions a graph into K pieces and produces no notion of outlyingness: a genuinely isolated point simply becomes its own component or its own cluster, which is a partition result rather than a score.
Optimization
AdaptedThis is the reference example of spectral relaxation. The normalized-cut objective is NP-hard over hard assignments; dropping the integrality constraint and keeping only orthogonality turns it into a trace minimization whose optimum is the bottom eigenvectors. The whole method is that one move, and the same move recurs in graph partitioning, in max-cut relaxations, and in embedding problems generally.
- The canonical demonstration that relaxing a combinatorial constraint can turn an NP-hard problem into an eigenproblem
- Graph partitioning for parallel workload distribution and circuit layout, where balanced cuts are the actual deliverable
- Load balancing and mesh decomposition, which are literally the normalized-cut problem under a different name
Worth studying because the gap between the relaxation and the original problem is unusually visible and unusually honest: there is no approximation guarantee, the rounding step back to a hard assignment is a heuristic, and the method is still excellent in practice. That combination — a principled relaxation, an unprincipled rounding, and strong empirical results — is the shape of a great many applied optimization methods, and this is the clearest place to see it.
- Encode the real cost structure in the edge weights; the objective optimizes the graph you supply, not the problem you meant
- Balance constraints belong in the normalization — cut alone will happily isolate a single vertex
For a partitioning problem, evaluate the actual cut and the balance achieved against a combinatorial baseline such as METIS. The relaxed objective value is not comparable to the combinatorial one and should not be reported as if it were.
- Reporting the relaxed objective as though it bounded the true optimum, which it does not
- Ignoring the rounding step, where a substantial part of the quality is won or lost
- Applying the dense formulation at a scale where only a sparse or approximate method is viable
Computer Vision
AdaptedNormalized cuts for image segmentation: treat pixels as vertices, weight edges by a combination of colour similarity and spatial proximity, and cut the resulting graph. This is the formulation that made spectral methods well known in vision, and it segments regions that no thresholding or centroid method separates.
- Classical image and video segmentation, where region boundaries follow connectivity rather than colour alone
- Superpixel and region-merging pipelines that feed a downstream recognizer
- Motion segmentation, where trajectories are grouped by an affinity that is not a metric
The affinity formulation lets colour similarity and spatial adjacency be combined into one weight, which is exactly what image segmentation needs and what a distance in pixel space cannot express. Its limits are computational and were decisive: a modest image has hundreds of thousands of pixels, and an n-by-n affinity over those is impossible without approximation. Learned segmentation superseded it for recognition, and the graph formulation survives inside superpixel and region-merging steps.
- Combine colour distance and spatial distance into one affinity with separate bandwidths, since they are in different units
- Restrict edges to a spatial neighbourhood — a fully connected pixel graph is both intractable and worse, as distant pixels contribute noise
- Work on superpixels rather than raw pixels to reduce n by two orders of magnitude before the eigendecomposition
Boundary precision-recall or region-covering measures against annotated segmentations, not the value of the cut objective.
- Quadratic memory on raw pixels, which is why superpixels are effectively mandatory
- A spatial bandwidth that merges across a genuine boundary because a few bridging pixels connect the regions
- Choosing the number of segments by the eigengap when natural images rarely produce a clean one
Risk & Fraud
ViableCommunity detection on a transaction or entity graph. Vertices are accounts, devices, or merchants; edges are shared attributes or transfers; and a tightly-connected community that is weakly connected to everything else is the structural signature of a fraud ring or a layering pattern.
- Fraud-ring and mule-network detection, where the signal is connectivity rather than any individual account’s behaviour
- Anti-money-laundering layering detection across transfer graphs
- Collusion detection in marketplaces and ad networks, where colluding parties transact mostly with each other
It sees exactly what a per-entity model cannot: coordinated behaviour that looks unremarkable account by account and obvious as a group. That is the failure mode named under the isolation forest, where a cluster of anomalies masks itself, and the graph view is the complement to it. Against it: the affinity graph must be built from real relational data, the computation is offline and quadratic, and communities are found whether or not any are fraudulent — the method proposes structure, and the investigation supplies the meaning.
- Weight edges by relationship strength — shared device, shared beneficiary, transfer volume — rather than treating all links as equal
- Prune hub vertices such as large exchanges before partitioning, or every community connects through them and the graph will not separate
- Work on a sparse relational graph, never a dense affinity over engineered features, which is both slower and less meaningful here
Precision of investigated communities against confirmed cases, and the size distribution of what is returned — a partition producing two enormous communities has told you nothing actionable.
- Hub vertices bridging every community into one, which is the most common way this fails on real financial graphs
- Treating every detected community as suspicious when most are legitimate structure
- Scale: real transaction graphs need sparse or streaming community methods rather than a dense spectral decomposition
What it costs to run
Dominated by two quadratic-or-worse steps: building the affinity and decomposing the Laplacian. A few thousand points is comfortable, tens of thousands needs a sparse graph and a Lanczos solver, and beyond that the honest answer is Nyström or a landmark approximation rather than the exact method.
There is no inference. The embedding is defined only for the points it was fitted on, so a new point cannot be assigned without recomputing — which is the property that keeps this an offline analysis. Out-of-sample extension exists (the Nyström formula) and is an approximation, not a lookup.
Batch, offline, and periodic. Because there is no incremental update and no cheap assignment, the practical pattern is to run it as an analysis, extract stable labels or centroids, and serve something else — often a k-means or a classifier fitted on the resulting labels.
- Track the eigengap across refits; a gap that closes means the structure being clustered no longer exists in the data
- Watch the number of near-zero eigenvalues, which counts connected components — a rise means the bandwidth has become too small for the current data
- Monitor the degree distribution of the affinity graph, since a shift there changes the normalization and therefore the embedding
- Compare partitions across refits with a set-overlap measure rather than by cluster id, which is meaningless between runs
- No out-of-sample assignment: a new point cannot be clustered without redoing the affinity and the decomposition, which surprises teams who expected a predict method
- The bandwidth is the model, and there is no default that works across datasets — it must be re-selected whenever the feature scale changes
- Disconnected components produce trivial eigenvectors that consume the bottom of the spectrum, so the clustering silently recovers connectivity instead of structure
- Quadratic memory is reached without warning as the dataset grows between runs; the affinity matrix, not the algorithm, is what fails
- The final k-means is stochastic, so the same eigenvectors give different labels on different runs unless the seed is pinned
Pros and cons, in context
Finds clusters that are not convex blobs
Interlocking crescents, rings, and elongated manifolds are separated cleanly, which is precisely what k-means and Gaussian mixtures cannot do. This is the entire reason to accept the cost.
Needs only a similarity, not a metric or a mean
Any affinity works — DTW, graph distance, a domain kernel — including ones with no meaningful average. That opens problems where centroid-based methods are not merely worse but undefined.
The eigenvalue spectrum is a genuine diagnostic
The eigengap is the closest thing to a principled way to choose the number of clusters anywhere in this section. It is often ambiguous on real data, and it is still more than k-means offers.
The embedding is useful on its own
Stopping before k-means gives a nonlinear embedding that preserves local structure and feeds any downstream model. One computation, two deliverables.
Quadratic memory and cubic decomposition
A hard ceiling in the low tens of thousands of points without approximation. This decides most spectral-versus-k-means arguments before cluster quality is discussed.
No way to assign a new point
The embedding exists only for the fitted points. Teams routinely discover this after building a pipeline around it, and the Nyström extension is an approximation rather than the missing method.
Results hinge on a bandwidth with no good default
Too small fragments the graph, too large erases the structure, and the usable range can be narrow. More sensitive to this one number than k-means is to anything, which is an uncomfortable trade for the extra flexibility.
The relaxation carries no guarantee
Nothing bounds the eigenvector solution against the true minimum cut, and the rounding to hard labels is a separate heuristic. It works well empirically, and any claim stronger than that is unsupported.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Spectral clustering - affinity, Laplacian, eigenvectors. Transcribed.
Build the graph, form the Laplacian, take its smallest eigenvectors, and run
k-means on the rows. The smallest eigenvectors are found here by power
iteration on a SHIFTED matrix - power iteration converges to the largest
eigenvector, so subtracting L from a multiple of the identity flips the
spectrum and turns the smallest into the largest.
"""
import math
import random
def rbf_affinity(X, sigma):
"""W_ij = exp(-||x_i - x_j||^2 / (2 sigma^2)), with no self-loops.
Sigma is the model. Too small and the graph fragments into isolated
points; too large and everything connects to everything and the structure
disappears. There is no default that works across datasets.
"""
n = len(X)
W = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(i + 1, n):
squared = 0.0
for k in range(len(X[i])):
difference = X[i][k] - X[j][k]
squared += difference * difference
weight = math.exp(-squared / (2.0 * sigma * sigma))
W[i][j] = weight
W[j][i] = weight
return W
def normalized_laplacian(W):
"""L_sym = I - D^(-1/2) W D^(-1/2).
The normalization is what stops a high-degree vertex dominating the
embedding purely because it has many neighbours.
"""
n = len(W)
degree = [sum(row) for row in W]
inverse_root = [0.0 if d == 0.0 else 1.0 / math.sqrt(d) for d in degree]
L = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(n):
L[i][j] = (1.0 if i == j else 0.0) - inverse_root[i] * W[i][j] * inverse_root[j]
return L
def power_iteration(matrix, iterations=1000, tol=1e-10):
"""Largest eigenvector of a symmetric matrix, by repeated multiplication."""
n = len(matrix)
vector = [1.0 / math.sqrt(n)] * n
for _ in range(iterations):
product = [0.0] * n
for i in range(n):
total = 0.0
for j in range(n):
total += matrix[i][j] * vector[j]
product[i] = total
norm = math.sqrt(sum(value * value for value in product))
if norm == 0.0:
break
product = [value / norm for value in product]
shift = sum(abs(product[i] - vector[i]) for i in range(n))
vector = product
if shift < tol:
break
eigenvalue = 0.0
for i in range(n):
for j in range(n):
eigenvalue += vector[i] * matrix[i][j] * vector[j]
return vector, eigenvalue
def bottom_eigenvectors(L, k):
"""The k SMALLEST eigenvectors of L, via the shifted matrix 2I - L.
Power iteration finds the largest; L_sym has eigenvalues in [0, 2], so
subtracting it from 2I reverses the order without changing the vectors.
"""
n = len(L)
shifted = [[(2.0 if i == j else 0.0) - L[i][j] for j in range(n)] for i in range(n)]
vectors = []
for _ in range(k):
vector, eigenvalue = power_iteration(shifted)
vectors.append(vector)
# Deflate so the next iteration finds the next one.
for i in range(n):
for j in range(n):
shifted[i][j] -= eigenvalue * vector[i] * vector[j]
return vectors
def fit(X, k, sigma=1.0, seed=0):
W = rbf_affinity(X, sigma)
L = normalized_laplacian(W)
vectors = bottom_eigenvectors(L, k)
# Rows of the embedding, normalized to unit length. Small step, large
# effect: it is what turns elongated eigenvector geometry into the
# spherical shape k-means requires.
n = len(X)
embedding = []
for i in range(n):
row = [vectors[c][i] for c in range(k)]
norm = math.sqrt(sum(value * value for value in row))
embedding.append([value / norm for value in row] if norm > 0 else row)
return simple_kmeans(embedding, k, seed)
def simple_kmeans(points, k, seed):
rng = random.Random(seed)
centroids = [list(points[i]) for i in rng.sample(range(len(points)), k)]
assignment = [0] * len(points)
for _ in range(100):
for i, point in enumerate(points):
best, best_distance = 0, None
for c, centroid in enumerate(centroids):
distance = sum((point[j] - centroid[j]) ** 2 for j in range(len(point)))
if best_distance is None or distance < best_distance:
best, best_distance = c, distance
assignment[i] = best
for c in range(k):
members = [points[i] for i in range(len(points)) if assignment[i] == c]
if members:
centroids[c] = [sum(m[j] for m in members) / len(members)
for j in range(len(members[0]))]
return assignmentO(n^2 d) for the affinity and O(k n^2) per power iteration, over dense lists of lists — the n-by-n matrix is the ceiling.
How this language actually wants it written.
"""Spectral clustering - typed, kNN affinity, symmetric Laplacian, eigh."""
from dataclasses import dataclass
from typing import Literal
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
Labels = NDArray[np.int64]
Affinity = Literal["knn", "rbf"]
@dataclass(frozen=True)
class SpectralEmbedding:
"""The embedding and the spectrum that produced it.
The eigenvalues are kept because they are the diagnostic: a clear gap
after position K supports K clusters, and no gap means the data does not
support the number being requested.
"""
embedding: Matrix # (n, K), rows normalized to unit length
eigenvalues: Vector # (K,)
@property
def eigengap(self) -> Vector:
"""Successive differences. The largest is the suggested K."""
return np.diff(self.eigenvalues)
def _knn_affinity(X: Matrix, n_neighbors: int, sigma: float | None) -> Matrix:
"""Symmetric k-nearest-neighbour graph with RBF weights.
Preferred to a fully-connected affinity for two reasons, not one: it is
far cheaper, and it is usually BETTER, because distant pairs contribute
weight that is noise rather than structure.
"""
squared = (
np.einsum("ij,ij->i", X, X)[:, None]
- 2.0 * (X @ X.T)
+ np.einsum("ij,ij->i", X, X)[None, :]
)
np.maximum(squared, 0.0, out=squared) # cancellation guard
np.fill_diagonal(squared, np.inf) # no self-loops
# Median heuristic: a defensible default where none really exists.
scale = sigma if sigma is not None else float(np.sqrt(np.median(squared[np.isfinite(squared)])))
neighbors = np.argpartition(squared, n_neighbors, axis=1)[:, :n_neighbors]
affinity = np.zeros_like(squared)
rows = np.repeat(np.arange(X.shape[0]), n_neighbors)
affinity[rows, neighbors.ravel()] = np.exp(
-squared[rows, neighbors.ravel()] / (2.0 * scale * scale)
)
# Symmetrize: i being a neighbour of j does not make j a neighbour of i.
return np.maximum(affinity, affinity.T)
def fit_embedding(
X: Matrix,
n_clusters: int,
n_neighbors: int = 10,
sigma: float | None = None,
) -> SpectralEmbedding:
"""Build the graph and embed. Raises ValueError on malformed input."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if not 1 <= n_clusters <= X.shape[0]:
raise ValueError(f"n_clusters must lie in [1, {X.shape[0]}]")
if not 1 <= n_neighbors < X.shape[0]:
raise ValueError(f"n_neighbors must lie in [1, {X.shape[0] - 1}]")
affinity = _knn_affinity(X, n_neighbors, sigma)
degree = affinity.sum(axis=1)
if (degree == 0.0).any():
raise ValueError(
"affinity graph has isolated vertices - sigma is too small, or "
"n_neighbors too low for this data"
)
inverse_root = 1.0 / np.sqrt(degree)
laplacian = np.eye(X.shape[0]) - inverse_root[:, None] * affinity * inverse_root[None, :]
# eigh, not eig: the Laplacian is symmetric, and the symmetric solver is
# both faster and guaranteed to return real, ordered eigenvalues.
eigenvalues, eigenvectors = np.linalg.eigh(laplacian)
embedding = eigenvectors[:, :n_clusters]
norms = np.linalg.norm(embedding, axis=1, keepdims=True)
norms[norms == 0.0] = 1.0
return SpectralEmbedding(
embedding=embedding / norms,
eigenvalues=eigenvalues[:n_clusters],
)
def fit(X: Matrix, n_clusters: int, n_neighbors: int = 10, seed: int = 0) -> Labels:
"""Embed, then cluster the rows. The k-means is stochastic; pin the seed."""
from sklearn.cluster import KMeans # the last step is ordinary k-means
result = fit_embedding(X, n_clusters, n_neighbors)
return KMeans(n_clusters=n_clusters, n_init=10, random_state=seed).fit_predict(
result.embedding
)What changed: The affinity changes from a fully-connected RBF to a symmetrized k-nearest-neighbour graph, which is not only cheaper but usually better — distant pairs contribute weight that is noise rather than structure, and a dense affinity lets that noise into the spectrum. Power iteration with deflation is replaced by a symmetric eigensolver, which returns every eigenvector at once without accumulating each one’s error into the next. The eigenvalues become part of the returned object because they are the only principled diagnostic available here, and the isolated-vertex case — which the previous stage would have divided straight through — becomes an explicit error naming its actual cause.
O(n^2 d) for the affinity and O(n^3) for the dense eigendecomposition — still quadratic in memory.
Faster, and the justification for each change.
"""Spectral clustering - sparse kNN graph, Lanczos for the bottom K only."""
import numpy as np
from numpy.typing import NDArray
from scipy import sparse
from scipy.sparse.linalg import eigsh
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def sparse_knn_affinity(
X: Matrix, n_neighbors: int, chunk: int = 2048
) -> sparse.csr_matrix:
"""Build the graph without ever materializing an n-by-n matrix.
The dense formulation needs O(n^2) memory, which is the real ceiling on
spectral clustering - not the eigendecomposition. Computing distances in
chunks and keeping only the k nearest per row bounds memory at O(n * k),
which moves the practical limit up by two orders of magnitude.
"""
design = np.ascontiguousarray(X, dtype=np.float64)
n = design.shape[0]
norms = np.einsum("ij,ij->i", design, design)
indptr = np.arange(0, (n + 1) * n_neighbors, n_neighbors) # allocated once
indices = np.empty(n * n_neighbors, dtype=np.int32)
values = np.empty(n * n_neighbors, dtype=np.float64)
collected = []
for start in range(0, n, chunk):
stop = min(start + chunk, n)
block = design[start:stop]
# ||a - b||^2 = ||a||^2 - 2 a.b + ||b||^2: one GEMM, built in place so
# the (chunk, n) block is the only intermediate that ever exists.
squared = block @ design.T
squared *= -2.0
squared += norms[None, :]
squared += np.einsum("ij,ij->i", block, block)[:, None]
np.maximum(squared, 0.0, out=squared)
squared[np.arange(stop - start), np.arange(start, stop)] = np.inf
nearest = np.argpartition(squared, n_neighbors, axis=1)[:, :n_neighbors]
block_rows = np.arange(stop - start)[:, None]
collected.append(squared[block_rows, nearest])
indices[start * n_neighbors : stop * n_neighbors] = nearest.ravel()
distances = np.concatenate(collected, axis=0)
scale = float(np.sqrt(np.median(distances)))
np.exp(-distances.ravel() / (2.0 * scale * scale), out=values)
graph = sparse.csr_matrix((values, indices, indptr), shape=(n, n))
return graph.maximum(graph.T) # symmetrize
def fit_embedding(X: Matrix, n_clusters: int, n_neighbors: int = 10) -> tuple[Matrix, Vector]:
"""Bottom K eigenvectors by shift-invert Lanczos.
A dense eigensolver computes all n eigenvectors and discards all but K.
Lanczos computes only what is asked for, and shift-invert around zero is
what makes it converge to the SMALLEST eigenvalues rather than the largest
- which is the standard trap, since the bottom of the spectrum is exactly
where Lanczos is slowest without it.
"""
affinity = sparse_knn_affinity(X, n_neighbors)
degree = np.asarray(affinity.sum(axis=1)).ravel()
if (degree == 0.0).any():
raise ValueError("affinity graph has isolated vertices; increase n_neighbors")
inverse_root = sparse.diags(1.0 / np.sqrt(degree))
# Fused: D^-1/2 W D^-1/2 is built by two sparse products, and the identity
# is subtracted by negating in place rather than forming a dense I.
normalized = inverse_root @ affinity @ inverse_root
laplacian = sparse.identity(X.shape[0], format="csr") - normalized
eigenvalues, eigenvectors = eigsh(
laplacian, k=n_clusters, sigma=0.0, which="LM"
)
norms = np.linalg.norm(eigenvectors, axis=1, keepdims=True)
norms[norms == 0.0] = 1.0
return eigenvectors / norms, eigenvaluesWhat changed: Two ceilings are removed, and they are different ceilings. The affinity is built in chunks and stored sparsely, so memory is O(n·k) rather than O(n²) — that is the constraint that actually stops spectral clustering, not the eigendecomposition. Then the dense symmetric solver is replaced by shift-invert Lanczos, which computes only the K eigenvectors wanted instead of all n; the shift-invert is the part that is easy to get wrong, because plain Lanczos converges fastest at the large end of the spectrum and this method needs the small end. The Laplacian itself is assembled through fused sparse products so no dense intermediate appears anywhere.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
Chunked distances come from one GEMM via the squared-norm expansion, and the norms are einsum reductions computed once per side rather than per pair.
Trade-off: The expansion is less numerically stable than direct subtraction, so cancellation can misrank two nearly-equidistant neighbours — which changes the graph and therefore the clustering, not just a distance value.
Fuse adjacent operations so intermediates are never materialized
The normalized Laplacian is assembled by two sparse diagonal products with the identity subtracted directly, so no dense n-by-n matrix exists at any point in the pipeline.
Trade-off: The sparse expression is far less legible than the dense formula it implements, and a sparse-times-dense slip anywhere in the chain silently reintroduces the O(n^2) allocation this exists to avoid.
Batch work to amortize interpreter overhead
Distances are computed a chunk of rows at a time, so the peak intermediate is (chunk, n) rather than (n, n) and BLAS is entered once per chunk instead of per row.
Trade-off: The chunk size is an untuned parameter, and the two-pass structure — collect distances, then compute the median bandwidth — means the data is traversed twice rather than once.
O(n^2 d) distances in BLAS with O(n*k) memory, plus Lanczos for K eigenvectors. Illustrative, not a measured benchmark.