# 2.12. t-SNE
`TSNE` is RustyML's implementation of t-distributed Stochastic Neighbor Embedding (t-SNE). It is a common method to turn high-dimensional data into a 2-D or 3-D plot. `TSNE` lives in `machine_learning::manifold`. Unlike the linear reducers in [2.10. Principal Component Analysis](./2.10._Principal_Component_Analysis.md) and [2.11. Kernel PCA](./2.11._Kernel_PCA.md), it learns no reusable projection. It embeds only the points you give it.
If you know scikit-learn, this is `sklearn.manifold.TSNE`. The API mirrors it on purpose. It exposes a `perplexity`, a `learning_rate`, an `n_iter`, a choice between an exact and a Barnes-Hut gradient, and PCA or random initialization.
## 2.12.1. What t-SNE is for, and what it is not
t-SNE is a *visualization* tool. Its only job is to place similar high-dimensional points near each other on a low-dimensional canvas. This lets a human eye see the cluster structure in the data. Most misuse of t-SNE comes from treating it as a general-purpose dimensionality reducer. It is not one.
The design makes this restriction concrete. In [the estimator traits](../Chapter-01/1.5._The_Prelude_and_Imports.md), PCA and Kernel PCA implement both `Transform` (project new, unseen data through a fitted model) and `FitTransform`. `TSNE` implements only `FitTransform`. It has no `transform` method, no stored fitted state, and no way to add a point to an existing embedding. The embedding coordinates are the optimization variables themselves, not the output of a learned function.
Use PCA to project new samples later. Use t-SNE to plot the data you already have. The 2 tools complement each other. A common pipeline runs PCA first and feeds the result to t-SNE (see [2.12.9](#2129-cost-scaling-and-preprocessing-wide-data)).
The inherent method takes `&self` and returns a fresh embedding, so it never mutates the model. You can call it directly, or through the `FitTransform` trait. The trait method takes `&mut self`. It just forwards to the inherent method.
## 2.12.2. How t-SNE works
t-SNE measures similarity twice. It measures once in the original space and once in the embedding. Then it moves the embedding until the two agree.
In the high-dimensional space, t-SNE builds a set of **pairwise affinities**. For each point `i`, it centers a Gaussian on that point. It converts squared distances to the neighbors of `i` into conditional probabilities `p_{j|i}`. Near points get a high probability. Far points get a probability close to zero.
The width of the Gaussian is not fixed. RustyML solves it *per point*, with a binary search over `sigma`. The search matches the entropy of each point's neighbor distribution to a target you set. That target is the **perplexity**. Perplexity reads as an effective neighbor count, roughly how many neighbors each point should feel.
RustyML runs the search for up to 50 bisection steps, to a tolerance of `1e-5`. The self-distance always maps to exactly zero, so a point is never its own neighbor. The conditional probabilities are then symmetrized into joint probabilities `p_ij`. These joint probabilities sum to 1 over all pairs.
In the embedding, t-SNE uses a different, heavier-tailed kernel. This kernel is the **Student-t distribution with 1 degree of freedom**, where `q_ij` is proportional to `(1 + ||y_i - y_j||^2)^-1`. This is the trick that gives the "t" its name, and it fixes the **crowding problem**. A Gaussian in 2-D cannot make room for all the points that were moderately distant neighbors in, say, 50-D.
There is simply not enough area, and everything collapses into a blob. The t-distribution's fat tail lets moderately similar points sit far apart in the map, without paying much probability cost. Clusters separate cleanly instead of crushing together.
The 2 distributions are matched by minimizing the **Kullback-Leibler divergence** `KL(P || Q)`, with gradient descent on the embedding coordinates. KL divergence is asymmetric on purpose. It heavily penalizes placing a high-`p` (truly near) pair far apart in the map. It barely penalizes placing a low-`p` (truly far) pair close together. That asymmetry is why t-SNE preserves *local* neighborhoods faithfully, and treats *global* distances as expendable. It also drives all the plot-reading caveats in [2.12.8](#2128-reading-a-t-sne-plot).
## 2.12.3. Constructing a model
`TSNE::new` takes the 4 core hyperparameters. It validates them up front. On a bad value, it returns [`Error::InvalidParameter`](../Chapter-01/1.6._Error_Handling.md) instead of failing later at fit time.
```rust
use ndarray::array;
use rustyml::machine_learning::manifold::t_sne::{TSNE, TSNEMethod};
fn main() {
// Two loose groups of points in 3-D feature space.
let x = array![
[0.0, 0.0, 0.0],
[0.2, 0.1, -0.1],
[-0.1, 0.2, 0.1],
[0.1, -0.2, 0.0],
[5.0, 5.0, 5.0],
[5.2, 4.9, 5.1],
[4.8, 5.1, 4.9],
[5.1, 5.0, 4.8],
];
// new(n_components, perplexity, learning_rate, n_iter) -> Result<TSNE, Error>.
let tsne = TSNE::new(2, 3.0, 200.0, 300)
.unwrap()
.with_method(TSNEMethod::Exact)
.unwrap();
// fit_transform takes &self and returns an (n_samples, n_components) array.
let embedding = tsne.fit_transform(&x).unwrap();
assert_eq!(embedding.shape(), &[8, 2]);
println!("embedding shape: {:?}", embedding.shape());
}
```
The constructor's validation rules are narrow but strict:
| Parameter | Type | Constraint | On violation |
| --- | --- | --- | --- |
| `n_components` | `usize` | greater than 0 (use 2, or 3 for a rotatable plot) | `InvalidParameter` |
| `perplexity` | `f64` | strictly positive and finite | `InvalidParameter` |
| `learning_rate` | `f64` | strictly positive and finite | `InvalidParameter` |
| `n_iter` | `usize` | greater than 0 | `InvalidParameter` |
`TSNE::default()` gives `new(2, 30.0, 200.0, 1000)`, with PCA initialization and Barnes-Hut. This is a reasonable starting point for a real dataset, not a toy-sized one. The builder methods set everything else. Each builder returns `Self` for chaining, except `with_method`. `with_method` returns `Result`, because Barnes-Hut must validate its angle and its dimensionality:
| Builder | Sets | Default |
| --- | --- | --- |
| `with_method(TSNEMethod)` returns `Result` | exact or Barnes-Hut gradient | Barnes-Hut when `n_components` is 3 or less, else Exact |
| `with_init(Init)` | `Init::PCA` or `Init::Random` | `Init::PCA` |
| `with_random_state(u64)` | seed for the random-init path | `None` |
| `with_min_grad_norm(f64)` | early-stopping gradient threshold | `1e-7` |
Every stored field has a matching getter: `get_n_components`, `get_perplexity`, `get_learning_rate`, `get_n_iter`, `get_random_state`, `get_init`, `get_method`, and `get_min_grad_norm`. Use these to check how a chain of builder calls resolved.
## 2.12.4. Exact versus Barnes-Hut
`TSNEMethod` sets the cost tradeoff. RustyML's default is **not** the exact method. This differs from what many users assume.
`TSNEMethod::Exact` builds the full dense `n x n` joint-probability matrix. It computes the gradient over every pair on each iteration. This costs `O(n^2)` per step, in both time and memory. `TSNEMethod::Exact` supports any `n_components`.
`TSNEMethod::BarnesHut { angle }` keeps the affinities sparse. Each point talks only to its `k = ceil(3 * perplexity) + 1` nearest neighbors. It summarizes the repulsive forces with a space-partitioning tree. This gives roughly `O(n log n)` per iteration.
The `angle` (theta) is in `[0, 1)` and trades accuracy for speed. A larger angle opens tree cells sooner, so it runs faster but coarser. The value `0.5` is the standard balance. The tree lives in the embedding space, so Barnes-Hut supports only `n_components` of 3 or less.
`TSNE::new` picks Barnes-Hut with `angle = 0.5` whenever `n_components` is 3 or less. This covers every visualization case. Otherwise it falls back to Exact. `with_method` enforces the same constraints. It rejects an `angle` outside `[0, 1)` with `InvalidParameter`. It also rejects Barnes-Hut paired with more than 3 components, with the same error.
```rust
use ndarray::array;
use rustyml::machine_learning::{TSNE, TSNEMethod};
fn main() {
let x = array![
[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0],
[8.0, 8.0], [9.0, 8.0], [8.0, 9.0], [9.0, 9.0],
];
// Default for 2 components: Barnes-Hut with angle 0.5.
let bh = TSNE::new(2, 3.0, 200.0, 300).unwrap();
assert_eq!(bh.get_method(), TSNEMethod::BarnesHut { angle: 0.5 });
// Coarser, faster tree.
let coarse = TSNE::new(2, 3.0, 200.0, 300)
.unwrap()
.with_method(TSNEMethod::BarnesHut { angle: 0.8 })
.unwrap();
// Exact O(n^2) gradient, the only option once n_components > 3.
let exact = TSNE::new(2, 3.0, 200.0, 300)
.unwrap()
.with_method(TSNEMethod::Exact)
.unwrap();
for model in [bh, coarse, exact] {
let emb = model.fit_transform(&x).unwrap();
assert_eq!(emb.ncols(), 2);
}
}
```
The Barnes-Hut gradient is scaled to match the exact path. The factor of 4 is folded in. The same `learning_rate` fits both methods, so you can switch between them without retuning. Both methods are bit-reproducible. The tree build is deterministic. The normalizer is summed in a fixed order, so the result does not depend on thread scheduling.
## 2.12.5. Perplexity and the sample-count rules
Perplexity is the hyperparameter that most changes the shape of your plot. It sets the effective neighbor count that each point calibrates its Gaussian to. A low perplexity emphasizes local structure and fractures the data into many small islands. A high perplexity blends larger neighborhoods and can smear distinct clusters together.
Values from 5 to 50 cover almost every use. The right value rises as your dataset grows.
There are 2 limits on perplexity, and they differ. `fit_transform`, not `TSNE::new`, enforces the first one. `perplexity` must be strictly less than the number of samples, or `fit_transform` returns `InvalidParameter`.
The second limit is about Barnes-Hut sparsity, not correctness. The Barnes-Hut path keeps `k = ceil(3 * perplexity) + 1` neighbors per point, capped at `n - 1`. Some t-SNE implementations reject a perplexity above roughly `n / 3` as an error. RustyML does not. Past that point, the neighbor cap has already saturated at `n - 1`, so every point already lists every other point as a neighbor. The embedding stays correct, but Barnes-Hut then loses its usual speed advantage over Exact for that data.
The enforced `perplexity < n` check is the one that matters, and it is the same for both methods. This is also why t-SNE on a handful of points gives little insight. With 20 samples, the enforced rule already caps perplexity below 20. A run near that ceiling has little real neighbor structure to find, and the clusters in the plot are mostly noise.
```rust
use ndarray::array;
use rustyml::error::Error;
use rustyml::machine_learning::{TSNE, TSNEMethod};
fn main() {
let x = array![[0.0, 0.0], [1.0, 1.0], [2.0, 0.5]]; // 3 samples
// perplexity must be < n_samples.
// Here 3.0 is not < 3.
let model = TSNE::new(2, 3.0, 200.0, 100)
.unwrap()
.with_method(TSNEMethod::Exact)
.unwrap();
match model.fit_transform(&x) {
Err(Error::InvalidParameter { .. }) => {
println!("perplexity too large for this sample count");
}
other => panic!("expected InvalidParameter, got {other:?}"),
}
}
```
## 2.12.6. Initialization, seeding, and reproducibility
`Init` chooses the starting point for optimization. This choice affects reproducibility.
`Init::PCA` is the default. It starts the embedding from the top principal components of the input, rescaled to a small spread. When PCA succeeds, this method is **deterministic** and ignores `random_state`. The same data always gives the same starting layout, and so the same final embedding. It also gives the optimizer a well-spread starting layout. This is why `Init::PCA` is the default.
PCA init can fail in 2 cases. The input has fewer features than `n_components`, or its leading component is degenerate (zero spread). When either happens, `Init::PCA` falls back to the same random init that `Init::Random` uses. That fallback does consult `random_state` (or the global seed), so determinism returns with it.
`Init::Random` starts from tiny random noise. It seeds this noise through the crate's central RNG. Apart from the PCA fallback above, this is the only path that consults `random_state`. Set a seed with `with_random_state`, and the run becomes bit-reproducible. Leave it unset, and the start comes from entropy, so every run differs.
```rust
use ndarray::array;
use rustyml::machine_learning::{Init, TSNE, TSNEMethod};
use rustyml::set_global_seed;
fn main() {
let x = array![
[0.0, 0.0], [1.0, 0.5], [0.5, 1.0],
[6.0, 6.0], [7.0, 6.5], [6.5, 7.0],
];
// Explicit per-model seed: reproducible regardless of any global seed.
let seeded = TSNE::new(2, 2.0, 200.0, 250)
.unwrap()
.with_init(Init::Random)
.with_random_state(42)
.with_method(TSNEMethod::Exact)
.unwrap();
let e1 = seeded.fit_transform(&x).unwrap();
let e2 = seeded.fit_transform(&x).unwrap();
assert_eq!(e1, e2); // bit-identical
// An unseeded random-init model becomes reproducible only because
// this code sets the thread-local global seed first.
set_global_seed(7);
let unseeded = TSNE::new(2, 2.0, 200.0, 250)
.unwrap()
.with_init(Init::Random)
.with_method(TSNEMethod::Exact)
.unwrap();
let _emb = unseeded.fit_transform(&x).unwrap();
}
```
The `random_state` field feeds `make_rng`. Every randomized component in the crate resolves its seed through this same path. An explicit `Some(seed)` is independent and never touches the global stream. A `None` seed defers to whatever [`set_global_seed`](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md) set on the current thread.
This has 2 practical effects. With the default PCA init, the seed does not matter. Do not expect `with_random_state` to change anything, unless you also switch to `Init::Random`. The parallel reductions in the optimizer are also order-stable. Because of this, a fixed configuration reproduces bit for bit on a given machine, no matter the thread count. Chapter [7.1. Reproducibility and Random Seeds](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md) covers the full mechanics, including the thread-local nature of the global seed.
## 2.12.7. Inside the optimizer
The gradient descent loop is not plain SGD. Its phases explain most of what `n_iter` and `learning_rate` do.
The first 250 iterations run under **early exaggeration** (or all iterations, if `n_iter` is smaller). During this phase, t-SNE multiplies every joint probability `p_ij` by 12. This inflates the attractive forces. Tight clusters pull together and carve out empty space between groups before the layout settles. When the phase ends, the exaggeration factor drops back to 1, and the map fine-tunes.
If `n_iter` is too low, the run may stop inside the exaggeration phase. This leaves an over-contracted embedding that has not yet fine-tuned. This is one reason the defaults use 1000 iterations. Real runs rarely want fewer than a few hundred iterations.
Momentum follows the same schedule. It is `0.5` during exaggeration and `0.8` afterward. This gives the refinement phase more inertia through the flat parts of the KL landscape. Each coordinate also has an **adaptive gain** (Jacobs' delta-bar-delta). This gain grows when the gradient keeps its sign, and decays when the step oscillates. The gain is floored at `0.01`.
So `learning_rate` sets only a base step size. The effective per-parameter rate adapts as the run proceeds. After every step, the embedding is re-centered on the origin, to stop it from drifting. This is why the output columns always come back with a mean near zero.
`min_grad_norm` controls early stopping. After the exaggeration phase ends, the loop checks the largest absolute gradient entry on each iteration. It stops as soon as that value drops below the threshold (the default is `1e-7`). This saves iterations once the map has converged. The loop skips this check during exaggeration, because the inflated gradients there would never trip it. Set `min_grad_norm` to `0.0` to disable early stopping and always run the full `n_iter`.
```rust,ignore
// The full builder surface, for reference.
let tsne = TSNE::new(2, 30.0, 200.0, 1000)? // n_components, perplexity, lr, n_iter
.with_init(Init::PCA) // or Init::Random
.with_random_state(0) // only affects Init::Random
.with_method(TSNEMethod::Exact)? // or BarnesHut { angle }
.with_min_grad_norm(0.0); // disable early stopping
```
Enable the `show_progress` feature ([1.2. Installation and Feature Flags](../Chapter-01/1.2._Installation_and_Feature_Flags.md)) to make the loop compute and show the live KL divergence on each iteration. Without the feature, the loop skips that pass entirely. It costs nothing in a production build.
## 2.12.8. Reading a t-SNE plot
t-SNE optimizes local neighborhoods, and it discards global geometry. Because of this, the plot misrepresents several things. Read it with that in mind.
The distance between 2 clusters means nothing. 2 blobs that land far apart on the canvas are not more different than 2 blobs that land close together. The gap between clusters is an artifact of the KL objective and the optimizer's random walk. It is not a measurement.
The same holds for cluster *size*. How spread out a cluster's points look is not a measure of its true variance. Dense regions get inflated, and sparse regions get compressed, to keep local neighborhoods intact. Never treat a between-cluster distance or a cluster diameter from a t-SNE map as a quantitative result.
The embedding depends on initialization, and the optimization is non-convex. So a single run is only 1 sample from a distribution of possible layouts. Run several. Vary the seed under `Init::Random`, or sweep perplexity across values such as 5, 30, and 50. Trust only the structure that survives across runs. A cluster that appears at one perplexity and dissolves at another is a likely artifact, not a finding.
Do not feed t-SNE coordinates into a downstream model as features. They are a picture, not a reusable representation. There is no `transform`, so they do not apply to new data. They also change across runs, and carry no meaning as distances. For reusable, projectable features, use [PCA](./2.10._Principal_Component_Analysis.md) instead.
A common and valid workflow runs the reverse order. Cluster in the original space, for example with [KMeans](./2.7._KMeans_Clustering.md). Then color a t-SNE plot by the cluster label, to check the clustering by eye.
## 2.12.9. Cost, scaling, and preprocessing wide data
Both methods pay an `O(n^2)` setup cost, once, up front. This setup computes pairwise distances and calibrates the per-point sigmas. Barnes-Hut still runs an `O(n^2)` neighbor search before its cheaper `O(n log n)` iterations start. So the per-iteration cost differs, `O(n^2)` for Exact against `O(n log n)` for Barnes-Hut. Neither method escapes the quadratic setup term. The exact path also needs a full `n x n` matrix in memory.
In round numbers, exact t-SNE stays comfortable into the low thousands of points. Beyond that, Barnes-Hut lets you push toward tens of thousands of points. Let the default method choice handle this.
The other lever is the *width* of your data. The distance computation costs `O(n^2 * d)`, for `d` input features. High-dimensional Euclidean distances are also noisy. The standard remedy reduces wide data to about 50 dimensions with PCA first, then runs t-SNE on that. This helps both speed and quality.
This step differs from `Init::PCA`. `Init::PCA` only uses the top 2 or 3 components, to seed the *starting layout*. Here, PCA compresses the *input itself*, before t-SNE ever runs.
```rust
use ndarray::array;
use rustyml::machine_learning::{PCA, TSNE, TSNEMethod};
fn main() {
// 10 samples, 8 raw features (imagine hundreds in a real problem).
let x = array![
[0.10, 0.21, 0.05, 0.30, 0.11, 0.02, 0.22, 0.14],
[0.12, 0.19, 0.08, 0.28, 0.09, 0.05, 0.20, 0.10],
[0.09, 0.23, 0.03, 0.31, 0.13, 0.01, 0.24, 0.16],
[0.11, 0.20, 0.06, 0.29, 0.10, 0.03, 0.21, 0.12],
[0.80, 0.70, 0.90, 0.10, 0.60, 0.85, 0.15, 0.75],
[0.82, 0.68, 0.88, 0.12, 0.62, 0.83, 0.17, 0.73],
[0.78, 0.72, 0.91, 0.09, 0.58, 0.87, 0.14, 0.77],
[0.81, 0.69, 0.89, 0.11, 0.61, 0.84, 0.16, 0.74],
[0.40, 0.45, 0.42, 0.55, 0.38, 0.47, 0.52, 0.44],
[0.42, 0.43, 0.44, 0.53, 0.36, 0.49, 0.50, 0.46],
];
// Step 1: compress with PCA (here to 4 dims, about 50 for wide data).
let mut pca = PCA::new(4).unwrap();
let reduced = pca.fit_transform(&x).unwrap();
// Step 2: run t-SNE on the compact representation.
let tsne = TSNE::new(2, 3.0, 200.0, 300)
.unwrap()
.with_method(TSNEMethod::Exact)
.unwrap();
let embedding = tsne.fit_transform(&reduced).unwrap();
assert_eq!(embedding.shape(), &[10, 2]);
}
```
The precompute pass and the per-iteration passes switch to RustyML's parallel primitives once the pairwise work is large enough. The GEMM calls parallelize on their own. The reductions stay in a fixed order, so the result does not depend on thread count. To tune throughput, see [7.3. Performance Tuning and Parallelism](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md), which covers the gating thresholds.
## 2.12.10. Errors from fit_transform
Construction validates the 4 hyperparameters. The rest of the checks run at `fit_transform` time. All of them return the crate-wide [`Error`](../Chapter-01/1.6._Error_Handling.md) type:
| Condition | Error variant |
| --- | --- |
| Zero rows | `EmptyInput` |
| A `NaN` or infinite entry in the input | `NonFinite` |
| Fewer than 2 samples | `InvalidInput` |
| `perplexity` not strictly less than the sample count | `InvalidParameter` |
Every failure is a typed `Error`, not a panic. You can match on the cause and react. Retry with a smaller perplexity, clean non-finite rows, or stop, as the example above shows.