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 and 2.11. Kernel PCA, 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, 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).

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.

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 instead of failing later at fit time.

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:

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:

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.

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.

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.

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 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 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.

// 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) 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 instead.

A common and valid workflow runs the reverse order. Cluster in the original space, for example with KMeans. 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.

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, 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 type:

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.