Manifold Learning: t-SNE, UMAP, and Reading the Plots

A manifold is a space that can be described locally by a small number of coordinates. A rolled sheet, for example, has two surface coordinates even though it sits in three-dimensional space. Two points can be close through the air while far apart along the sheet. Manifold learning uses relations among sampled points to build lower-dimensional representations. A curved surface may require a nonlinear map to recover its intrinsic distances; curvature alone does not rule out every useful linear projection. The linear baseline is developed in PCA and Matrix Factorization. Here we explain t-SNE and UMAP, then test what their layouts reveal about neighborhoods, distances, and downstream predictions.

The family

These are objectives and modeling choices, not promises that every listed property survives the fit. Too few output dimensions, a poor input metric, sampling gaps, and incomplete optimization can all matter. Conversely, a method can preserve useful information beyond what its objective explicitly names. The measurements below separate neighborhood agreement, between-group distances, and predictive usefulness.

Both t-SNE and UMAP start by measuring relationships among observations in the input space. The Python examples limit numerical-library threads to one for controlled comparisons; coordinates and derived scores can still vary with library versions and numerical implementations. Feature scaling and the distance metric influence those input relationships before the layout is made. UMAP is available in the separate umap-learn package; its worked example below also shows how a fitted embedding can place new observations.

How t-SNE builds the map

For each observation \(x_i\), t-SNE converts distances to other observations into conditional affinities: \(p_{j\mid i}=\exp[-\lVert x_i-x_j\rVert^2/(2\sigma_i^2)]/\sum_{l\ne i}\exp[-\lVert x_i-x_l\rVert^2/(2\sigma_i^2)]\), with \(p_{i\mid i}=0\). Nearby points receive more weight. The bandwidth \(\sigma_i\) adapts to local density through the perplexity setting. Symmetrizing gives \(p_{ij}=(p_{j\mid i}+p_{i\mid j})/(2n)\), where \(n\) is the number of observations.

In a two-dimensional layout, coordinates \(z_i\) define \(q_{ij}\) proportional to \((1+\lVert z_i-z_j\rVert^2)^{-1}\), normalized over all distinct pairs. This Student-t form has heavier tails than a Gaussian, helping accommodate separation when high-dimensional neighbors must fit into a small map. The algorithm adjusts the coordinates to reduce \(\sum_{i\ne j}p_{ij}\log(p_{ij}/q_{ij})\), a KL divergence. Pairs with large input affinity are costly to separate; the normalization also produces interactions among distant points. This is an affinity-matching objective, not a guarantee of exact neighbor preservation.

Between-group distances can change

import numpy as np
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    rng = np.random.default_rng(0)
    centers = np.array([[0, 0], [3, 0], [30, 0]], float)   # generating center distances: 3, 30, 27
    X = np.vstack([rng.normal(centers[i], 0.5, (300, 2)) for i in range(3)])
    X = np.hstack([X, rng.normal(0, 0.1, (900, 8))])
    y = np.repeat([0, 1, 2], 300)

    def centroid_distances(Z):
        m = np.array([Z[y == k].mean(0) for k in range(3)])
        return np.array([np.linalg.norm(m[0] - m[1]), np.linalg.norm(m[0] - m[2]),
                         np.linalg.norm(m[1] - m[2])])

    for name, Z in (("original", X), ("PCA", PCA(2).fit_transform(X)),
                    ("t-SNE p=10", TSNE(2, perplexity=10, random_state=0,
                                        init="pca").fit_transform(X)),
                    ("t-SNE p=30", TSNE(2, perplexity=30, random_state=0,
                                        init="pca").fit_transform(X))):
        d = centroid_distances(Z)
        print(f"{name:12s} d01,d02,d12 = {np.round(d, 2)}   ratio d02/d01 = {d[1] / d[0]:.2f}")
# original     d01,d02,d12 = [ 3.02 30.09 27.07]   ratio d02/d01 = 9.98
# PCA          d01,d02,d12 = [ 3.02 30.09 27.07]   ratio d02/d01 = 9.98
# t-SNE p=10   d01,d02,d12 = [ 62.35 122.35  63.31]   ratio d02/d01 = 1.96
# t-SNE p=30   d01,d02,d12 = [49.96 71.49 51.74]   ratio d02/d01 = 1.43

