TLDR

• A projection is a claim about what to preserve. It is not a neutral window into high-dimensional reality, and visual separation alone does not establish clusters or a manifold.
• PCA, t-SNE and UMAP optimise different losses. Those losses reward linear variance, probabilistic neighbourhoods and a fuzzy neighbour graph respectively, so honest maps of one dataset can disagree.
• State the relation you need before choosing a method, then measure it in the original space and after projection across seeds, scales, metrics, resamples and credible nulls.
• The method stops helping when the reference metric, sampling process or feature representation is itself wrong. No projection test supplies missing causal, temporal, semantic or authoritative evidence.
• Use a map to generate a bounded hypothesis. Operationalise the hypothesis only after a predeclared invariant, independent evidence and a negative control survive the projection claim gate.

Part IEvery map chooses a loss

Imagine fitting a relief sculpture into a postcard. You may preserve the outline seen from above, the shortest walking paths between nearby points, or the relative height of major peaks. You cannot preserve every property after removing a dimension. A high-dimensional dataset presents the same constraint more severely. Two output coordinates cannot generally retain all pairwise distances, angles, neighbourhoods, densities, connected components and directions from twenty, two hundred or two thousand coordinates.

The useful question is therefore not whether a map is distorted. It is which distortion it makes cheap. A method defines a loss function, explicit or implicit. Optimisation moves points until violations that the loss prices heavily become small. Relations priced weakly may deform drastically. The picture is the visible residue of that economic choice.

A minimal worked example: one point, two errors

Take a focal record A. In the declared high-dimensional metric, its two nearest neighbours are B and C. On the map, the two nearest visible points are B and D. B is preserved. D is a false neighbour: the map makes it look locally related when it was not. C is a missing neighbour: an original local relation has been torn apart.

These two errors answer different questions. Trustworthiness asks whether neighbours displayed by the map were genuine in the reference space. Continuity asks whether neighbours from the reference space remain visible. A map can score well on one and poorly on the other. Venna and Kaski developed this distinction explicitly for nonlinear projection quality [7]; scikit-learn exposes the trustworthiness measure and its rank penalty [8].

Thought experiment 1: three cartographers, one cloud

Give three cartographers the same point cloud and the same two square centimetres of paper. Tell the first to preserve the directions with greatest variation. Tell the second to keep each point beside the points it most strongly resembles. Tell the third to preserve a weighted local graph built with a separate scale around every point. Hold the data fixed. Only the contractual promise changes.

The first draws a broad arc, the second makes compact islands, and the third draws a chain with a narrow break. Asking which picture is the true shape confuses an optimisation target with an observation. The discriminating question is which pictured relation survives when checked against the original measurements. That is the causal role of the objective: change the objective while holding the data constant, and the map changes.

This paper builds on the earlier geometric distinction that coordinates are reports relative to a basis and metric, rather than self-interpreting facts. The adjacent article Vectors as Directions, not Lists [16] establishes that foundation. Here the new burden is controlled loss: even with a valid reference geometry, reducing dimension chooses which relations to damage.

The reference space is already a model

Projection quality is always relative to something. If Euclidean distance in the input space treats a one-standard-deviation change in every feature as comparable, then trustworthiness rewards the map for preserving that choice. Change from raw units to standardised units and the nearest-neighbour graph can change before PCA, t-SNE or UMAP is fitted. Change from Euclidean to cosine distance and radial information stops contributing. Remove a duplicated feature and a direction that previously counted twice loses influence.

Consider a case matrix with twenty timing fields and one verified harm indicator. Standardisation gives every column unit variance, but the timing family still has twenty votes. A projection can faithfully preserve the resulting neighbourhoods while burying the harm relation. Conversely, giving the harm indicator a large manual weight may make it dominate every local map. Neither choice is neutral. Each is a measurement hypothesis about which differences should make two cases dissimilar.

Fidelity to an undeclared metric is not scientific validity. The reference contract should name feature provenance, units, transformations, missing-value treatment, duplicated signals, temporal validity and the reason the metric matches the claim. When several contracts are plausible, run them as competing hypotheses. A structure that exists only under one unowned scaling rule should be reported as metric-sensitive, not as a discovered population.

flowchart TD
  A["The map that found two businesses"] --> B["Part I Every map chooses a loss"]
  B --> C["Part II Three preservation<br/>contracts"]
  C --> D["Part III When a manifold claim<br/>outruns the map"]
  D -. "re-test" .-> A
  classDef first fill:#0a1424,stroke:#5B5CF6,stroke-width:2px,color:#e6f3ff
  classDef step fill:#0a1424,stroke:#3f7da0,stroke-width:1px,color:#e6f3ff
  class A first
  class B,C,D step

Part IIThree preservation contracts

PCA, t-SNE and UMAP are often placed side by side as visualisation options. That catalogue view hides the important difference. Each method creates a different evidential object. The method name does not tell us whether the resulting object supports a cluster, ordering, distance or density claim. Its loss and assumptions do.

PCA preserves the best linear reconstruction

Start with centred observations in a matrix X, with n rows and d features. PCA chooses orthonormal directions W whose two columns span the plane that minimises squared reconstruction error. Equivalently, those directions maximise captured variance. Pearson formulated the closest-fit geometry at the beginning of the twentieth century [1], and modern implementations compute it through eigendecomposition or singular value decomposition [2].

