Manifold
Dimension reduction and embedding: classical scaling, kernel PCA, Isomap, truncated SVD, non-negative factorisation, random projection and t-SNE.
Prelude.Manifold describes a cloud of points from fewer dimensions. It covers classical
multidimensional scaling (from a distance matrix), kernel principal components (from a kernel
matrix), Isomap (scaling on distances measured along a neighbourhood graph), truncated SVD for
latent semantic analysis, non-negative matrix factorisation, Johnson–Lindenstrauss random projection
and the gradient of t-distributed stochastic neighbour embedding.
Conventions: a data matrix x is n×p with one point per row; a distance matrix d is n×n and
holds plain (not squared) Euclidean distances; a kernel matrix k is n×n. The eigen-based methods
return every component as a column, ordered by descending eigenvalue magnitude, so “reduce to two
dimensions” means reading the first two columns. Their eigensolver runs a fixed 40 turns, which is
exact to machine precision for a well-separated spectrum but less accurate when two eigenvalues are
nearly equal. Nothing here draws random numbers: the projection matrix of a random projection, the
starting factors of NMF and the starting layout of t-SNE are arguments you supply.
Notation used below: for a symmetric \(B \in \mathbb{R}^{n \times n}\), \(\operatorname{eigenvectors}(B) = V\)
is the matrix whose columns are the eigenvectors found by the 40-turn QR iteration of
Prelude.LinearAlgebra, and \(\operatorname{eigenvalues}(B) = \lambda\) with
\(\lambda_j = (V^{\top} B V)_{jj}\). Squared distances \(\lVert x_i - x_j \rVert^2\) are computed as
\(\lVert x_i \rVert^2 + \lVert x_j \rVert^2 - 2 x_i^{\top} x_j\) (squared in Prelude.Kernels).
Not provided: locally linear embedding, factor analysis and FastICA (each needs a batch of small solves or a growing basis), and UMAP. The t-SNE here uses one bandwidth for every point rather than a per-point perplexity search.
open Prelude.Manifold
Classical multidimensional scaling
Given distances \(D\), double-centre the squared distances and halve them to get a Gram matrix,
\[B = -\tfrac{1}{2} \, J \, D^{(2)} J, \qquad J = I - \tfrac{1}{n} \mathbf{1}\mathbf{1}^\top\]where \(D^{(2)}\) squares each entry. The eigenvectors of \(B\) scaled by the square roots of its eigenvalues are coordinates that reproduce the distances exactly whenever they were Euclidean to begin with. The double-centring is the same operation kernel PCA applies to a kernel.
mdsCoords
mdsCoords(d) · d distance matrix n×n → coordinates n×n
Let \(D \in \mathbb{R}^{n \times n}\) (the argument d) be a symmetric matrix of distances, let
\(B = \operatorname{mdsGram}(D)\), \(V = \operatorname{eigenvectors}(B)\) and
\(\lambda = \operatorname{eigenvalues}(B)\). Then
The classical scaling coordinates: row \(i\) is point \(i\), column \(j\) is the \(j\)-th axis, axes ordered by eigenvalue magnitude. Keep the first few columns. A negative eigenvalue (distances that are not Euclidean) gives an axis of zeros rather than an imaginary coordinate. The coordinates are determined up to the sign of each axis.
Intuition Draws a map from a table of distances: it finds positions whose straight-line distances reproduce the table. The first axis follows the direction of greatest spread and the next the greatest remaining one, as in PCA. Distances measured between points in p dimensions give at most p non-zero axes.
Use it when all you have are pairwise dissimilarities (road distances between cities, similarity ratings, distances between sequences) and you want a 2-D or 3-D picture, or coordinates for a method that needs them. Keep the first two or three columns and check with mdsStrength how much they capture. It solves an n×n eigenproblem. If the points follow a curved surface, use isomap.
visualisationdimensionality-reduction
def mdsCoords(d) ≔ eigenvectors(mdsGram(d)) · sqrt(max(eigenvalues(mdsGram(d)), 0))
mdsStrength
mdsStrength(d) · d distance matrix n×n → vector n
Let \(D \in \mathbb{R}^{n \times n}\) be a symmetric matrix of distances and \(B = \operatorname{mdsGram}(D)\). Then
\[\operatorname{mdsStrength}(D) = \operatorname{eigenvalues}(B) = \lambda \in \mathbb{R}^n .\]The eigenvalues of the Gram matrix \(B\), in the same order as the columns of mdsCoords. Each is
the variance an axis carries; significantly negative values show the distances are not Euclidean.
Intuition How much spread each axis of mdsCoords carries. A sharp drop after the first few values says the points really live in that many dimensions.
Use it when you decide how many scaling axes to keep: the share of the positive values covered by the first few plays the role of “variance explained” in PCA. Large negative values mean the dissimilarities cannot be drawn exactly in any flat space, so every map will distort some of them.
dimensionality-reductionhyperparameter-choicemodel-evaluation
def mdsStrength(d) ≔ eigenvalues(mdsGram(d))
Building blocks
mdsGram(d)— Let \(D \in \mathbb{R}^{n \times n}\) and \(S = D \odot D\), with row means \(\bar{s}_{i \cdot} = \frac{1}{n}\sum_j S_{ij}\), column means \(\bar{s}_{\cdot j} = \frac{1}{n}\sum_i S_{ij}\) and grand mean \(\bar{s}\). Then \(\operatorname{mdsGram}(D)_{ij} = -\tfrac{1}{2}\big(S_{ij} - \bar{s}_{i \cdot} - \bar{s}_{\cdot j} + \bar{s}\big)\), which is \(-\tfrac{1}{2} J S J\). The double-centred Gram matrix of the squared distances, negated and halved.
Kernel principal components
Principal component analysis in the feature space of a kernel, without ever building the features:
centring the kernel matrix centres the feature map, so the method is the eigenproblem of the centred
kernel. With a linear kernel it gives the same answer as ordinary PCA in Prelude.Multivariate.
Build the kernel with Prelude.Kernels (for example rbf(x, x, g)).
In this section let \(K \in \mathbb{R}^{n \times n}\) (the argument k) be a kernel matrix and
\(\tilde{K} = J K J\) its centred form, \(\tilde{K}_{ij} = K_{ij} - \bar{k}_{i \cdot} - \bar{k}_{\cdot j} + \bar{k}\)
(row, column and grand means, as computed by centreKernel in Prelude.Kernels).
kpcaAxes
kpcaAxes(k) · k kernel matrix n×n → matrix n×n
Let \(K \in \mathbb{R}^{n \times n}\) be a symmetric kernel matrix and \(\tilde{K}\) its centred form. Then
\[\operatorname{kpcaAxes}(K) = V = \operatorname{eigenvectors}(\tilde{K}) .\]The eigenvectors of the centred kernel, one per column, ordered by eigenvalue magnitude.
Intuition PCA carried out after a non-linear map into the feature space the kernel defines, without ever computing that map. Each column gives every training point’s loading on one non-linear direction of variation.
Use it when you need the components themselves: to inspect them, or to project new points, which also needs the eigenvalues and the kernel between new and training points (this module has no function for that projection). For the coordinates of the training points, kpcaScores is what you usually want.
dimensionality-reductionfeature-extraction
def kpcaAxes(k) ≔ eigenvectors(centreKernel(k))
kpcaScores
kpcaScores(k) · k kernel matrix n×n → scores n×n
Let \(K \in \mathbb{R}^{n \times n}\) be a symmetric kernel matrix, \(V = \operatorname{kpcaAxes}(K)\) and \(\lambda = \operatorname{eigenvalues}(\tilde{K})\). Then
\[\operatorname{kpcaScores}(K)_{ij} = V_{ij} \sqrt{\max(\lambda_j, 0)} .\]The kernel principal component scores of the training points: the axes scaled by the square roots of the (non-negative part of the) eigenvalues. Row \(i\) is point \(i\); keep the first few columns.
Intuition The training points’ coordinates along the kernel principal components: each axis scaled by the square root of its eigenvalue, so axes with more variance spread the points further. With a linear kernel these are ordinary PCA scores, up to sign.
Use it when you want to unfold non-linear structure for a plot or as features: concentric rings, curved clusters. An RBF kernel rbf(x, x, g) is the usual choice, and its width changes the picture a great deal, so try several. The eigenproblem is n×n, so cost grows as \(n^3\). For straight-line structure use ordinary PCA in Prelude.Multivariate; for points on one curved sheet, compare with isomap.
dimensionality-reductionfeature-extractionvisualisation
def kpcaScores(k) ≔ kpcaAxes(k) · sqrt(max(eigenvalues(centreKernel(k)), 0))
Isomap
Straight-line distance says that two points on opposite ends of a curved sheet are close; Isomap
measures distance along the sheet instead. It links each point to its k nearest neighbours, finds
shortest paths over that graph, and runs classical scaling on the path lengths.
The shortest paths use the min-plus form of matrix squaring,
\[d_{ij} \leftarrow \min_k \,( d_{ik} + d_{kj} )\]Each turn doubles the number of hops a path may take, so \(\lceil \log_2 n \rceil\) turns give every geodesic, at a cost of \(n^3\) per turn. A missing edge is represented as \(10^{12}\) (not \(10^{300}\), which would overflow when two are added), so this is meant for data whose real distances are far below that; a pair of points in disconnected parts of the graph keeps a distance near \(10^{12}\).
In this section let \(X \in \mathbb{R}^{n \times p}\) be the data matrix with rows \(x_i^{\top}\) and \(k \in \mathbb{N}\) the number of neighbours per point.
geodesics
geodesics(x, k) · x data n×p · k neighbours per point scalar → distance matrix n×n
Let \(X \in \mathbb{R}^{n \times p}\) and \(k \in \mathbb{N}\), and let \(G = \operatorname{isoGraph}(X, k)\). Let \(G_0 = G\) and \((G_{r+1})_{ij} = \min_{l} \big( (G_r)_{il} + (G_r)_{lj} \big)\) for \(r = 0, 1, \dots\). Then, with \(R = \lceil \log_2 \max(n, 2) \rceil\),
\[\operatorname{geodesics}(X, k) = G_R .\]Shortest-path distances over the neighbourhood graph. Point \(i\) is linked to its k nearest other
points (ties at the k-th distance are all kept) with an edge as long as their Euclidean distance.
Links are taken from each point’s own neighbour list and are not symmetrised.
Intuition Distance measured by walking only between near neighbours, like driving along roads instead of flying. On a rolled-up sheet, two points close through the air but far apart along the sheet get a large geodesic distance.
Use it when the data lies on a curved, low-dimensional surface (a swiss roll, a sequence of poses, a sensor trajectory) and you need distances along it, usually for mdsCoords. Choose k so the graph is connected: disconnected pieces show up as distances near \(10^{12}\). Too large a k adds shortcuts across folds. Each turn builds an n×n×n tensor.
dimensionality-reductionnetwork-distances
def geodesics(x, k) ≔ pathLoop(hops(distance(x, x)), isoGraph(x, k))
isomap
isomap(x, k) · x data n×p · k neighbours per point scalar → coordinates n×n
Let \(X \in \mathbb{R}^{n \times p}\) and \(k \in \mathbb{N}\). Then
\[\operatorname{isomap}(X, k) = \operatorname{mdsCoords}\big(\operatorname{geodesics}(X, k)\big) .\]The Isomap embedding: mdsCoords applied to geodesics(x, k). Keep the first few columns. In
practice, call geodesics and mdsCoords from two separate doorways, as in the example below: a
single doorway that names isomap currently expands into very large generated code.
Intuition Classical scaling applied to geodesic distances: it flattens a curved sheet onto a plane while keeping distances measured along the sheet, so a swiss roll unrolls into a rectangle and the U in the example straightens into a line.
Use it when you want to unfold one connected, curved manifold into a few coordinates: recovering an angle or pose from images, or parametrising a trajectory or curve. It needs a connected neighbour graph, and noise that creates shortcuts or holes in the sheet distort the result. For pictures of clusters rather than of one surface, embedded (t-SNE) is the better tool. Call geodesics and mdsCoords separately, as noted above.
dimensionality-reductionvisualisation
def isomap(x, k) ≔ mdsCoords(geodesics(x, k))
Building blocks
isoLinks(x, k)— Let \(X \in \mathbb{R}^{n \times p}\), \(k \in \mathbb{N}\), \(\tilde{D}_{ij} = \lVert x_i - x_j \rVert^2 + 10^{300} [\, i = j \,]\) and \(N_{ij} = \big[\, \#\{ l : \tilde{D}_{il} < \tilde{D}_{ij} \} < k \,\big]\) (the maskneighboursIn(x, k)ofPrelude.Neighbours). Then \(\operatorname{isoLinks}(X, k)_{ij} = [\, N_{ij} + [i = j] > 0 \,]\). n×n mask of each point’sknearest neighbours, plus the point itself.isoGraph(x, k)— Let \(L = \operatorname{isoLinks}(X, k)\). Then \(\operatorname{isoGraph}(X, k)_{ij} = \lVert x_i - x_j \rVert\, L_{ij} + (1 - L_{ij}) \cdot 10^{12}\), where the distance is \(\sqrt{\max(\lVert x_i - x_j \rVert^2, 0)}\). Edge lengths: the Euclidean distance on a link, 1e12 where there is none.viaK(d)— Let \(D \in \mathbb{R}^{n \times n}\). Then \(\operatorname{viaK}(D) \in \mathbb{R}^{n \times n \times n}\) with \(\operatorname{viaK}(D)_{ikj} = D_{ik} + D_{kj}\). The n×n×n tensor of path lengths through every intermediate point.relaxed(d)— Let \(D \in \mathbb{R}^{n \times n}\). Then \(\operatorname{relaxed}(D)_{ij} = \min_{k} \big( D_{ik} + D_{kj} \big)\). One min-plus squaring turn: the shortest route through any one intermediate point.hops(d)— Let \(D \in \mathbb{R}^{n \times n}\). Then \(\operatorname{hops}(D) = \left\lceil \ln \max(n, 2) / \ln 2 \right\rceil = \lceil \log_2 \max(n, 2) \rceil\). The number of turns needed, the ceiling of the base-2 logarithm of the row count (at least 1).pathLoop(t, d)— Let \(t \in \mathbb{N}\) and \(D \in \mathbb{R}^{n \times n}\). Then \(\operatorname{pathLoop}(0, D) = D\) and \(\operatorname{pathLoop}(t, D) = \operatorname{pathLoop}(t - 1, \operatorname{relaxed}(D))\). Appliesrelaxedttimes.
Truncated SVD
Latent semantic analysis is the singular value decomposition of a document-term matrix \(A = U \Sigma V^\top\) (documents in rows, terms in columns). The right singular vectors are the topics, the singular values their weights, and \(AV = U\Sigma\) places each document among the topics. “Truncated” is only a choice of how many leading columns to read. The singular vectors come from the eigenproblem of \(A^\top A\), so the smallest singular values of a badly conditioned matrix are less accurate than a dedicated SVD routine would give.
topics
topics(a) · a document-term matrix n×p → matrix p×p
Let \(A \in \mathbb{R}^{n \times p}\) (the argument a) be a document-term matrix. Then
The right singular vectors \(V\), one topic per column, ordered by singular value.
Intuition Each column is a weighted mix of terms that tend to occur in the same documents. The first captures the strongest co-occurrence pattern and each later one the strongest remaining pattern at right angles to the earlier ones, so words used in similar documents load on the same topics.
Use it when you do latent semantic analysis: finding themes in a document-term matrix, building compact document vectors for search, or matching a query to documents that use different words for the same idea. Weighting terms (for example by tf-idf) before the decomposition usually helps. Topics may mix positive and negative weights and so can be hard to read; for additive, readable parts use nmfMix. It works on a p×p matrix, so vocabulary size drives the cost.
topic-modellingsemantic-searchdimensionality-reduction
def topics(a) ≔ rightVectors(a)
topicWeights
topicWeights(a) · a document-term matrix n×p → vector p
Let \(A \in \mathbb{R}^{n \times p}\) be a document-term matrix and \(\mu = \operatorname{eigenvalues}(A^{\top} A)\). Then
\[\operatorname{topicWeights}(A)_j = \sigma_j = \sqrt{\max(\mu_j, 0)}, \qquad j = 1, \dots, p .\]The singular values \(\Sigma\), in descending order.
Intuition How much of the matrix each topic accounts for: the squares of these values add up to the sum of squares of every entry, so their fall-off shows how many topics carry real signal.
Use it when you decide how many columns of topics and documentsIn to keep: look for an elbow, or keep enough to cover a chosen share of the sum of squares. The number of values clearly above zero is the matrix’s effective rank.
topic-modellinghyperparameter-choicedimensionality-reduction
def topicWeights(a) ≔ singular(a)
documentsIn
documentsIn(a) · a document-term matrix n×p → matrix n×p
Let \(A \in \mathbb{R}^{n \times p}\) be a document-term matrix and \(V = \operatorname{topics}(A)\). Then
\[\operatorname{documentsIn}(A) = A V .\]Each document’s coordinates in topic space, \(AV = U\Sigma\): row \(i\) is document \(i\), column \(j\) its weight on topic \(j\) (scaled by that topic’s singular value).
Intuition Places each document among the topics: row \(i\) gives document \(i\)’s coordinates, and a large magnitude means that topic features strongly in it. Documents on the same theme land close together even when they share few exact words.
Use it when you compare, cluster or search documents in a compact space: keep the first few columns and compare rows by cosine similarity. A new document is placed the same way, by multiplying its term vector by the kept columns of topics. Columns are scaled by the singular values, which gives strong topics more influence on distances.
semantic-searchsimilarity-searchtopic-modelling
def documentsIn(a) ≔ mul(a, rightVectors(a))
Non-negative matrix factorisation
Approximates a non-negative matrix \(V \approx WH\) with non-negative factors, using Lee and Seung’s multiplicative updates, taken alternately (the update of \(H\) uses the already-updated \(W\)):
\[W \leftarrow W \odot \frac{V H^\top}{W H H^\top}, \qquad H \leftarrow H \odot \frac{W^\top V}{W^\top W H}\]Every factor in an update is non-negative, so entries that start positive stay positive with no projection step, and the squared error does not increase from turn to turn. Unlike singular vectors, the parts cannot cancel one another, so they read as parts of the whole. Divisors are guarded at \(10^{-12}\).
Shapes: v is n×p, the starting parts w are n×r and the starting mixing weights h are r×p,
for a chosen rank r. Start from strictly positive values (for example your own uniform draws). The
result depends on the start, and the fixed number of turns t is not a convergence test.
In this section let \(V \in [0, \infty)^{n \times p}\) (the argument v), \(W_0 \in (0, \infty)^{n \times r}\)
(the argument w) and \(H_0 \in (0, \infty)^{r \times p}\) (the argument h), let \(t \in \mathbb{N}\),
and define for \(s = 0, 1, \dots\)
nmfParts
nmfParts(t, v, w, h) · t turns scalar · v data n×p · w starting parts n×r · h starting weights r×p → matrix n×r
Let \(t \in \mathbb{N}\), \(V \in [0, \infty)^{n \times p}\), \(W_0 \in (0, \infty)^{n \times r}\) and \(H_0 \in (0, \infty)^{r \times p}\), with \(W_s, H_s\) as above. Then
\[\operatorname{nmfParts}(t, V, W_0, H_0) = W_t .\]The factor \(W\) after t alternating updates.
Intuition \(V \approx WH\) with nothing negative: each column of v is rebuilt by adding up the r columns of \(W\), in amounts given by the matching column of \(H\). Since nothing can be subtracted, the columns of \(W\) come out as additive parts, such as image fragments or groups of co-occurring items.
Use it when you want the factor attached to the rows of v. With documents in rows and terms in columns, row \(i\) of \(W\) is document \(i\)’s topic mixture; with images as columns, the columns of \(W\) are the part images. Choose a small rank, start from strictly positive random values (a zero stays zero forever), and compare several starts, since results depend on them.
topic-modellingsignal-separationrecommendation
def nmfParts(t, v, w, h) ≔ partsLoop(t, v, w, h)
nmfMix
nmfMix(t, v, w, h) · t turns scalar · v data n×p · w starting parts n×r · h starting weights r×p → matrix r×p
Let \(t\), \(V\), \(W_0\), \(H_0\) and \(W_s, H_s\) be as above. Then
\[\operatorname{nmfMix}(t, V, W_0, H_0) = H_t .\]The factor \(H\) after the same t alternating updates, from the same run as nmfParts.
Intuition The other factor, \(H\): each of its r rows is a non-negative pattern over the columns of v, and each row of v is rebuilt as a sum of those patterns in amounts given by the matching row of \(W\).
Use it when you want the factor attached to the columns of v. With documents in rows and terms in columns, each row of \(H\) is a topic; read its largest entries as key terms. In a user-by-item matrix the rows are item profiles of hidden tastes. Pass the same t, v, w and h as to nmfParts so both factors come from one run. Any positive scale can move between \(W\) and \(H\), so normalise rows of \(H\) before comparing them.
topic-modellingkeyword-extractionrecommendation
def nmfMix(t, v, w, h) ≔ mixLoop(t, v, w, h)
nmfFit
nmfFit(t, v, w, h) · t turns scalar · v data n×p · w starting parts n×r · h starting weights r×p → matrix n×p
Let \(t\), \(V\), \(W_0\), \(H_0\) and \(W_s, H_s\) be as above. Then
\[\operatorname{nmfFit}(t, V, W_0, H_0) = W_t H_t .\]The reconstruction \(WH\) after t turns; compare it with v to judge the fit.
Intuition What the data looks like when rebuilt from only r additive parts. The gap between it and v is what those parts could not explain.
Use it when you judge a factorisation: compute the squared error against v, compare it across ranks to find where adding parts stops helping, or across turn counts to see whether t was enough (nothing else checks convergence). It also gives a smoothed, low-rank version of the matrix, such as predicted scores in a small recommender. Entries missing from v are not skipped: whatever value stands in for them (say 0) is fitted like real data.
model-evaluationhyperparameter-choicerecommendation
def nmfFit(t, v, w, h) ≔ approx(nmfParts(t, v, w, h), nmfMix(t, v, w, h))
Building blocks
approx(w, h)— Let \(W \in \mathbb{R}^{n \times r}\) and \(H \in \mathbb{R}^{r \times p}\). Then \(\operatorname{approx}(W, H) = W H\). The product of the two factors.stepW(v, w, h)— Let \(V \in \mathbb{R}^{n \times p}\), \(W \in \mathbb{R}^{n \times r}\) and \(H \in \mathbb{R}^{r \times p}\). Then \(\operatorname{stepW}(V, W, H)_{ij} = W_{ij} \dfrac{(V H^{\top})_{ij}}{\max\big((W H H^{\top})_{ij}, 10^{-12}\big)}\). One multiplicative update of the parts.stepH(v, w, h)— Let \(V\), \(W\) and \(H\) be as forstepW. Then \(\operatorname{stepH}(V, W, H)_{ij} = H_{ij} \dfrac{(W^{\top} V)_{ij}}{\max\big((W^{\top} W H)_{ij}, 10^{-12}\big)}\). One multiplicative update of the mixing weights.nextW(v, w, h)— Let \(V\), \(W\) and \(H\) be as forstepW. Then \(\operatorname{nextW}(V, W, H) = \operatorname{stepW}(V, W, H)\). The parts after one alternating turn.nextH(v, w, h)— Let \(V\), \(W\) and \(H\) be as forstepW. Then \(\operatorname{nextH}(V, W, H) = \operatorname{stepH}\big(V, \operatorname{nextW}(V, W, H), H\big)\). The weights after one alternating turn, stepped from the updated parts.partsLoop(t, v, w, h)— Let \(t \in \mathbb{N}\) and \(V\), \(W\), \(H\) be as forstepW. Then \(\operatorname{partsLoop}(0, V, W, H) = W\) and \(\operatorname{partsLoop}(t, V, W, H) = \operatorname{partsLoop}\big(t - 1, V, \operatorname{nextW}(V, W, H), \operatorname{nextH}(V, W, H)\big)\). Runstturns and answers the parts.mixLoop(t, v, w, h)— Let \(t \in \mathbb{N}\) and \(V\), \(W\), \(H\) be as forstepW. Then \(\operatorname{mixLoop}(0, V, W, H) = H\) and \(\operatorname{mixLoop}(t, V, W, H) = \operatorname{mixLoop}\big(t - 1, V, \operatorname{nextW}(V, W, H), \operatorname{nextH}(V, W, H)\big)\). Runs the sametturns and answers the weights.
Random projection
The Johnson–Lindenstrauss lemma says that multiplying by a random matrix nearly preserves pairwise distances, whatever the data. It is the one reduction that looks at nothing. You supply the random matrix, typically with independent standard normal entries.
projected
projected(x, r) · x data n×p · r random matrix p×k → matrix n×k
Let \(X \in \mathbb{R}^{n \times p}\) be the data matrix and \(R \in \mathbb{R}^{p \times k}\) (the
argument r) a matrix of draws, for example standard normal draws. Then
The points projected to k dimensions, \(X R / \sqrt{k}\). The \(1/\sqrt{k}\) keeps distances on
the original scale when r has unit-variance entries.
Intuition A random linear map to k dimensions nearly preserves every pairwise distance with high probability, and the k needed grows only with the logarithm of the number of points, not with the original dimension p.
Use it when p is huge (hashed text features, genomic data, large embeddings) and you want a fast reduction that needs no fitting before nearest-neighbour search, clustering or kernel methods. To keep distortion within a factor \(1 \pm \varepsilon\), k must be on the order of \(\ln n / \varepsilon^2\), often hundreds. Fill r with independent standard normal draws; the \(1/\sqrt{k}\) keeps scale only for unit-variance entries. PCA or topics preserve more at a given k, at far higher cost.
dimensionality-reductionhigh-dimensional-datasimilarity-search
def projected(x, r) ≔ einsum("np,pk->nk", x, r) / sqrt(columns(r))
Stochastic neighbour embedding
t-SNE matches neighbourhoods expressed as probabilities. In the original space, a Gaussian affinity gives the joint probabilities \(p_{ij}\); in the low-dimensional layout, a heavy-tailed Student-t kernel gives
\[q_{ij} = \frac{(1 + \lVert y_i - y_j \rVert^2)^{-1}}{\sum_{k \ne l} (1 + \lVert y_k - y_l \rVert^2)^{-1}}\]and gradient descent moves the layout to lower the Kullback–Leibler divergence \(\mathrm{KL}(P \,\Vert\, Q)\), whose gradient is
\[\frac{\partial C}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij}) (1 + \lVert y_i - y_j \rVert^2)^{-1} (y_i - y_j)\]The heavy tail is what lets distant points sit far apart cheaply, so clusters separate instead of
crowding together. This version uses one bandwidth g for every point (no perplexity calibration),
plain gradient steps (no momentum or early exaggeration), and a starting layout you supply.
In this section let \(X \in \mathbb{R}^{n \times p}\) be the data matrix with rows \(x_i^{\top}\), \(g > 0\) the bandwidth, and \(Y \in \mathbb{R}^{n \times d}\) a layout with rows \(y_i^{\top}\).
nearProbs
nearProbs(x, g) · x data n×p · g bandwidth scalar → matrix n×n
Let \(X \in \mathbb{R}^{n \times p}\) and \(g > 0\), and let \(a_{ij} = \operatorname{rawNear}(X, g)_{ij}\). Then
\[\operatorname{nearProbs}(X, g)_{ij} = p_{ij} = \frac{a_{ij}}{\max\big(\sum_{k,l} a_{kl},\, 10^{-300}\big)}, \qquad a_{ij} = [\, i \ne j \,] \exp\!\big(-g \lVert x_i - x_j \rVert^2\big) .\]The joint neighbour probabilities \(p_{ij} \propto \exp(-g \lVert x_i - x_j \rVert^2)\) with a zero diagonal, normalised so the whole matrix sums to one. Unlike standard t-SNE, the matrix is normalised once as a whole rather than per point and then symmetrised.
Intuition Turns distances in the original space into a probability distribution over pairs: close pairs get most of the mass, distant pairs almost none. This is the neighbourhood structure t-SNE tries to reproduce in the layout.
Use it when you compute the target p for sneGrad; embedded builds it for you. g decides what counts as near. With one global g, points in sparse regions can end up with almost no neighbour mass, which real t-SNE avoids by tuning a width per point. Standardise features, and pick g so each point’s few nearest neighbours keep noticeable affinity, for example \(1/g\) near typical nearest-neighbour squared distances.
visualisationdimensionality-reduction
def nearProbs(x, g) ≔ rawNear(x, g) / max(einsum("ij->", rawNear(x, g)), 1e-300)
farProbs
farProbs(y) · y layout n×d → matrix n×n
Let \(Y \in \mathbb{R}^{n \times d}\) be a layout and \(b_{ij} = \operatorname{rawFar}(Y)_{ij}\). Then
\[\operatorname{farProbs}(Y)_{ij} = q_{ij} = \frac{b_{ij}}{\max\big(\sum_{k,l} b_{kl},\, 10^{-300}\big)}, \qquad b_{ij} = \frac{[\, i \ne j \,]}{1 + \lVert y_i - y_j \rVert^2} .\]The Student-t probabilities \(q_{ij}\) of the layout, zero on the diagonal, summing to one.
Intuition The same kind of pair distribution, measured in the layout with a Student-t kernel. Its heavy tail leaves moderately distant layout points some probability, so pushing them further apart costs little, and that is what lets clusters separate on the page.
Use it when you monitor an embedding, comparing it with nearProbs (for example by their KL divergence) to see how well the layout keeps neighbourhoods, or write your own t-SNE variant. It has no width parameter: the layout’s scale is learnt instead.
visualisationmodel-evaluation
def farProbs(y) ≔ rawFar(y) / max(einsum("ij->", rawFar(y)), 1e-300)
sneGrad
sneGrad(p, y) · p neighbour probabilities n×n · y layout n×d → matrix n×d
Let \(P \in \mathbb{R}^{n \times n}\) (the argument p) be neighbour probabilities and
\(Y \in \mathbb{R}^{n \times d}\) a layout, and let \(q = \operatorname{farProbs}(Y)\). Then
The gradient of \(\mathrm{KL}(P \,\Vert\, Q)\) with respect to every row of the layout.
Intuition A force on every layout point. Pairs that are too far apart in the layout (\(p_{ij} > q_{ij}\)) attract and pairs that are too close repel, each weighted by the Student-t factor. A descent step moves points toward their true neighbours and away from impostors.
Use it when you write your own optimiser around t-SNE, for instance adding momentum (momentum and heavyStep in Prelude.Learning) or early exaggeration (multiplying p by a constant for the first turns), neither of which embedded does. It builds an n×n×d tensor, so each step costs on the order of \(n^2\).
visualisationgradient-based-training
def sneGrad(p, y) ≔ 4 · einsum("ij,ijd->id", (p − farProbs(y)) · farWeight(y), spreadY(y))
embedded
embedded(t, r, x, g, y0) · t steps scalar · r learning rate scalar · x data n×p · g bandwidth scalar · y0 starting layout n×d → layout n×d
Let \(t \in \mathbb{N}\), let \(\eta > 0\) (the argument r) be the learning rate, \(X \in \mathbb{R}^{n \times p}\),
\(g > 0\) and \(Y_0 \in \mathbb{R}^{n \times d}\) (the argument y0), and let
\(P = \operatorname{nearProbs}(X, g)\). Let \(Y_{s+1} = Y_s - \eta \operatorname{sneGrad}(P, Y_s)\) for
\(s = 0, 1, \dots\). Then
Runs t gradient steps of size r from y0 and returns the layout. Start from a small random
layout of your own; the answer depends on it.
Intuition Gradient descent that arranges points in a few dimensions, usually two, so their neighbourhoods match those of the original data. Local groupings survive; distances between clusters and cluster sizes in the picture mean little.
Use it when you want to see clusters in high-dimensional data: embeddings, cell profiles, image features. Start from a small random layout, never from identical points, whose gradient is exactly zero so they never move. The probabilities sum to one over all pairs, so gradients are tiny and the rate must be large by ordinary standards. Compare a few seeds. Do not read global distances from the result or use it as model features; mdsCoords or isomap keep global geometry.
visualisationclusteringdimensionality-reduction
def embedded(t, r, x, g, y0) ≔ sneLoop(t, r, nearProbs(x, g), y0)
Building blocks
notSelf(a)— Let \(A \in \mathbb{R}^{n \times n}\). Then \(\operatorname{notSelf}(A)_{ij} = 1 - [\, i = j \,]\). A matrix of ones with a zero diagonal, the shape ofa.rawNear(x, g)— Let \(X \in \mathbb{R}^{n \times p}\) and \(g > 0\). Then \(\operatorname{rawNear}(X, g)_{ij} = [\, i \ne j \,] \exp\big(-g \lVert x_i - x_j \rVert^2\big)\). Gaussian affinities with the diagonal removed, before normalisation.rawFar(y)— Let \(Y \in \mathbb{R}^{n \times d}\). Then \(\operatorname{rawFar}(Y)_{ij} = \dfrac{[\, i \ne j \,]}{1 + \lVert y_i - y_j \rVert^2}\). Student-t affinities of the layout with the diagonal removed, before normalisation.farWeight(y)— Let \(Y \in \mathbb{R}^{n \times d}\). Then \(\operatorname{farWeight}(Y)_{ij} = \dfrac{[\, i \ne j \,]}{1 + \lVert y_i - y_j \rVert^2}\), the same matrix as \(\operatorname{rawFar}(Y)\). The Student-t factor that weights each pair in the gradient.spreadY(y)— Let \(Y \in \mathbb{R}^{n \times d}\). Then \(\operatorname{spreadY}(Y) \in \mathbb{R}^{n \times n \times d}\) with \(\operatorname{spreadY}(Y)_{ijc} = Y_{ic} - Y_{jc}\). The n×n×d tensor of differences between layout rows.sneLoop(t, r, p, y)— Let \(t \in \mathbb{N}\), \(\eta > 0\) (the argumentr), \(P \in \mathbb{R}^{n \times n}\) and \(Y \in \mathbb{R}^{n \times d}\). Then \(\operatorname{sneLoop}(0, \eta, P, Y) = Y\) and \(\operatorname{sneLoop}(t, \eta, P, Y) = \operatorname{sneLoop}\big(t - 1, \eta, P, Y - \eta \operatorname{sneGrad}(P, Y)\big)\). Takestgradient steps of sizer; used byembedded.
Examples
Every example below is one file, Shapes.hb, so the C# class is Shapes.
Place points on a line from their distances
Three points at 0, 1 and 3 on a line, given only their distance matrix.
open Prelude.Manifold
fn scaledCoords(d) ↦ mdsCoords(d)
fn scaledStrength(d) ↦ mdsStrength(d)
using Examples;
var d = new double[,] { { 0, 1, 3 }, { 1, 0, 2 }, { 3, 2, 0 } };
var coords = Shapes.ScaledCoords(d);
var strength = Shapes.ScaledStrength(d);
Output (verified, rounded to 6 significant digits, negative zeros shown as 0):
mdsCoords = [[1.33333, 0, 0], [0.333333, 0, 0], [-1.66667, 0, 0]]
mdsStrength = [4.66667, 0, 0]
The first column is the line centred on its mean 4/3, with its sign flipped (axes have no preferred sign). Only one eigenvalue is non-zero, because the points are one-dimensional.
Unroll a U-shaped path with Isomap
Seven points walking up, across and down a U. The two ends are 2 apart in a straight line but 6 apart along the path. Isomap is called as two doorways.
open Prelude.Manifold
fn straightLine(x) ↦ distance(x, x)
fn alongManifold(x, k) ↦ geodesics(x, k)
fn scaledCoords(d) ↦ mdsCoords(d)
using Examples;
var u = new double[,] { { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 2 }, { 2, 2 }, { 2, 1 }, { 2, 0 } };
var straight = Shapes.StraightLine(u);
var along = Shapes.AlongManifold(u, 1.0); // one neighbour per point
var unrolled = Shapes.ScaledCoords(along);
Output (verified, first row of each distance matrix and first column of the coordinates, rounded):
straight row 0 = [0, 1, 2, 2.23607, 2.82843, 2.23607, 2]
along row 0 = [0, 1, 2, 3, 4, 5, 6]
unrolled col 0 = [3, 2, 1, 0, -1, -2, -3]
Every other column of the coordinates is zero: measured along the path, the U is a straight segment of length 6.
Read topics from a document-term matrix
Three documents over two terms.
open Prelude.Manifold
fn topicStrength(a) ↦ topicWeights(a)
fn topicAxes(a) ↦ topics(a)
fn placedDocuments(a) ↦ documentsIn(a)
using Examples;
var a = new double[,] { { 2, 0 }, { 0, 1 }, { 0, 0 } };
var weights = Shapes.TopicStrength(a);
var axes = Shapes.TopicAxes(a);
var docs = Shapes.PlacedDocuments(a);
Output (verified):
topicWeights = [2, 1]
topics = [[1, 0], [0, 1]]
documentsIn = [[2, 0], [0, 1], [0, 0]]
Factor a rank-one matrix into non-negative parts
A 3×2 matrix whose rows are multiples of [1, 2], factored at rank 1 from a start of all ones, for
200 turns.
open Prelude.Manifold
fn partsOf(t, v, w, h) ↦ nmfParts(t, v, w, h)
fn mixOf(t, v, w, h) ↦ nmfMix(t, v, w, h)
fn rebuilt(t, v, w, h) ↦ nmfFit(t, v, w, h)
using Examples;
var v = new double[,] { { 1, 2 }, { 2, 4 }, { 3, 6 } };
var w0 = new double[,] { { 1 }, { 1 }, { 1 } };
var h0 = new double[,] { { 1, 1 } };
var w = Shapes.PartsOf(200.0, v, w0, h0);
var h = Shapes.MixOf(200.0, v, w0, h0);
var fit = Shapes.Rebuilt(200.0, v, w0, h0);
Output (verified, rounded to 6 significant digits):
nmfParts = [[1.5], [3], [4.5]]
nmfMix = [[0.666667, 1.33333]]
nmfFit = [[1, 2], [2, 4], [3, 6]]
The split between \(W\) and \(H\) is not unique (any positive scale moves between them), but their
product reproduces v.
Project points with a supplied random matrix
Two points in four dimensions projected to two with a fixed ±1 matrix (in real use, draw it at random).
open Prelude.Manifold
fn squeezed(x, r) ↦ projected(x, r)
using Examples;
var x = new double[,] { { 1, 2, 3, 4 }, { 0, 0, 1, 1 } };
var r = new double[,] { { 1, 0 }, { 1, 0 }, { 0, 1 }, { 0, -1 } };
var low = Shapes.Squeezed(x, r);
Output (verified, rounded to 6 significant digits):
projected = [[2.12132, -0.707107], [0, 0]]
The first point maps to \((1 + 2, \; 3 - 4)/\sqrt{2}\).