The generating center distances from group 0 are 3 and 30, a ratio of ten. The sampled centroids give 9.98, and PCA agrees at the displayed precision because the separation lies almost entirely in the retained directions. t-SNE gives 1.96 and 1.43 for the two perplexities. These layouts substantially compress this particular ratio; PCA does not generally preserve all centroid distances exactly either.

t-SNE matches affinity distributions rather than the original numerical distances. Widely separated groups have weak cross-group input affinities, which leaves their arrangement less constrained by attraction. They still interact through the normalized output similarities and repulsive forces, so they are not free to move anywhere without changing the loss. There is no general calibration converting the distance between two islands into a distance between populations in the input space.

A visible group suggests a neighborhood pattern to investigate; its membership is not automatically a valid cluster assignment. The number of plotted observations can be counted, but the area occupied by a group is not a reliable measure of its original variance or density. Verify claims about group separation or relatedness using a relevant input-space metric, graph, or independent evidence.

A silhouette score depends on the representation

import numpy as np
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    rng = np.random.default_rng(0)
    noise = rng.normal(size=(1000, 50))
    def best_partition(M):
        candidates = []
        for k in range(2, 8):
            labels = KMeans(k, n_init=10, random_state=0).fit_predict(M)
            candidates.append((silhouette_score(M, labels), k, labels))
        return max(candidates, key=lambda result: result[0])
    raw_score, raw_k, _ = best_partition(noise)
    print(f"raw 50-D: best silhouette {raw_score:.4f}, k={raw_k}")
    for name, Z in [("PCA-2", PCA(2).fit_transform(noise)),
                    ("2-D Gaussian", rng.normal(size=(1000, 2)))]:
        value, k, _ = best_partition(Z)
        print(f"{name:12s}: best silhouette {value:.4f}, k={k}")
    for perp in (5, 30, 100):
        Z = TSNE(2, perplexity=perp, random_state=0, init="pca").fit_transform(noise)
        value, k, labels = best_partition(Z)
        print(f"perplexity {perp:3d}: embedding silhouette {value:.4f}, k={k},"
              f" same labels in raw space {silhouette_score(noise, labels):.4f}")
# raw 50-D: best silhouette 0.0185, k=2
# PCA-2       : best silhouette 0.3374, k=6
# 2-D Gaussian: best silhouette 0.3247, k=3
# perplexity   5: embedding silhouette 0.3534, k=7, same labels in raw space -0.0048
# perplexity  30: embedding silhouette 0.3503, k=3, same labels in raw space 0.0079
# perplexity 100: embedding silhouette 0.3440, k=3, same labels in raw space 0.0128

The input is sampled from a single isotropic Gaussian population, with no planted mixture groups. The best raw-space silhouette among \(k=2,\ldots,7\) is 0.0185, close to zero. The t-SNE partitions score about 0.34–0.35 in their own coordinates. However, a two-dimensional PCA projection and even a directly generated two-dimensional Gaussian also score above 0.32. A value in this range does not establish discrete groups.

Silhouette evaluates a partition in a chosen geometry. K-means divides even a continuous Gaussian cloud into regions, and a two-dimensional geometry can give those regions a moderate silhouette. The same labels selected in the t-SNE layouts have much lower scores in the original 50-dimensional data. This experiment establishes a mismatch between evaluations in the two spaces; it does not isolate a single mechanism or prove that every t-SNE layout contains artificial islands. The same noise sample is reused across perplexities so that the input is controlled.

Treat clusters found in an embedding as hypotheses and validate them in an appropriate representation, using stability checks, domain evidence, or an independent task. Clustering an embedding can be useful, especially when a learned representation removes irrelevant variation, but a high silhouette in that same embedding is not sufficient validation. Raw coordinates also need a sensible metric and preprocessing; they are not automatically the right reference for every dataset.

Perplexity