X is the centred n by d data matrix. W is a d by 2 orthonormal basis. XWWᵀ is the reconstruction obtained after projecting into two dimensions and lifting back. Adding a component cannot increase the optimal reconstruction error; changing feature scale can change W completely.

PCA therefore supplies interpretable linear axes and an explicit residual. It is an excellent baseline when a few global directions carry the relevant variance. It can fail visually when the data bend around a nonlinear curve: points far apart along the curve may overlap in the linear shadow. But the failure is informative. It tells us that no two-dimensional linear plane retains the measured variance, rather than pretending that the curve has been unfolded without cost.

PCA supports a statement about linear variance or reconstruction, not automatically about semantic importance. A low-variance direction can still govern an outcome. Standardisation, weighting and feature engineering determine what variance means. The first principal component is not “the most important factor” unless importance has independently been defined by variance under that measurement contract.

t-SNE prices pairwise probability errors

t-SNE begins locally. It converts high-dimensional distances into pairwise similarities and normalises them into probabilities p_ij. In the map it defines probabilities q_ij with a heavy-tailed Student distribution. Optimisation minimises Kullback-Leibler divergence from P to Q [3]. Official documentation also warns that the loss is non-convex, so different initialisations can yield different results [4].

p_ij is the high-dimensional similarity probability; q_ij is its two-dimensional counterpart. A large p with a small q contributes strongly, so separating a pair deemed close is expensive. Perplexity changes the effective neighbourhood scale from which p is built.

This explains both t-SNE's power and its visual temptations. Local groups can become remarkably legible because the loss spends much of its budget on near relationships and uses repulsion to avoid crowding. Yet the white space between groups is not a calibrated distance ruler. Island size is not a reliable density estimate. Early exaggeration, perplexity, learning rate, initialisation and random seed can alter the visible arrangement. Practical guides have repeatedly demonstrated that cluster sizes and inter-cluster distances can be misleading when read literally [14].

There is no contradiction between a good local map and a poor global geometry. In the sparse-bridge laboratory, t-SNE at perplexity 30 attains trustworthiness 0.998 at k = 10, while the Spearman correlation between sampled original and map distances is 0.911. At perplexity 5, trustworthiness remains 0.997 but distance correlation falls to 0.603. The local question is still answered well; the global question changed.

UMAP turns local scales into a fuzzy graph

UMAP begins with a k-nearest-neighbour graph, but it does not treat every neighbourhood with one global distance scale. Around each point it estimates a local connectivity offset and scale, converts neighbour relations into directed membership strengths, then combines opposite directions through a fuzzy union. The original formulation connects this construction to Riemannian geometry and fuzzy simplicial sets [5].

μ_ij is the symmetric high-dimensional edge membership. ν_ij is the low-dimensional membership determined by mapped coordinates y and curve parameters a and b. Attraction reproduces strong edges; sampled repulsion separates non-edges.

The parameter n_neighbors changes the scale at which the graph is built. Smaller values give more authority to very local relations; larger values incorporate a wider context. min_dist changes how tightly the optimiser may pack points in the map [6]. These are not cosmetic controls. They alter the preservation contract and can change whether a bridge remains visible, whether small variations become islands and how compact a group looks.

UMAP's manifold language must be read with its assumptions. The software documentation names a uniformly sampled Riemannian manifold, a locally constant metric and local connectedness as foundations. Real data may violate all three. Local normalisation is useful when sampling density should not dominate similarity, but it also means standard UMAP is not a density-preserving instrument. densMAP was introduced precisely to add density preservation to the objective [11].

UMAP preserves evidence about a constructed fuzzy neighbourhood graph, not the full topology or density of an unknown generating process. The graph depends on features, metric, scaling, sample density and k. If the graph tears a sparse region, the optimiser cannot recover a connection it was never given.

Scale parameters are causal interventions

Perplexity and neighbourhood size are often tuned until a plot looks legible. That practice reverses the evidential order. Each value asks what should count as local. A small setting allows fine-grained relations to dominate; a larger setting forces each point to negotiate with a wider population. If the claimed group exists only at one narrow scale, that scale dependence is part of the result.

A disciplined atlas therefore shows a range rather than a winner. Hold preprocessing, metric, sample and initialisation fixed. Sweep the neighbourhood scale. Record quality measures and the fate of the named relation, such as whether a bridge remains connected or whether membership is stable. Then vary initialisation while holding scale fixed. This factorial design distinguishes a scale-sensitive phenomenon from an optimiser-sensitive drawing.

There is no universal parameter at which “the manifold” appears. Real datasets can contain structure at several scales: households inside regions, cell states inside lineages, or transaction motifs inside customer journeys. A useful map may expose one level while erasing another. The claim must state its scale and denominator rather than borrowing certainty from a default value.

The baseline is part of the explanation

A nonlinear map looks impressive partly because it is viewed without a matched linear baseline. PCA reveals whether the dramatic separation was already present in a simple global projection. A map of shuffled or simulated null data reveals whether the algorithm and human visual system can manufacture equally persuasive groupings without the proposed phenomenon. A no-reduction model tests whether compressing to two coordinates helps the actual task at all.

