Clustering
k-means, DBSCAN, single-linkage and spectral clustering over a data matrix, written as masks and one closure with no loop over points.
Prelude.Clustering groups the rows of a data matrix. It covers four methods: k-means (Lloyd’s
algorithm), DBSCAN, single-linkage hierarchical clustering cut at a height, and the embedding that
spectral clustering runs k-means on. None of them loops over points or clusters. “Which centroid is
nearest” is a 0/1 mask, and “which points are connected” is a matrix closed under repeated squaring,
so the only loops are over turns.
Conventions: the data x is a matrix n×p whose rows are the points. Distances are Euclidean.
Cluster labels are positions counted from zero, and a label names the lowest-numbered point that
anchors the cluster, so labels are not consecutive (a result such as [0, 0, 3, 3, 5] is three
clusters). k-means takes no hidden randomness: its starting centroids are chosen by a selection
matrix you build, so the same seed always gives the same answer. The number of Lloyd turns is fixed
by the caller; nothing checks for convergence.
Throughout this page, for a data matrix \(X \in \mathbb{R}^{n \times p}\) with rows \(x_i^{\top}\) and
centroids \(C \in \mathbb{R}^{k \times p}\) with rows \(c_l^{\top}\), \(D(X, C) \in \mathbb{R}^{n \times k}\)
is the squared-distance matrix computed by squared in Prelude.Kernels,
Limitations: complete and average linkage are not provided, because they do not reduce to a closure
the way single linkage does. A full-covariance Gaussian mixture is not provided either; the
one-dimensional mixture lives in Prelude.Fitting.
open Prelude.Clustering
k-means
Lloyd’s algorithm alternates two steps: assign every point to its nearest centroid, then move each centroid to the mean of the points assigned to it. The quantity it lowers at every turn is the inertia,
\[J = \sum_{i} \sum_{k} m_{ik} \, \lVert x_i - c_k \rVert^2\]where \(m_{ik}\) is 1 when point \(i\) is assigned to centroid \(k\). The inertia never rises from one turn to the next, which is the one guarantee k-means offers; the result is a local minimum that depends on the seed.
seeded
seeded(s, x) · s selection matrix k×n · x data n×p → matrix k×p
Let \(S \in \{0, 1\}^{k \times n}\) (the argument s) be a selection matrix whose rows are one-hot,
and let \(X \in \mathbb{R}^{n \times p}\) be the data matrix. Then
Reads starting centroids out of the data. Each row of s is one-hot: a 1 in column \(i\) picks
point \(i\) as that centroid’s start. Build s from your own random draws (or pick points by hand)
to make k-means reproducible.
Intuition Picks existing data points as starting centroids by multiplying with a one-hot selector. The randomness lives in the matrix you build, so any run can be repeated exactly.
Use it when you start kmeans. Choose k distinct points that are well spread out, whether picked by hand or from random draws, and compare several selections with inertia. Never select the same point twice: identical centroids tie for every point they are nearest to, receive the same points, and never separate. Keep the selection matrix of a good run to reproduce it.
reproducibilityclustering
def seeded(s, x) ≔ einsum("ki,ij->kj", s, x)
assigned
assigned(x, c) · x data n×p · c centroids k×p → mask n×k
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(C \in \mathbb{R}^{k \times p}\) the centroids, and let \(D = D(X, C)\). Then \(\operatorname{assigned}(X, C) = M \in \{0, 1\}^{n \times k}\) with
\[M_{il} = \big[\, D_{il} = \min_{l'} D_{il'} \,\big] .\]The assignment mask: entry \([i, k]\) is 1 when centroid \(k\) is the nearest to point \(i\) by squared Euclidean distance, 0 otherwise. Hilbert returns this mask, not a vector of indices. A point exactly equidistant from two centroids is marked for both.
Intuition Each point belongs to the centroid it is closest to, which carves space into one cell per centroid. This is the assignment half of every k-means turn, and the rule for labelling any new point against fixed centroids.
Use it when you read cluster memberships after kmeans, label new data against fitted centroids (placing new customers in existing segments), or quantise vectors by mapping each to its nearest codebook entry. The result is a mask; a tie marks several centroids, which is rare with continuous data. Scale features first, since distance weighs every column equally.
customer-segmentationclusteringdata-compression
def assigned(x, c) ≔ lowest(squared(x, c))
kmeans
kmeans(t, s, x) · t number of turns scalar · s selection matrix k×n · x data n×p → centroids k×p
Let \(t \in \mathbb{N}\), let \(S \in \{0, 1\}^{k \times n}\) be a one-hot selection matrix and \(X \in \mathbb{R}^{n \times p}\) the data matrix. Let \(C^{(0)} = S X\) and \(C^{(r+1)} = \operatorname{recentred}(X, C^{(r)})\) for \(r = 0, 1, \dots\). Then
\[\operatorname{kmeans}(t, S, X) = C^{(t)} .\]Runs t turns of Lloyd’s algorithm starting from the centroids seeded(s, x) and returns the final
centroids. Read the cluster of each point with assigned(x, kmeans(t, s, x)). A cluster that loses
all its points is not an error: its centroid moves to the origin (the numerator is zero and the
divisor is guarded to one). If that happens, reseed.
Intuition Alternate “assign each point to its nearest centre” and “move each centre to the mean of its points”. Neither half can raise the inertia, so the centres settle into a local minimum made of compact, roughly spherical groups.
Use it when you want k groups of similar points and have a guess for k: customer segmentation, colour quantisation, grouping embedded documents, or the last step of spectral clustering. Standardise features, run several seeds and keep the lowest inertia, and give enough turns that the centroids stop moving (nothing checks this). It handles elongated shapes, uneven densities and outliers badly; use dbscan or clustersAt there. Read labels with assigned.
customer-segmentationclustering
def kmeans(t, s, x) ≔ lloyd(t, x, seeded(s, x))
inertia
inertia(x, c) · x data n×p · c centroids k×p → scalar
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(C \in \mathbb{R}^{k \times p}\) the centroids, and let \(M = \operatorname{assigned}(X, C)\) and \(D = D(X, C)\). Then
\[\operatorname{inertia}(X, C) = J = \sum_{i=1}^{n} \sum_{l=1}^{k} M_{il} \, D_{il} .\]The within-cluster sum of squared distances \(J\) for the given centroids: every point’s squared distance to its nearest centroid, summed. A point tied between several nearest centroids is counted once for each of them. Use it to compare seeds or numbers of clusters.
Intuition The total squared distance from points to their nearest centroid: how tight the clusters are. Lower is tighter, but it always falls as k grows (with one centroid per point it is 0), so on its own it cannot pick k.
Use it when you compare seeds for the same k (keep the lowest), check that more turns have stopped helping, or choose k by the elbow method: plot it against k and pick where the drop flattens. Compare only runs on the same, identically scaled data.
model-evaluationhyperparameter-choiceclustering
def inertia(x, c) ≔ einsum("ik,ik->", assigned(x, c), squared(x, c))
Building blocks
clusterMass(x, c)— Let \(X \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{k \times p}\) and \(M = \operatorname{assigned}(X, C)\). Then \(\operatorname{clusterMass}(X, C)_{lj} = \sum_{i=1}^{n} M_{il}\) for every column \(j = 1, \dots, p\). Number of points assigned to each centroid, laid out k×p; the divisor of the recentring step.recentred(x, c)— Let \(X \in \mathbb{R}^{n \times p}\), \(C \in \mathbb{R}^{k \times p}\), \(M = \operatorname{assigned}(X, C)\) and \(N_l = \sum_i M_{il}\). Then \(\operatorname{recentred}(X, C)_{lj} = \dfrac{\sum_{i=1}^{n} M_{il} X_{ij}}{N_l + [N_l = 0]}\). One recentring step: each centroid becomes the mean of its assigned points (the origin if it has none).lloyd(t, x, c)— Let \(t \in \mathbb{N}\), \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{k \times p}\). Then \(\operatorname{lloyd}(0, X, C) = C\) and \(\operatorname{lloyd}(t, X, C) = \operatorname{lloyd}(t - 1, X, \operatorname{recentred}(X, C))\). Repeatsrecentredttimes from centroidsc; used bykmeans.
The closure
DBSCAN and single-linkage clustering share one piece of machinery: the transitive closure of a 0/1 relation on the points. Two points belong together when a chain of links joins them, and the update
\[A \leftarrow \big[\, A A + A > 0 \,\big]\]doubles the length of the chains it has found at each turn, so \(\lceil \log_2 n \rceil\) turns close any relation on \(n\) points. A component is then named by the lowest position it contains.
joined
joined(x, a) · x data n×p (only its row count is read) · a relation mask n×n → mask n×n
Let \(X \in \mathbb{R}^{n \times p}\) (only \(n\) is read) and let \(A \in \{0, 1\}^{n \times n}\) be a relation mask. Let \(A_0 = A\) and \(A_{r+1} = [\, A_r A_r + A_r > 0 \,]\) entrywise, for \(r = 0, 1, \dots\), so that
\[(A_r)_{ij} = \Big[\, \textstyle\sum_{\ell=1}^{2^r} (A^{\ell})_{ij} > 0 \,\Big] .\]Then, with \(R = \lceil \log_2 \max(n, 2) \rceil\),
\[\operatorname{joined}(X, A) = A_R .\]The transitive closure of the relation a: entry \([i, j]\) is 1 when a chain of links in a
leads from \(i\) to \(j\). For the closure to mark each point as reaching itself, a should have
ones on its diagonal.
Intuition If \(i\) links to \(j\) and \(j\) links to \(l\), then \(i\) reaches \(l\). Each squaring turn doubles the length of the chains found, so after \(\lceil \log_2 n \rceil\) turns every chain is known and the rows mark whole connected groups.
Use it when you have any pairwise “belongs with” relation and need the groups it implies transitively: records sharing an email or phone number (entity resolution), connectivity in a network, pixels joined to similar neighbours. Put ones on the diagonal, and make the relation symmetric if links have no direction. Name each group by the lowest position in its row, as dbscan and clustersAt do. It costs about \(\log_2 n\) products of n×n matrices.
entity-resolutionconnectivity-analysisclustering
def joined(x, a) ≔ joinLoop(rounds(x), a)
Building blocks
rounds(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rounds}(X) = \left\lceil \ln \max(n, 2) / \ln 2 \right\rceil = \lceil \log_2 \max(n, 2) \rceil\). The number of squaring turns needed, the ceiling of the base-2 logarithm of the row count (at least 1).joinLoop(t, a)— Let \(t \in \mathbb{N}\) and \(A \in \{0, 1\}^{n \times n}\). Then \(\operatorname{joinLoop}(0, A) = A\) and \(\operatorname{joinLoop}(t, A) = \operatorname{joinLoop}(t - 1, [\, A A + A > 0 \,])\). Applies the squaring updatettimes.indexOf(x)— Let \(X \in \mathbb{R}^{n \times p}\). Then \(\operatorname{indexOf}(X) = (0, 1, \dots, n - 1) \in \mathbb{R}^n\). The positions 0, 1, …, n−1 of the rows ofx, as a vector.spanIndex(m, x)— Let \(M \in \mathbb{R}^{n \times n}\) (the argumentm) and \(X \in \mathbb{R}^{n \times p}\), with positions \(0, 1, \dots, n-1\). Then \(\operatorname{spanIndex}(M, X)_{ij} = j\). The position vector copied into every row of an n×n matrix.nameOf(m, x)— Let \(M \in \{0, 1\}^{n \times n}\) (the argumentm) and \(X \in \mathbb{R}^{n \times p}\), with positions \(0, 1, \dots, n-1\). Then \(\operatorname{nameOf}(M, X)_i = \min_{j} \big( M_{ij}\, j + (1 - M_{ij}) \cdot 10^{300} \big)\). For each row of the maskm, the lowest column position holding a 1; a row with no 1 gets 1e300.
DBSCAN
DBSCAN finds clusters as dense regions. A point is a core point when at least m points,
itself included (the scikit-learn convention), lie within distance eps of it. Core points within
eps of each other are in the same cluster, taken transitively; a non-core border point joins the
cluster of any core point within eps; everything else is noise. Noise is a result, not an error.
The radius parameter is called eps rather than e, because e is Euler’s number in this
language.
In this section let \(X \in \mathbb{R}^{n \times p}\) be the data matrix, \(\varepsilon > 0\) (the
argument eps) the radius and \(m > 0\) the minimum neighbourhood size; write
\(N = \operatorname{near}(X, \varepsilon)\), \(\kappa = \operatorname{coreOf}(X, \varepsilon, m)\) and
\(B = \operatorname{attached}(X, \varepsilon, m)\), and count positions \(0, 1, \dots, n-1\).
coreOf
coreOf(x, eps, m) · x data n×p · eps radius scalar · m minimum neighbourhood size scalar → mask n
Let \(X \in \mathbb{R}^{n \times p}\), \(\varepsilon > 0\) and \(m > 0\), and let \(N = \operatorname{near}(X, \varepsilon)\). Then
\[\operatorname{coreOf}(X, \varepsilon, m)_i = \kappa_i = \Big[\, \sum_{j=1}^{n} N_{ij} \ge m \,\Big] .\]1 for each core point: a point with at least m points (counting itself) within distance eps.
Intuition A core point stands in a crowd: at least m points, itself included, within radius eps. Core points form the interior of dense regions, and clusters grow outward from them.
Use it when you inspect DBSCAN’s decisions or tune its parameters. If almost nothing is core, eps is too small or m too large; if nearly everything is core, clusters will merge. A common start is m about twice the number of features (at least 3), with eps read from the elbow of sorted kDistance(x, m − 1) values in Prelude.Neighbours, square-rooted because those are squared distances.
hyperparameter-choiceanomaly-detectionclustering
def coreOf(x, eps, m) ≔ einsum("ij->i", near(x, eps)) ≥ m
dbscan
dbscan(x, eps, m) · x data n×p · eps radius scalar · m minimum neighbourhood size scalar → labels n
Let \(X \in \mathbb{R}^{n \times p}\), \(\varepsilon > 0\) and \(m > 0\), and let \(B = \operatorname{attached}(X, \varepsilon, m)\). Then, with positions \(j = 0, 1, \dots, n-1\),
\[\operatorname{dbscan}(X, \varepsilon, m)_i = \min_{j} \big( B_{ij}\, j + (1 - B_{ij}) \cdot 10^{300} \big) .\]The DBSCAN cluster label of every point: the position of the lowest-numbered core point in its
cluster. A border point within reach of two clusters takes the lower label. Noise points get the
label 1e300; test for them with dbNoise rather than by comparing labels.
Intuition Clusters are dense regions separated by sparse gaps. Core points within eps of one another chain into a cluster of whatever shape the chain traces, and points in no dense region are called noise instead of being forced into a cluster.
Use it when the number of clusters is unknown, clusters have irregular shapes (geographic hot spots, stops in GPS traces), or outliers should stay out of every group. One eps must suit all clusters, so it struggles when their densities differ a lot. Scale features first. Labels are not consecutive and noise is labelled 1e300; detect it with dbNoise. Memory is n×n and the closure takes about \(\log_2 n\) matrix products.
clusteringgeospatial-analysisanomaly-detection
def dbscan(x, eps, m) ≔ nameOf(attached(x, eps, m), x)
dbNoise
dbNoise(x, eps, m) · x data n×p · eps radius scalar · m minimum neighbourhood size scalar → mask n
Let \(X \in \mathbb{R}^{n \times p}\), \(\varepsilon > 0\) and \(m > 0\), and let \(B = \operatorname{attached}(X, \varepsilon, m)\). Then
\[\operatorname{dbNoise}(X, \varepsilon, m)_i = \Big[\, \sum_{j=1}^{n} B_{ij} = 0 \,\Big] .\]1 for every point that is neither a core point nor within eps of one.
Intuition The points DBSCAN declines to cluster: not core, and not within eps of any core point. They are isolated relative to the density you asked for.
Use it when you remove outliers before another analysis, flag anomalies, or tune: a large noise fraction means eps is too small or m too large. It is safer than comparing labels with the 1e300 sentinel.
anomaly-detectiondata-cleaning
def dbNoise(x, eps, m) ≔ einsum("ij->i", attached(x, eps, m)) = 0
Building blocks
near(x, eps)— Let \(X \in \mathbb{R}^{n \times p}\), \(\varepsilon > 0\) and \(D = D(X, X)\). Then \(\operatorname{near}(X, \varepsilon)_{ij} = [\, D_{ij} \le \varepsilon^2 \,]\). n×n mask of pairs at distance at mosteps(the diagonal included).coreEdges(x, eps, m)— Let \(N\) and \(\kappa\) be as above. Then \(\operatorname{coreEdges}(X, \varepsilon, m)_{ij} = N_{ij}\, \kappa_i\, \kappa_j\). Thenearmask restricted to pairs of core points.coreGroups(x, eps, m)— Let \(E = \operatorname{coreEdges}(X, \varepsilon, m)\). Then \(\operatorname{coreGroups}(X, \varepsilon, m) = \operatorname{joined}(X, E)\). The closure ofcoreEdges: which core points share a cluster.toCore(x, eps, m)— Let \(N\) and \(\kappa\) be as above. Then \(\operatorname{toCore}(X, \varepsilon, m)_{ij} = N_{ij}\, \kappa_j\). n×n mask from every point to each core point withineps.attached(x, eps, m)— Let \(T = \operatorname{toCore}(X, \varepsilon, m)\) and \(G = \operatorname{coreGroups}(X, \varepsilon, m)\). Then \(\operatorname{attached}(X, \varepsilon, m)_{ik} = \big[\, (T G)_{ik} > 0 \,\big] = \big[\, \sum_j T_{ij} G_{jk} > 0 \,\big]\). n×n mask from every point to every core point of the cluster it belongs to.
Single linkage
Cutting a single-linkage dendrogram at height \(h\) gives exactly the connected components of the graph that joins every pair of points at distance at most \(h\). So these definitions need no merge loop and never build a dendrogram: they close that graph and read the components off. To see the hierarchy, call them at several heights.
clustersAt
clustersAt(x, h) · x data n×p · h cut height scalar → labels n
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(h \ge 0\) the cut height, and let \(G = \operatorname{groupsAt}(X, h)\). Then, with positions \(j = 0, 1, \dots, n-1\),
\[\operatorname{clustersAt}(X, h)_i = \min_{j} \big( G_{ij}\, j + (1 - G_{ij}) \cdot 10^{300} \big) .\]The single-linkage cluster of every point when the dendrogram is cut at height h. A label is the
position of the lowest-numbered point in the cluster.
Intuition Two points share a group whenever a chain of hops, none longer than h, connects them. Cutting at h therefore returns the groups separated by gaps wider than h, whatever their shape.
Use it when groups are separated by clear gaps and may be elongated or curved: splitting a track into segments, grouping locations within walking distance, or merging near-identical records with a small h. It suffers from chaining: a thin trail of points can bridge two real clusters into one, which dbscan resists by requiring density. Call clusterCount over several heights to see the hierarchy.
clusteringduplicate-detection
def clustersAt(x, h) ≔ nameOf(groupsAt(x, h), x)
clusterCount
clusterCount(x, h) · x data n×p · h cut height scalar → scalar
Let \(X \in \mathbb{R}^{n \times p}\) and \(h \ge 0\), and let \(\ell = \operatorname{clustersAt}(X, h)\), with points at positions \(i = 0, 1, \dots, n-1\). Then
\[\operatorname{clusterCount}(X, h) = \sum_{i=0}^{n-1} [\, \ell_i = i \,] .\]How many clusters the cut at height h gives: each cluster is counted once, at the point that names
it.
Intuition The number of groups at gap threshold h. As h grows, groups merge and the count falls to 1, so the count against h is a compact summary of the dendrogram.
Use it when you choose a cut height: a long flat stretch of the count over a range of h marks a grouping that does not depend on the exact threshold. It also finds the smallest h that yields a target number of groups. Every call recomputes the whole closure, so sweep a handful of heights rather than thousands.
hyperparameter-choiceclustering
def clusterCount(x, h) ≔ einsum("i->", nameOf(groupsAt(x, h), x) = indexOf(x))
Building blocks
closerThan(x, h)— Let \(X \in \mathbb{R}^{n \times p}\), \(h \ge 0\) and \(D = D(X, X)\). Then \(\operatorname{closerThan}(X, h)_{ij} = [\, D_{ij} \le h^2 \,]\). n×n mask of pairs at distance at mosth.groupsAt(x, h)— Let \(X \in \mathbb{R}^{n \times p}\) and \(h \ge 0\). Then \(\operatorname{groupsAt}(X, h) = \operatorname{joined}(X, \operatorname{closerThan}(X, h))\). The closure ofcloserThan: which points share a cluster at heighth.
Spectral
Spectral clustering embeds the points using eigenvectors of a normalised graph Laplacian built from a Gaussian affinity, then clusters the embedded rows with k-means. With affinity \(A\) and degrees \(d_i = \sum_j A_{ij}\),
\[L = I - D^{-1/2} A \, D^{-1/2}\]This section gives the affinity, the Laplacian and its eigenvectors. The eigenvectors come from
Prelude.LinearAlgebra, which orders them by descending eigenvalue magnitude, so the eigenvectors
of the smallest eigenvalues (the ones spectral clustering uses) are the last columns. Pick those
columns and pass the resulting matrix to kmeans. The eigensolver runs a fixed 40 turns, so nearly
equal eigenvalues are resolved less accurately than a shifted solver would.
affinity
affinity(x, g) · x data n×p · g bandwidth scalar → matrix n×n
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(g > 0\) the bandwidth, and let \(D = D(X, X)\). Then
\[\operatorname{affinity}(X, g)_{ij} = A_{ij} = \exp(-g\, D_{ij}) = \exp\!\big(-g \lVert x_i - x_j \rVert^2\big) .\]A larger g makes the affinity fall off faster with distance. The diagonal is 1.
Intuition Turns distances into similarities between 0 and 1: 1 for identical points, falling off like a bell curve. \(1/\sqrt{g}\) is the length scale: at that distance the affinity is \(e^{-1} \approx 0.37\), and pairs several times further apart are effectively unconnected.
Use it when you build the similarity graph for graphLaplacian and spectral clustering, or a similarity heat map. A common start sets \(1/g\) near the median squared distance between points; scale features first. Too large a g isolates every point, too small a g connects everything almost equally, and both hide the clusters.
similarity-searchclustering
def affinity(x, g) ≔ rbf(x, x, g)
graphLaplacian
graphLaplacian(a) · a symmetric affinity n×n → matrix n×n
Let \(A \in \mathbb{R}^{n \times n}\) be a symmetric non-negative affinity matrix and \(d = \operatorname{nodeDegrees}(A)\). Then
\[\operatorname{graphLaplacian}(A)_{ij} = L_{ij} = [\, i = j \,] - \frac{A_{ij}}{\sqrt{\max(d_i d_j,\, 10^{-300})}} .\]The normalised Laplacian \(L\) of the weighted graph a. The product of degrees \(d_i d_j\) is floored
at \(10^{-300}\), so an isolated node does not produce a division by zero; wherever that floor is not
reached this is exactly \(I - D^{-1/2} A D^{-1/2}\).
Intuition Measures how much a signal on the graph varies across strong links, normalised by each node’s degree. A group with strong links inside and weak links outside produces an eigenvalue near 0 whose eigenvector is nearly constant across the group (up to degree scaling), so the small-eigenvalue eigenvectors reveal the groups.
Use it when you do spectral clustering or community detection in a weighted network (social graphs, co-purchase graphs). Its eigenvalues lie in \([0, 2]\), and the number close to 0 suggests how many well-separated groups exist. a must be symmetric and non-negative.
community-detectionclustering
def graphLaplacian(a) ≔ identity(a) − a / sqrt(max(einsum("i,j->ij", nodeDegrees(a), nodeDegrees(a)), 1e-300))
spectralEmbedding
spectralEmbedding(x, g) · x data n×p · g bandwidth scalar → matrix n×n
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(g > 0\) the bandwidth, and let \(L = \operatorname{graphLaplacian}(\operatorname{affinity}(X, g))\). Then
\[\operatorname{spectralEmbedding}(X, g) = \operatorname{eigenvectors}(L) = V_{40},\]where \(L_0 = L\), \(V_0 = I_n\), and \(L_{r+1} = Q_r^{\top} L_r Q_r\), \(V_{r+1} = V_r Q_r\) with
\(Q_r\) the orthogonal factor of the QR decomposition of \(L_r\) (the unshifted QR algorithm of
Prelude.LinearAlgebra).
The eigenvectors of the normalised Laplacian of affinity(x, g), one per column. Row \(i\) is point
\(i\)’s embedding; keep the last few columns (the smallest eigenvalues) before clustering.
Intuition Replaces each point by its entries in the Laplacian’s eigenvectors. Points joined by chains of strong similarity land close together even when their clusters are curved or intertwined in the original space, so kmeans can then separate them with straight boundaries.
Use it when clusters are non-convex (rings, half-moons, spirals) and n is moderate, since it solves an n×n eigenproblem. For k clusters, keep the last k columns, optionally scale each row to unit length, then run kmeans. Results depend strongly on g. Clean clusters give several eigenvalues near 0, exactly the nearly equal case the fixed 40-turn solver resolves least accurately, so check the result on data whose grouping you know.
clusteringimage-segmentation
def spectralEmbedding(x, g) ≔ eigenvectors(graphLaplacian(affinity(x, g)))
Building blocks
nodeDegrees(a)— Let \(A \in \mathbb{R}^{n \times n}\). Then \(\operatorname{nodeDegrees}(A)_i = d_i = \sum_{j=1}^{n} A_{ij}\). The row sums of an affinity matrix: each node’s weighted degree.
Examples
Split six points into two groups with k-means
Two tight triples of points, seeded with the first two points as starting centroids, run for ten turns.
open Prelude.Clustering
fn centroidsAfter(t, s, x) ↦ kmeans(t, s, x)
fn spreadAround(x, c) ↦ inertia(x, c)
fn membership(x, c) ↦ assigned(x, c)
using Examples;
var x = new double[,] { { 0, 0 }, { 0, 1 }, { 1, 0 }, { 10, 10 }, { 10, 11 }, { 11, 10 } };
var s = new double[,] { { 1, 0, 0, 0, 0, 0 }, { 0, 1, 0, 0, 0, 0 } }; // start from points 0 and 1
var c = Groups.CentroidsAfter(10.0, s, x); // class = the .hb file name (Groups.hb)
var j = Groups.SpreadAround(x, c);
var m = Groups.Membership(x, c);
Output (verified, rounded to 6 significant digits):
centroids = [[0.333333, 0.333333], [10.3333, 10.3333]]
inertia = 2.66667
assigned = [[1, 0], [1, 0], [1, 0], [0, 1], [0, 1], [0, 1]]
Both starting centroids were in the first triple, yet the second one moves to the far group within a few turns. Each triple contributes \(2/9 + 5/9 + 5/9 = 4/3\) to the inertia.
Find dense groups and an outlier with DBSCAN
Two unit squares of points and one far-away point, with radius 1.5 and at least three points per neighbourhood.
open Prelude.Clustering
fn densityLabels(x, eps, m) ↦ dbscan(x, eps, m)
fn densityNoise(x, eps, m) ↦ dbNoise(x, eps, m)
using Examples;
var pts = new double[,] {
{ 0, 0 }, { 0, 1 }, { 1, 0 }, { 1, 1 },
{ 5, 5 }, { 5, 6 }, { 6, 5 }, { 6, 6 },
{ 20, 20 } };
var labels = Groups.DensityLabels(pts, 1.5, 3.0);
var noise = Groups.DensityNoise(pts, 1.5, 3.0);
Output (verified; the last label is written here in exponent form):
labels = [0, 0, 0, 0, 4, 4, 4, 4, 1e300]
noise = [0, 0, 0, 0, 0, 0, 0, 0, 1]
The first square is named after point 0 and the second after point 4. The lone point is noise, and its label is the 1e300 sentinel.
Cut a single-linkage dendrogram at several heights
Six points on a line at 0, 1, 2, 6, 7 and 15.
open Prelude.Clustering
fn linkageLabels(x, h) ↦ clustersAt(x, h)
fn linkageCount(x, h) ↦ clusterCount(x, h)
using Examples;
var line = new double[,] { { 0 }, { 1 }, { 2 }, { 6 }, { 7 }, { 15 } };
var labels = Groups.LinkageLabels(line, 1.5);
var at1_5 = Groups.LinkageCount(line, 1.5);
var at5 = Groups.LinkageCount(line, 5.0);
var at10 = Groups.LinkageCount(line, 10.0);
Output (verified):
clustersAt(1.5) = [0, 0, 0, 3, 3, 5]
clusterCount(1.5) = 3
clusterCount(5) = 2
clusterCount(10) = 1
At height 1.5 the groups are {0, 1, 2}, {6, 7} and {15}. The gap of 4 closes by height 5 and the gap of 8 by height 10.