Perplexity is an entropy-based effective neighbor count: \(\operatorname{Perp}(P_i)=\exp[-\sum_jp_{j\mid i}\log p_{j\mid i}]\). If four neighbors each receive probability \(1/4\), it equals four; unequal weights reduce the effective count below the number of nonzero entries. t-SNE adjusts a separate bandwidth \(\sigma_i\) at each point to match the requested perplexity. Perplexity is therefore neither a fixed radius nor a hard count of neighbors. It must be below the sample size, and initialization, learning rate, early exaggeration, iteration count, and distance metric also affect the result.

In the next example, the returned t is the generating roll parameter. We first embed all observations without their labels, then use cross-validation only for the regressor that predicts held-out t values on this fixed map. This is a transductive, held-label diagnostic: the test observations’ features already influenced the embedding. It does not estimate performance on newly arriving points. The final pipeline rows instead fit PCA or Isomap within each training fold and transform that fold’s held-out observations.

import numpy as np
from sklearn.manifold import TSNE, Isomap
from sklearn.decomposition import PCA
from sklearn.datasets import make_swiss_roll
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import cross_val_score
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    X, t = make_swiss_roll(n_samples=1500, noise=0.05, random_state=0)
    score = lambda Z: cross_val_score(KNeighborsRegressor(10), Z, t, cv=5).mean()

    for perp in (2, 5, 15, 50, 200, 500):
        Z = TSNE(2, perplexity=perp, random_state=0, init="pca").fit_transform(X)
        print(f"t-SNE perplexity {perp:4d}   fixed-map held-label R2 for roll parameter {score(Z):7.4f}")
    for nn in (5, 10, 30):
        print(f"Isomap n_neighbors {nn:3d}    kNN R2 {score(Isomap(n_neighbors=nn).fit_transform(X)):7.4f}")
    print(f"PCA                        kNN R2 {score(PCA(2).fit_transform(X)):7.4f}")
    print(f"raw 3-D                    kNN R2 {score(X):7.4f}")
    from sklearn.pipeline import make_pipeline
    for name, reducer in [("PCA", PCA(2)), ("Isomap", Isomap(n_neighbors=10))]:
        pipeline = make_pipeline(reducer, KNeighborsRegressor(10))
        print(f"{name:8s} fold-fitted projection, held-out points R2 "
              f"{cross_val_score(pipeline, X, t, cv=5).mean():.4f}")
# t-SNE perplexity    2   fixed-map held-label R2 for roll parameter  0.9685
# t-SNE perplexity    5   fixed-map held-label R2 for roll parameter  0.9986
# t-SNE perplexity   15   fixed-map held-label R2 for roll parameter  0.9988
# t-SNE perplexity   50   fixed-map held-label R2 for roll parameter  0.9993
# t-SNE perplexity  200   fixed-map held-label R2 for roll parameter  0.9924
# t-SNE perplexity  500   fixed-map held-label R2 for roll parameter  0.6393
# Isomap n_neighbors   5    kNN R2  0.9992
# Isomap n_neighbors  10    kNN R2  0.9996
# Isomap n_neighbors  30    kNN R2  0.9804
# PCA                        kNN R2  0.9995
# raw 3-D                    kNN R2  0.9997
# PCA      fold-fitted projection, held-out points R2 0.9763
# Isomap   fold-fitted projection, held-out points R2 0.9995

For this dataset and the tested settings, perplexities 5 through 200 all give held-label \(R^2>0.99\), while 2 and 500 score lower. These are six trials, not a guarantee for every value between them. A perplexity of 500 targets a broader affinity distribution than 5, but does not literally select one third of the sample. The score alone does not show whether the map fragmented or lost its topology.

PCA reaches 0.9995 in the fixed-map diagnostic. The prediction is made by a nonlinear k-nearest-neighbor regressor after the linear projection, so this does not imply the roll parameter is linear or monotone in either PCA coordinate. This projection retains enough information for that predictor on these points. The separate fold-fitted pipeline gives PCA 0.9763 and Isomap 0.9995 for observations excluded from fitting their projections. This changes the evaluation protocol, so it should not be mixed with the fixed-map scores.