Research comparing modern methods finds a persistent local-global trade-off rather than a universal ranking [9]. A separate controlled analysis showed that an apparent global advantage of UMAP over t-SNE could be attributed to different default initialisations; when initialisation was matched, both changed materially [10]. The practical implication is severe: method comparison must vary one causal feature at a time. Comparing package defaults compares objectives, initialisation, optimisation schedules and parameter conventions together.

Part IIIWhen a manifold claim outruns the map

A manifold claim says more than “these dots can be arranged attractively in two dimensions”. It proposes that observations concentrate near a lower-dimensional structure embedded in a higher-dimensional space, with some meaningful local geometry. Stronger claims may assert connected components, branches, loops, trajectories or an intrinsic dimension. A finite scatter plot establishes none of these by appearance alone.

Any finite set can be assigned two coordinates. The difficult question is whether the map preserves relations that correspond to the generating process rather than relations induced by sampling, measurement and optimisation. A neighbour graph is built from the observed sample, not from the unobserved manifold. Sparse coverage, noise and mixed metrics can change its topology before projection starts.

Intrinsic dimension is a separate empirical question

A successful two-dimensional visualisation does not show that the data have intrinsic dimension two. The optimiser is required to return two coordinates regardless of whether the relevant variation needs two, twelve or two hundred degrees of freedom. It can fold, overlap or tear relations to comply. Conversely, poor visual separation does not prove high intrinsic dimension; noise, a poor metric or uneven sampling can obscure a simple generator.

Intrinsic dimension estimators ask how neighbourhood volume, neighbour distances, reconstruction or local tangent rank scale as more observations are considered. They bring their own assumptions and finite-sample biases, but they address the dimensionality question directly. A responsible manifold claim triangulates such diagnostics with domain mechanism and projection stability. It does not infer dimension from the number of axes requested by an API.

Topology also needs a named invariant. “The map preserves topology” could mean only that close points tend to stay close, or it could mean connected components, loops and branch junctions survive. Those are different tests. If a decision depends on a loop, measure that loop across neighbourhood scales and resamples. If it depends on a branch, test the branch point against new observations. The word manifold should narrow the experiment, not decorate the caption.

Thought experiment 2: thin only the bridge

Imagine a continuous behavioural process indexed by t. At low t, cases are simple and slow. At high t, they are complex and urgent. Intermediate cases exist, so the generating curve is connected. First sample evenly along the curve. Then hold the curve, features, noise and metric fixed while discarding most observations from a narrow middle interval.

No ontology changed. Only sampling density changed. Yet a fixed-k neighbour graph may stop crossing the sparse interval. Local scales expand differently at the two dense ends. Attraction operates within each end, repulsion enlarges the gap, and the map becomes two islands. A viewer can now invent two populations from a process that remained continuous.

The 2025 map-continuity study makes this concern formal. It argues that neighbour embeddings lack a sample-independent continuous map in the ordinary formulation and documents discontinuities that can exaggerate separation or create spurious local structures [12]. Earlier broad evaluations likewise warn that low-quality projections can display structure absent from the high-dimensional data [13].

Three common promotions therefore need separate evidence:

• Cluster claim: a gap or island is treated as a population boundary. Required evidence includes stability under resampling, metrics and seeds, separation in a declared reference space, and an external variable or mechanism not used to create the map.
• Trajectory claim: a curved ordering is treated as time, progression or development. Required evidence includes temporal or intervention data, directionality and branch stability. Colour gradients alone can be circular when the colour variable helped form the embedding.
• Manifold claim: local adjacency is treated as intrinsic geometry or topology. Required evidence includes intrinsic-dimension diagnostics, sampling adequacy and tests of connectedness, loops or geodesic relations at multiple scales.

Density, area and distance are especially easy to overread

Point density in a map is produced jointly by original sampling density, local distance normalisation, attraction, repulsion and the minimum-distance rule. A large island may contain more observations, a more variable population, or merely an optimiser that spread it. A compact island may be homogeneous or simply permitted to pack tightly. Standard t-SNE and UMAP were not designed to make island area proportional to high-dimensional variance or probability mass.

Likewise, the distance between island centroids has no universal interpretation. It can change under rotation harmlessly, under initialisation materially, and under local objectives without representing a comparable change in original distance. An analyst may still use such maps for exploration, but the caption should say what is and is not metrically meaningful.

A negative control and a serious counterexample

A sceptical account must explain both false structure and real structure. The lab therefore includes an isotropic Gaussian null with 360 observations in 24 dimensions. There are no planted groups. K-means in the original standardised space obtains silhouette 0.036. After t-SNE, a three-way partition obtains silhouette 0.366. The map is more partitionable, but partitions across three random initialisations have mean adjusted Rand agreement only 0.096. Visual clarity increased while evidential stability collapsed.

The counterexample uses three genuinely separated Gaussian populations. Their ground-truth silhouette is 0.781 in the original standardised space. PCA, t-SNE and the UMAP reference all preserve clear separation, with map silhouettes 0.937, 0.887 and 0.879. Here the visual claim is not dismissed. It is supported because separation existed before projection and survives objectives with different biases.

Part IVStress-test the claim