Predicting the generating roll parameter is different from preserving distances along the sheet. In the noiseless generator, \(t=\sqrt{x_0^2+x_2^2}\): retaining those two coordinates preserves the parameter while discarding the sheet’s height coordinate. Recovering this one target therefore cannot establish that both intrinsic dimensions or their distances were preserved. The comparison in Exercise 1 therefore adds a direct intrinsic-distance measurement. A PCA baseline remains useful even when the intended geometric representation requires a nonlinear map.

UMAP: a neighborhood graph and a reusable map

UMAP connects nearby observations in a graph and assigns each edge a membership strength between zero and one, using locally adjusted distance scales. It combines the two directed strengths \(a\) and \(b\) as \(a+b-ab\); strengths 0.2 and 0.5 therefore become 0.6. A low-dimensional layout is fitted through attractive and repulsive updates motivated by a graph cross-entropy objective. The implementation samples edges and non-neighbor relationships to make optimization practical. These strengths encode graph relationships; they do not identify class membership or guarantee that the fitted layout recovers the original topology. The UMAP paper develops the construction.

n_neighbors controls the neighborhood scale used to build the graph. min_dist shapes the low-dimensional similarity curve and influences how tightly points pack; it is not an enforced minimum separation for every pair. metric defines input-space similarity, and n_components sets the output dimension, which can exceed two for downstream modeling. The parameter guide illustrates these choices. Neither a smaller min_dist nor a larger n_neighbors guarantees a better representation.

The following example fits UMAP on training digits without passing their labels, then transforms held-out observations into that reference map. A nearest-neighbor classifier tests one use of the representation. Training trustworthiness penalizes embedding neighbors that were distant in the original ranking; one is best. It checks local false neighbors, not global geometry or predictive generalization. The first UMAP call may also compile numerical routines, so it is not used as a speed benchmark. This code requires umap-learn in the Python environment. The recorded UMAP run uses version 0.5.12 with Numba 0.61.2.

import numpy as np
from umap import UMAP
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.manifold import trustworthiness
from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    X, y = load_digits(return_X_y=True)
    X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25,
                                             stratify=y, random_state=0)
    for neighbors, min_dist in ((15, 0.0), (15, 0.5), (50, 0.1)):
        reducer = UMAP(n_neighbors=neighbors, min_dist=min_dist, n_components=2,
                       n_jobs=1, random_state=0, transform_seed=0)
        Z_tr = reducer.fit_transform(X_tr)
        Z_te = reducer.transform(X_te)
        accuracy = KNeighborsClassifier(10).fit(Z_tr, y_tr).score(Z_te, y_te)
        print(f"UMAP neighbors={neighbors:2d} min_dist={min_dist:.1f}: "
              f"training trustworthiness {trustworthiness(X_tr, Z_tr, n_neighbors=10):.4f}, "
              f"held-out accuracy {accuracy:.4f}")
    pca = PCA(2).fit(X_tr)
    pca_accuracy = KNeighborsClassifier(10).fit(pca.transform(X_tr), y_tr).score(pca.transform(X_te), y_te)
    raw_accuracy = KNeighborsClassifier(10).fit(X_tr, y_tr).score(X_te, y_te)
    print(f"PCA-2 held-out accuracy {pca_accuracy:.4f}, raw 64-D {raw_accuracy:.4f}")
# UMAP neighbors=15 min_dist=0.0: training trustworthiness 0.9885, held-out accuracy 0.9778
# UMAP neighbors=15 min_dist=0.5: training trustworthiness 0.9850, held-out accuracy 0.9711
# UMAP neighbors=50 min_dist=0.1: training trustworthiness 0.9860, held-out accuracy 0.9511
# PCA-2 held-out accuracy 0.6356, raw 64-D 0.9756

The reduced and raw-space classifiers are evaluated on the same held-out rows, and the reducer never sees those labels. The highest recorded UMAP accuracy is 0.9778 (440 of 450 points), against 0.9756 (439 points) in the raw space. A one-observation difference on this split does not establish a general accuracy advantage; the UMAP representation also uses only two coordinates. If used to select settings, these rows become validation data and a further test set is needed. UMAP’s transform operation places observations relative to the fitted training representation; it is not a joint refit and does not ensure sensible extrapolation to a new population.

Practical rules

• Choose the input features, scaling, and distance metric to express relevant similarity. PCA preprocessing can reduce cost on high-dimensional inputs, but can also discard useful directions.
• Use a recorded initialization and seed, and compare several reasonable neighborhood settings. Agreement is evidence of stability under those choices, not proof of a population structure.
• Check specific plot-based claims with quantitative measures. Distances can be measured in the embedding, but their interpretation must be justified; apparent islands and occupied areas do not validate clusters or density.
• Scikit-learn’s TSNE has fit_transform but no new-point transform. This is a limitation of that implementation, not of every t-SNE extension. UMAP supports transforming new points relative to a fitted map.
• For prediction, fit preprocessing and any learned embedding on training rows only. Use validation data to select settings and separate test data to evaluate the selected pipeline.

Measuring the instability on real data

The handwritten digits provide 1,797 observations in 64 dimensions. For each layout, compute the 45 distances between the ten labeled class centroids and divide them by that layout’s largest centroid distance. This removes overall scale and is unaffected by rigid rotations or translations. Changing perplexity changes the neighborhood objective; changing the seed with a fixed objective probes optimization sensitivity.

import numpy as np
from sklearn.datasets import load_digits
from sklearn.manifold import TSNE
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    X, y = load_digits(return_X_y=True)

    def norm_class_dists(Z):
        C = np.array([Z[y == d].mean(0) for d in range(10)])
        D = np.linalg.norm(C[:, None] - C[None, :], axis=-1)
        return D[np.triu_indices(10, 1)] / D.max()

    print("perplexity changes the between-class distances")
    base = None
    for perp in (5, 30, 100):
        v = norm_class_dists(TSNE(2, perplexity=perp, init="pca",
                                  random_state=0).fit_transform(X))
        if base is None:
            base = v
            print(f"  perplexity {perp:3d}  (reference)")
        else:
            print(f"  perplexity {perp:3d}  max change {np.abs(v - base).max():.3f}"
                  f"   correlation with reference {np.corrcoef(v, base)[0, 1]:.3f}")

    print("random initialization changes them too")
    runs = [norm_class_dists(TSNE(2, perplexity=30, init="random",
                                  random_state=s).fit_transform(X)) for s in (0, 1, 2)]
    for i in range(3):
        for j in range(i + 1, 3):
            print(f"  seed {i} vs {j}  max change {np.abs(runs[i] - runs[j]).max():.3f}"
                  f"   correlation {np.corrcoef(runs[i], runs[j])[0, 1]:.3f}")
# perplexity changes the between-class distances
#   perplexity   5  (reference)
#   perplexity  30  max change 0.320   correlation with reference 0.891
#   perplexity 100  max change 0.293   correlation with reference 0.874
# random initialization changes them too
#   seed 0 vs 1  max change 0.655   correlation 0.622
#   seed 0 vs 2  max change 0.145   correlation 0.980
#   seed 1 vs 2  max change 0.628   correlation 0.625

Changing perplexity from 5 to 30 changes one normalized class distance by as much as 0.320. These are distances between centroids defined by the known labels; a class need not form a single island in the embedding. The comparison shows sensitivity to the selected neighborhood scale, not a change in the observations or labels.

With random initialization, seeds 0 and 1 differ by up to 0.655 in normalized class-centroid distance, and their distance profiles have correlation 0.622. The observed large-scale arrangement is sensitive to the starting layout in this comparison. The maximum difference summarizes the largest change among all class pairs, rather than establishing that every pair moved by that amount.

PCA initialization can improve consistency, but it does not universally remove seed dependence. Whether the PCA step is randomized depends on the implementation and solver policy. In the installed scikit-learn 1.9 implementation, PCA uses its automatic solver; these tall, low-feature-count examples take the covariance-eigendecomposition route. Other shapes or versions can use a randomized solver and produce seed-dependent starting coordinates. Tiny numerical differences can also propagate through optimization. Exercise 2 distinguishes the built-in PCA initialization from an explicitly fixed initial coordinate array.

A seed-sensitive conclusion needs a stability qualification; seed agreement alone does not make it true. Local groups, relative positions, and apparent density should each be checked against the particular claim being made. Repeating optimization addresses one source of uncertainty, while resampling observations and changing meaningful preprocessing choices address others.