A good projection review begins with a sentence that can fail: “Cases adjacent in this map should have interchangeable resolution paths”, or “The bridge represents a continuous progression in elapsed risk”, or “The three islands are stable populations with different outcomes”. The map is then an instrument for that sentence. Without the sentence, reviewers can admire shapes indefinitely without knowing what would count as error.

Production-shaped worked example: source system masquerading as segment

The opening case is generated, not a deployment. It contains 540 synthetic cases with continuous urgency, complexity and age factors. A binary source-system field changes the extraction convention for six of 36 observed features. Later outcome is sampled from urgency, complexity and age; source does not directly determine it. The UMAP reference map is held fixed below while only the annotation changes.

The right intervention is to harmonise or remove the source-specific fields, rerun the map, and test whether outcome or process differences remain in held-out data. A second intervention reweights features according to a domain measurement contract. A third resamples each source to equalise coverage. If the islands disappear while the business factors remain continuous, the segment claim weakens. If independent outcomes, interventions and operational behaviours remain separated, the claim strengthens.

The changed decision is to treat batch-aligned geometry as evidence about the data pipeline before treating it as evidence about the world. That ordering prevents a visualisation tool from laundering a source-system boundary into a customer ontology.

Use four families of evidence

No scalar quality score can certify a projection for every claim. Use a small portfolio whose members fail differently:

1. Neighbour evidence. Measure trustworthiness, continuity and direct k-neighbour overlap at several k values. Report the denominator and reference metric. High values support local fidelity only at those scales.
2. Global evidence. Compare rank correlation of pairwise distances, reconstruction error for linear methods, landmark distances and connected components. Sampling pairs is acceptable when the sampling rule is recorded.
3. Stability evidence. Vary seed, initialisation, neighbourhood size, perplexity, minimum distance, feature scaling, metric and resample. Align maps only for coordinate comparison; compare neighbourhoods and claims directly when rotation is irrelevant.
4. External evidence. Test outcomes, labels, time, interventions or mechanisms not used to fit the projection. Hold these out when possible. A label-coloured map is descriptive, not predictive validation.

Trustworthiness illustrates why the denominator matters. For n observations and neighbourhood size k, the score penalises displayed neighbours whose rank in the original space lies beyond k. The penalty grows with rank error. It approaches one when false neighbours are few or only just outside the boundary. It does not test whether original neighbours were lost, whether density was preserved, or whether the reference distance was meaningful.

Do not operationalise the two coordinates by accident

A common shortcut fits a visually pleasing map on all available observations, then feeds its two coordinates into clustering, rules or a classifier. This changes the role of the projection. It is no longer an explanatory display; it is a learned feature transformation whose randomness, fitting population and out-of-sample behaviour enter the production model. Validation must be repeated inside each training fold, and new points must be transformed without using future observations.

t-SNE in its classical form does not provide a simple fixed parametric mapping for unseen points. UMAP can transform new observations relative to a fitted model, but that does not guarantee that a future population shares the training graph or local scales. Re-fitting can move every old coordinate. An identifier such as “island 2” therefore lacks stable meaning unless membership is defined independently of map position and versioned across releases.

The safer pattern is to test the substantive relation in the original or a task-validated representation. Use the map to inspect errors, coverage and candidate structure. If a cluster becomes operational, define it through reproducible features or a separately validated model, assign versioned membership, and monitor transitions. Keep the projection as evidence about the model, not as an untyped source of authority.

The projection claim card

The following card is practitioner synthesis. It is not a published standard. Its purpose is to stop an exploratory picture from acquiring more authority than its evidence. Fill it before the final map is selected, because post-hoc choice encourages the clearest picture to define the claim.

Executable artefact: PCA, t-SNE and UMAP projection claim lab

The companion Python programme creates four synthetic environments: a sparse continuous bridge, three separated populations, an isotropic Gaussian null and the service-case estate. It fits PCA and t-SNE with scikit-learn. When umap-learn is installed it uses the official implementation; otherwise it runs a compact reference path that constructs UMAP's fuzzy graph and optimises an attraction-repulsion objective for small datasets.

The assumptions are explicit: Euclidean distance after feature standardisation, synthetic public data, fixed seeds, no claim that two dimensions preserve every relation, and a pedagogical fallback that omits approximate-neighbour search and production optimisations. The positive assertion requires planted populations to remain separated. The negative assertion requires a t-SNE partition of the Gaussian null to look stronger than its original partition while changing substantially across seeds.

#!/usr/bin/env python3
"""Projection Claim Lab: stress-test PCA, t-SNE and UMAP claims on synthetic data.

The script uses umap-learn when available. Without it, it falls back to a compact
reference implementation that follows UMAP's fuzzy-neighbour graph and
attraction/repulsion objective. The fallback is for transparent small-data
experiments, not a speed-compatible replacement for umap-learn.

Requirements: numpy, scipy and scikit-learn. umap-learn is optional.
"""
from __future__ import annotations

import argparse
import json
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, Tuple

import numpy as np
from numpy.typing import NDArray
from scipy.optimize import curve_fit
from scipy.spatial.distance import pdist
from scipy.stats import spearmanr
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, SpectralEmbedding, trustworthiness
from sklearn.metrics import adjusted_rand_score, pairwise_distances, silhouette_score
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