Exercises

1. Neighborhoods, prediction, and intrinsic distance. Embed a noiseless Swiss roll with PCA, Isomap, LLE, and t-SNE. Compare neighbor overlap, prediction of its roll parameter, and rank agreement with distances along the original sheet.

High prediction accuracy does not imply faithful geometry. Evaluate the three columns as different properties of the same representation.

import numpy as np
from sklearn.manifold import TSNE, Isomap, LocallyLinearEmbedding
from sklearn.decomposition import PCA
from sklearn.datasets import make_swiss_roll
from sklearn.neighbors import NearestNeighbors, KNeighborsRegressor
from sklearn.model_selection import cross_val_score
from scipy.spatial.distance import pdist
from scipy.stats import spearmanr
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    X, t = make_swiss_roll(n_samples=1200, noise=0.0, random_state=0)
    arc = 0.5 * (t * np.sqrt(1 + t*t) + np.arcsinh(t))
    intrinsic_distances = pdist(np.column_stack((arc, X[:, 1])))
    nn_orig = NearestNeighbors(n_neighbors=11).fit(X).kneighbors(X)[1][:, 1:]

    def neighbour_overlap(Z):
        nn_new = NearestNeighbors(n_neighbors=11).fit(Z).kneighbors(Z)[1][:, 1:]
        return np.mean([len(set(a) & set(b)) / 10 for a, b in zip(nn_orig, nn_new)])

    for name, Z in (
            ("PCA", PCA(2).fit_transform(X)),
            ("Isomap", Isomap(n_neighbors=10, n_components=2).fit_transform(X)),
            ("LLE", LocallyLinearEmbedding(n_neighbors=12, n_components=2,
                                           random_state=0).fit_transform(X)),
            ("t-SNE", TSNE(2, perplexity=30, random_state=0, init="pca").fit_transform(X))):
        print(f"{name:8s} 10-NN overlap {neighbour_overlap(Z):.4f}"
              f"   held-label R2 {cross_val_score(KNeighborsRegressor(10), Z, t, cv=5).mean():.4f}"
              f"   intrinsic-distance rank correlation {spearmanr(pdist(Z), intrinsic_distances).statistic:.4f}")
# PCA      10-NN overlap 0.5430   held-label R2 0.8936   intrinsic-distance rank correlation 0.3794
# Isomap   10-NN overlap 0.8508   held-label R2 0.9994   intrinsic-distance rank correlation 0.9993
# LLE      10-NN overlap 0.6425   held-label R2 0.9974   intrinsic-distance rank correlation 0.7340
# t-SNE    10-NN overlap 0.8714   held-label R2 0.9988   intrinsic-distance rank correlation 0.6346

2. Reproducibility of a t-SNE layout. Compare random initialization, built-in PCA initialization, and an explicit fixed PCA start across seeds. Separate agreement between runs from agreement with the input geometry.

Use neighbor overlap and normalized centroid-distance variation. Rounded summary statistics cannot establish identical coordinate arrays.

import numpy as np
from itertools import combinations
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
from sklearn.datasets import make_blobs
from sklearn.neighbors import NearestNeighbors
from threadpoolctl import threadpool_limits

with threadpool_limits(limits=1):

    X, y = make_blobs(n_samples=900, centers=5, cluster_std=1.2, n_features=20,
                      random_state=0)

    def overlap(A, B, k=10):
        na = NearestNeighbors(n_neighbors=k + 1).fit(A).kneighbors(A)[1][:, 1:]
        nb = NearestNeighbors(n_neighbors=k + 1).fit(B).kneighbors(B)[1][:, 1:]
        return np.mean([len(set(a) & set(b)) / k for a, b in zip(na, nb)])

    def centroid_shape(Z):                          # pairwise centroid distances, scaled
        m = np.array([Z[y == c].mean(0) for c in range(5)])
        d = np.array([np.linalg.norm(m[i] - m[j]) for i, j in combinations(range(5), 2)])
        return d / d.max()

    for init in ("random", "pca"):
        runs = [TSNE(2, perplexity=30, init=init, random_state=s).fit_transform(X)
                for s in range(4)]
        loc = np.mean([overlap(a, b) for a, b in combinations(runs, 2)])
        shapes = np.array([centroid_shape(r) for r in runs])
        print(f"init={init:7s} mean 10-NN overlap between runs {loc:.4f}"
              f"   centroid-geometry spread {shapes.std(0).mean():.4f}")
    fixed = PCA(n_components=2, svd_solver="full").fit_transform(X).astype(np.float32)
    fixed = fixed / fixed[:, 0].std() * 1e-4
    controlled = [TSNE(2, perplexity=30, init=fixed.copy(), random_state=s).fit_transform(X)
                  for s in range(4)]
    print(f"explicit fixed PCA start: max coordinate difference "
          f"{max(np.abs(z-controlled[0]).max() for z in controlled[1:]):.6g}")