Array = NDArray[np.float64]

try:
    import umap  # type: ignore
except ImportError:
    umap = None


@dataclass(frozen=True)
class EmbeddingScore:
    trust_10: float
    continuity_10: float
    local_overlap_10: float
    distance_rho: float


def zscore(x: Array) -> Array:
    return StandardScaler().fit_transform(x).astype(np.float64)


def make_sparse_bridge(n: int = 520, d: int = 24, seed: int = 7) -> Tuple[Array, Array]:
    """Continuous one-dimensional process sampled densely at its ends and sparsely in its middle."""
    rng = np.random.default_rng(seed)
    n_bridge = max(24, n // 13)
    n_end = (n - n_bridge) // 2
    t = np.concatenate([
        rng.uniform(-3.2, -0.65, n_end),
        rng.uniform(-0.65, 0.65, n_bridge),
        rng.uniform(0.65, 3.2, n - n_end - n_bridge),
    ])
    t.sort()
    base = np.column_stack([
        t,
        np.sin(1.35 * t),
        np.cos(1.35 * t),
        0.20 * t**2,
        np.sin(0.55 * t + 0.4),
    ])
    mix = rng.normal(size=(base.shape[1], d))
    x = base @ mix + rng.normal(scale=0.12, size=(n, d))
    return zscore(x), t


def make_three_populations(n: int = 480, d: int = 20, seed: int = 11) -> Tuple[Array, Array]:
    rng = np.random.default_rng(seed)
    centres = np.vstack([
        np.full(d, -3.5),
        np.full(d, 3.5),
        np.where(np.arange(d) % 2 == 0, 3.5, -3.5),
    ])
    x, y = make_blobs(n_samples=n, centers=centres, cluster_std=0.78, random_state=seed)
    # Add a shared, low-amplitude nuisance direction so the control is not axis-trivial.
    x += rng.normal(scale=0.08, size=x.shape)
    return zscore(x.astype(np.float64)), y.astype(int)


def make_gaussian_null(n: int = 480, d: int = 24, seed: int = 19) -> Array:
    rng = np.random.default_rng(seed)
    return zscore(rng.normal(size=(n, d)).astype(np.float64))


def make_case_estate(n: int = 720, d: int = 36, seed: int = 23) -> Tuple[Array, Dict[str, Array]]:
    """Synthetic service-case estate with continuous business factors and an acquisition-system batch effect."""
    rng = np.random.default_rng(seed)
    urgency = rng.beta(2.2, 2.0, n)
    complexity = np.clip(0.20 + 0.68 * rng.beta(2.0, 2.7, n) + 0.12 * urgency + rng.normal(0, 0.05, n), 0, 1)
    age = np.clip(rng.gamma(shape=2.0, scale=0.16, size=n), 0, 1)
    source = (rng.random(n) < (0.46 + 0.13 * (age > 0.55))).astype(int)
    phase = np.clip(0.55 * urgency + 0.35 * complexity + rng.normal(0, 0.08, n), 0, 1)

    latent = np.column_stack([
        urgency,
        complexity,
        age,
        phase,
        urgency * complexity,
        np.sin(np.pi * phase),
        np.cos(np.pi * phase),
    ])
    mix = rng.normal(size=(latent.shape[1], d))
    x = latent @ mix
    # A source-system signature affects six measured fields without changing the underlying case continuum.
    x[:, :6] += source[:, None] * np.array([1.8, 1.3, 1.0, -1.4, 0.9, -1.1])
    x += rng.normal(scale=0.26, size=x.shape)
    logits = -2.5 + 2.9 * urgency + 1.5 * complexity + 0.9 * age
    outcome_p = 1 / (1 + np.exp(-logits))
    outcome = (rng.random(n) < outcome_p).astype(int)
    return zscore(x.astype(np.float64)), {
        "urgency": urgency.astype(np.float64),
        "complexity": complexity.astype(np.float64),
        "source": source.astype(np.float64),
        "outcome": outcome.astype(np.float64),
    }


def _find_ab(spread: float = 1.0, min_dist: float = 0.1) -> Tuple[float, float]:
    def curve(x: Array, a: float, b: float) -> Array:
        return 1.0 / (1.0 + a * x ** (2.0 * b))

    xv = np.linspace(0, spread * 3.0, 300)
    yv = np.where(xv < min_dist, 1.0, np.exp(-(xv - min_dist) / spread))
    (a, b), _ = curve_fit(curve, xv, yv, p0=(1.5, 0.9), maxfev=20_000)
    return float(a), float(b)


def _smooth_knn_weights(x: Array, n_neighbors: int) -> Array:
    nn = NearestNeighbors(n_neighbors=n_neighbors + 1).fit(x)
    distances, indices = nn.kneighbors(x)
    distances, indices = distances[:, 1:], indices[:, 1:]
    n = x.shape[0]
    target = math.log2(n_neighbors)
    directed = np.zeros((n, n), dtype=np.float64)

    for i in range(n):
        row = distances[i]
        positive = row[row > 0]
        rho = positive[0] if positive.size else 0.0
        lo, hi = 0.0, np.inf
        sigma = 1.0
        for _ in range(64):
            score = np.exp(-np.maximum(0.0, row - rho) / max(sigma, 1e-12)).sum()
            if abs(score - target) < 1e-5:
                break
            if score > target:
                hi = sigma
                sigma = (lo + hi) / 2.0
            else:
                lo = sigma
                sigma = sigma * 2.0 if np.isinf(hi) else (lo + hi) / 2.0
        values = np.exp(-np.maximum(0.0, row - rho) / max(sigma, 1e-12))
        directed[i, indices[i]] = values

    # Probabilistic union of the two directed local fuzzy sets.
    return directed + directed.T - directed * directed.T


def small_umap_reference(
    x: Array,
    n_neighbors: int = 15,
    min_dist: float = 0.1,
    seed: int = 42,
    epochs: int = 280,
    negative_sample_rate: int = 5,
) -> Array:
    """Small-data UMAP reference: fuzzy graph plus attraction/repulsion SGD.

    It intentionally omits approximate-neighbour search, epoch scheduling details,
    transform support and compiled kernels used by umap-learn.
    """
    rng = np.random.default_rng(seed)
    w = _smooth_knn_weights(x, n_neighbors)
    rows, cols = np.where(np.triu(w, 1) > 1e-5)
    weights = w[rows, cols]
    weights = weights / weights.sum()
    a, b = _find_ab(min_dist=min_dist)

    # A graph-informed initialisation. PCA is the safe fallback if the graph is awkward.
    try:
        y = SpectralEmbedding(
            n_components=2,
            affinity="precomputed",
            random_state=seed,
            eigen_solver="arpack",
        ).fit_transform(w)
    except Exception:
        y = PCA(n_components=2, random_state=seed).fit_transform(x)
    y = y / max(np.std(y), 1e-9) * 0.0001
    y += rng.normal(scale=1e-5, size=y.shape)

    batch = min(10_000, max(3_000, len(rows) * 2))
    for epoch in range(epochs):
        alpha = 1.0 * (1.0 - epoch / max(epochs, 1)) + 0.02
        pick = rng.choice(len(rows), size=batch, replace=True, p=weights)
        i, j = rows[pick], cols[pick]
        diff = y[i] - y[j]
        d2 = np.einsum("ij,ij->i", diff, diff) + 1e-7
        attraction = -2.0 * a * b * np.power(d2, b - 1.0) / (1.0 + a * np.power(d2, b))
        grad = np.clip(attraction[:, None] * diff, -4.0, 4.0)
        np.add.at(y, i, alpha * grad / batch * 180.0)
        np.add.at(y, j, -alpha * grad / batch * 180.0)

        neg_n = batch * negative_sample_rate
        ni = rng.integers(0, len(y), neg_n)
        nj = rng.integers(0, len(y), neg_n)
        keep = ni != nj
        ni, nj = ni[keep], nj[keep]
        diff = y[ni] - y[nj]
        d2 = np.einsum("ij,ij->i", diff, diff) + 1e-3
        repulsion = 2.0 * b / (d2 * (1.0 + a * np.power(d2, b)))
        grad = np.clip(repulsion[:, None] * diff, -4.0, 4.0)
        np.add.at(y, ni, alpha * grad / max(len(ni), 1) * 180.0)
        np.add.at(y, nj, -alpha * grad / max(len(ni), 1) * 180.0)
        y -= y.mean(axis=0, keepdims=True)

    return y.astype(np.float64)


def fit_umap(x: Array, seed: int, n_neighbors: int = 15, min_dist: float = 0.1) -> Tuple[Array, str]:
    if umap is not None:
        model = umap.UMAP(
            n_components=2,
            n_neighbors=n_neighbors,
            min_dist=min_dist,
            metric="euclidean",
            random_state=seed,
        )
        return model.fit_transform(x).astype(np.float64), "umap-learn"
    return small_umap_reference(x, n_neighbors, min_dist, seed), "small reference"


def fit_all(x: Array, seed: int = 42, perplexity: float = 30.0, n_neighbors: int = 15) -> Tuple[Dict[str, Array], str]:
    pca = PCA(n_components=2, random_state=seed).fit_transform(x)
    tsne = TSNE(
        n_components=2,
        perplexity=min(perplexity, (len(x) - 1) / 3),
        init="pca",
        learning_rate="auto",
        max_iter=1_000,
        random_state=seed,
    ).fit_transform(x)
    umap_y, implementation = fit_umap(x, seed=seed, n_neighbors=n_neighbors)
    return {"PCA": pca, "t-SNE": tsne, "UMAP": umap_y}, implementation


def _rank_matrix(x: Array) -> NDArray[np.int64]:
    distances = pairwise_distances(x, metric="euclidean")
    order = np.argsort(distances, axis=1)
    ranks = np.empty_like(order)
    rows = np.arange(len(x))[:, None]
    ranks[rows, order] = np.arange(len(x))[None, :]
    return ranks


def continuity(x: Array, y: Array, k: int = 10) -> float:
    n = len(x)
    if not 1 <= k < n / 2:
        raise ValueError("k must be at least 1 and smaller than n / 2")
    rx = _rank_matrix(x)
    ry = _rank_matrix(y)
    original = rx <= k
    embedded = ry <= k
    missing = original & ~embedded
    penalty = np.maximum(0, ry - k)[missing].sum()
    normaliser = n * k * (2 * n - 3 * k - 1)
    return float(1.0 - 2.0 * penalty / normaliser)


def local_overlap(x: Array, y: Array, k: int = 10) -> float:
    rx = _rank_matrix(x)
    ry = _rank_matrix(y)
    overlaps = []
    for i in range(len(x)):
        nx = set(np.where((rx[i] >= 1) & (rx[i] <= k))[0])
        ny = set(np.where((ry[i] >= 1) & (ry[i] <= k))[0])
        overlaps.append(len(nx & ny) / k)
    return float(np.mean(overlaps))


def distance_rho(x: Array, y: Array, max_pairs: int = 60_000, seed: int = 5) -> float:
    dx, dy = pdist(x), pdist(y)
    if len(dx) > max_pairs:
        rng = np.random.default_rng(seed)
        ix = rng.choice(len(dx), max_pairs, replace=False)
        dx, dy = dx[ix], dy[ix]
    return float(spearmanr(dx, dy).statistic)


def score_embedding(x: Array, y: Array, k: int = 10) -> EmbeddingScore:
    return EmbeddingScore(
        trust_10=float(trustworthiness(x, y, n_neighbors=k)),
        continuity_10=continuity(x, y, k=k),
        local_overlap_10=local_overlap(x, y, k=k),
        distance_rho=distance_rho(x, y),
    )


def neighbourhood_jaccard(a: Array, b: Array, k: int = 15) -> float:
    ra, rb = _rank_matrix(a), _rank_matrix(b)
    values = []
    for i in range(len(a)):
        na = set(np.where((ra[i] >= 1) & (ra[i] <= k))[0])
        nb = set(np.where((rb[i] >= 1) & (rb[i] <= k))[0])
        values.append(len(na & nb) / len(na | nb))
    return float(np.mean(values))


def kmeans_map_labels(y: Array, seed: int = 0, k: int = 3) -> NDArray[np.int64]:
    return KMeans(n_clusters=k, n_init=20, random_state=seed).fit_predict(y)


def run_lab(output: Path, quick: bool = False) -> Dict[str, object]:
    output.mkdir(parents=True, exist_ok=True)
    bridge_x, bridge_t = make_sparse_bridge(n=390 if quick else 520)
    blobs_x, blobs_label = make_three_populations(n=360 if quick else 480)
    null_x = make_gaussian_null(n=360 if quick else 480)
    cases_x, case_meta = make_case_estate(n=540 if quick else 720)

    bridge_maps, implementation = fit_all(bridge_x, seed=42)
    blob_maps, _ = fit_all(blobs_x, seed=42)
    null_maps, _ = fit_all(null_x, seed=42)
    case_maps, _ = fit_all(cases_x, seed=42)

    scores: Dict[str, Dict[str, Dict[str, float]]] = {}
    for name, x, maps in [
        ("sparse_bridge", bridge_x, bridge_maps),
        ("three_populations", blobs_x, blob_maps),
        ("gaussian_null", null_x, null_maps),
        ("case_estate", cases_x, case_maps),
    ]:
        scores[name] = {m: score_embedding(x, y).__dict__ for m, y in maps.items()}

    # Seed intervention on the null. Initialisation and optimiser change; the input stays fixed.
    null_tsne_seeds = [
        TSNE(n_components=2, perplexity=30, init="random", learning_rate="auto", max_iter=1_000, random_state=s).fit_transform(null_x)
        for s in (1, 2, 3)
    ]
    null_umap_seeds = [fit_umap(null_x, seed=s)[0] for s in (1, 2, 3)]
    tsne_seed_jaccard = np.mean([
        neighbourhood_jaccard(null_tsne_seeds[i], null_tsne_seeds[j])
        for i in range(3) for j in range(i + 1, 3)
    ])
    umap_seed_jaccard = np.mean([
        neighbourhood_jaccard(null_umap_seeds[i], null_umap_seeds[j])
        for i in range(3) for j in range(i + 1, 3)
    ])
    tsne_cluster_ari = np.mean([
        adjusted_rand_score(kmeans_map_labels(null_tsne_seeds[i]), kmeans_map_labels(null_tsne_seeds[j]))
        for i in range(3) for j in range(i + 1, 3)
    ])
    umap_cluster_ari = np.mean([
        adjusted_rand_score(kmeans_map_labels(null_umap_seeds[i]), kmeans_map_labels(null_umap_seeds[j]))
        for i in range(3) for j in range(i + 1, 3)
    ])

    null_original_labels = KMeans(n_clusters=3, n_init=20, random_state=0).fit_predict(null_x)
    null_original_sil = silhouette_score(null_x, null_original_labels)
    null_tsne_sil = silhouette_score(null_tsne_seeds[0], kmeans_map_labels(null_tsne_seeds[0]))
    null_umap_sil = silhouette_score(null_umap_seeds[0], kmeans_map_labels(null_umap_seeds[0]))

    blob_original_sil = silhouette_score(blobs_x, blobs_label)
    blob_map_sil = {m: float(silhouette_score(y, blobs_label)) for m, y in blob_maps.items()}

    # PCA loadings are stable up to sign here; nonlinear map axes have no comparable feature meaning.
    pca_model = PCA(n_components=2).fit(cases_x)
    explained = pca_model.explained_variance_ratio_.tolist()

    results: Dict[str, object] = {
        "implementation": implementation,
        "scores": scores,
        "seed_intervention": {
            "tsne_neighbour_jaccard": float(tsne_seed_jaccard),
            "umap_neighbour_jaccard": float(umap_seed_jaccard),
            "tsne_cluster_ari": float(tsne_cluster_ari),
            "umap_cluster_ari": float(umap_cluster_ari),
        },
        "null_control": {
            "original_kmeans_silhouette": float(null_original_sil),
            "tsne_map_silhouette": float(null_tsne_sil),
            "umap_map_silhouette": float(null_umap_sil),
        },
        "positive_control": {
            "original_label_silhouette": float(blob_original_sil),
            "map_label_silhouette": blob_map_sil,
        },
        "case_pca_explained_variance": explained,
    }

    # Positive and negative assertions: genuine separation survives; a pretty partition of a null is not stable evidence.
    assert blob_original_sil > 0.45
    assert min(blob_map_sil.values()) > 0.25
    assert null_tsne_sil > null_original_sil + 0.12
    assert tsne_cluster_ari < 0.65

    np.savez_compressed(
        output / "projection_claim_lab_arrays.npz",
        bridge_x=bridge_x,
        bridge_t=bridge_t,
        blobs_x=blobs_x,
        blobs_label=blobs_label,
        null_x=null_x,
        cases_x=cases_x,
        **{f"bridge_{k}": v for k, v in bridge_maps.items()},
        **{f"blobs_{k}": v for k, v in blob_maps.items()},
        **{f"null_{k}": v for k, v in null_maps.items()},
        **{f"cases_{k}": v for k, v in case_maps.items()},
        **{f"case_{k}": v for k, v in case_meta.items()},
    )
    (output / "projection_claim_lab_results.json").write_text(json.dumps(results, indent=2), encoding="utf-8")
    return results


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output", type=Path, default=Path("projection_claim_lab_output"))
    parser.add_argument("--quick", action="store_true", help="Use fewer synthetic points")
    args = parser.parse_args()
    results = run_lab(args.output, quick=args.quick)
    print(json.dumps(results, indent=2))


if __name__ == "__main__":
    main()

What would weaken the central thesis?

The thesis would weaken if a broad family of practically relevant datasets admitted two-dimensional embeddings that preserved all decision-relevant distances, neighbourhoods, densities and topology across methods without a measurable trade-off. It would also weaken if visual cluster judgements remained calibrated under credible nulls and irrelevant interventions. Existing theory and empirical results point the other way, but the claim is bounded: a particular dataset can be intrinsically two-dimensional, nearly linear or cleanly separated, as the positive control shows.

An open hypothesis follows. Projection-risk may be predicted better by a stability surface across metric, sampling and scale than by the highest quality score from one configuration. The lab provides the beginning of such an instrument, but a cross-domain study would be needed before treating it as established.

stateDiagram-v2
  [*] --> Claim
  state "PCA, preserve a linear shadow" as Claim
  state "t-SNE, preserve probable<br/>neighbours" as Test
  state "UMAP, preserve a fuzzy local graph" as Revision
  Claim --> Test: expose mechanism
  Test --> Decision: evidence survives
  Test --> Revision: rival survives
  Revision --> Claim: narrow or rebuild
  Decision --> [*]

Compact glossary

Projection artefact

A visible relation caused or materially shaped by the reduction pipeline rather than the proposed generating phenomenon.

Reference geometry

The features, scaling, metric and temporal interpretation against which preservation is judged.

Trustworthiness

How rarely the map introduces neighbours that were not near in the reference space.

Continuity

How rarely the map loses neighbours that were near in the reference space.

Manifold claim

A hypothesis that observations lie near a lower-dimensional structure with meaningful local geometry or topology.

Stability surface

The behaviour of a named claim across seeds, samples, metrics, scales and other controlled interventions.

Release the claim, not the picture

A two-dimensional map is useful precisely because it forgets. PCA forgets what lies outside its best linear plane. t-SNE spends most of its attention on probabilistic neighbours. UMAP turns locally scaled neighbour relations into a fuzzy graph and tries to reproduce that graph. Their differences are not defects to be averaged away. They are clues to what each picture can support.

The durable decision is to specify the preserved relation before choosing the projection. Then test that relation in the reference space, measure false and missing neighbours at several scales, vary sampling and optimisation, compare a simple baseline, and require an external observation before naming a cluster, trajectory or manifold. A beautiful map that fails these tests remains valuable as a hypothesis generator. It does not become a population model, a scientific mechanism or an operational segmentation.

For architecture and governance, this means two-dimensional coordinates should not silently become world state or decision authority. Store the feature and metric contract, projection configuration, seed, sample definition, quality curves, intervention results and the exact claim made. Preserve the original evidence so the claim can be replayed when data or policy changes.

The question after seeing an island is therefore no longer “What group is this?” It is “Which relationship made this island cheap to draw, and does that relationship survive outside the drawing?” That question turns dimensionality reduction from persuasive geometry into a disciplined measurement system.