# init=random  mean 10-NN overlap between runs 0.4705   centroid-geometry spread 0.1338
# init=pca     mean 10-NN overlap between runs 1.0000   centroid-geometry spread 0.0000
# explicit fixed PCA start: max coordinate difference 0

3. Runtime and PCA preprocessing. Compare median fit times at several sample sizes and dimensions. Include the cost of PCA in the preprocessing route and keep the t-SNE iteration budget fixed.

Measure the ordering on your machine. Runtime alone does not establish embedding quality, and a PCA preprocessing step can change both the data geometry and the work required.

import time
import numpy as np
from sklearn.manifold import TSNE, Isomap
from sklearn.decomposition import PCA
from threadpoolctl import threadpool_limits

rng = np.random.default_rng(0)
def embed(M, reduce_first=False):
    if reduce_first:
        M = PCA(50, svd_solver="randomized", random_state=0).fit_transform(M)
    return TSNE(2, perplexity=30, max_iter=500, random_state=0,
                init="pca").fit_transform(M)
def median_time(fn):
    elapsed = []
    for _ in range(3):
        start = time.perf_counter()
        fn()
        elapsed.append(time.perf_counter()-start)
    return np.median(elapsed)
with threadpool_limits(limits=1):
    warmup = rng.normal(size=(150, 60))
    embed(warmup)
    Isomap(n_neighbors=10).fit_transform(warmup)
    for n in (500, 1000, 2000):
        X = rng.normal(size=(n, 300))
        raw = median_time(lambda: embed(X))
        reduced = median_time(lambda: embed(X, True))
        iso = median_time(lambda: Isomap(n_neighbors=10).fit_transform(X))
        print(f"n={n:4d} d=300: raw {raw:.2f}s, PCA+t-SNE {reduced:.2f}s, Isomap {iso:.2f}s")
    for d in (300, 3000):
        X = rng.normal(size=(1000, d))
        raw = median_time(lambda: embed(X))
        reduced = median_time(lambda: embed(X, True))
        print(f"n=1000 d={d:4d}: raw {raw:.2f}s, PCA+t-SNE {reduced:.2f}s, speedup {raw/reduced:.2f}x")
for n in (10_000, 100_000):
    print(f"one float64 distance matrix at n={n}: {8*n*n/2**30:.2f} GiB")
# n= 500 d=300: raw 1.22s, PCA+t-SNE 1.32s, Isomap 0.05s # varies by machine
# n=1000 d=300: raw 2.83s, PCA+t-SNE 3.71s, Isomap 0.21s # varies by machine
# n=2000 d=300: raw 7.27s, PCA+t-SNE 5.45s, Isomap 0.88s # varies by machine
# n=1000 d= 300: raw 4.32s, PCA+t-SNE 3.21s, speedup 1.34x # varies by machine
# n=1000 d=3000: raw 3.24s, PCA+t-SNE 3.31s, speedup 0.98x # varies by machine
# one float64 distance matrix at n=10000: 0.75 GiB
# one float64 distance matrix at n=100000: 74.51 GiB

References

• van der Maaten, L., & Hinton, G. (2008). Visualizing Data using t-SNE. Journal of Machine Learning Research.
• Tenenbaum, J. B., de Silva, V., & Langford, J. C. (2000). A Global Geometric Framework for Nonlinear Dimensionality Reduction. Science.