UMAP

With default local connectivity and bandwidth settings, the calibration searches for a scale whose retained directed weights sum to $\log_2(k)$. In the umap-learn training convention reproduced here, the stored list has $k$ entries including the anchor itself. With $k=4$, three other candidates remain, and their target sum is $\log_2(4)=2$.

The search changes $\sigma_i$ until it reaches that target within numerical tolerance. Our five distinct inputs make this calculation straightforward. Duplicate points, tied nearest distances and numerical bandwidth floors require additional handling in a general implementation. The C++ companion states its narrower assumptions explicitly.

import numpy as np
from scipy.optimize import brentq
from sklearn.neighbors import NearestNeighbors

x = np.array([0, 1, 3, 6, 10.])
k = 4  # The stored training list includes self; three other candidates remain.
distances, indices = NearestNeighbors(n_neighbors=k).fit(x[:, None]).kneighbors(x[:, None])
rho = distances[:, 1]  # Nearest other point; all inputs here are distinct.
sigma = np.zeros(len(x))
W = np.zeros((len(x), len(x)))
for i in range(len(x)):
    excess = np.maximum(0, distances[i, 1:] - rho[i])
    # The bracket and target are suitable for this small example.
    sigma[i] = brentq(lambda s: np.exp(-excess/s).sum() - np.log2(k), 1e-6, 100)
    W[i, indices[i, 1:]] = np.exp(-excess/sigma[i])
U = W + W.T - W*W.T  # Elementwise product, not matrix multiplication.
print(W[0])  # About [0, 1, 0.654045, 0.345955, 0].
print(U[0, 2])  # About 0.827023: combined A--C weight.
np.testing.assert_allclose(W.sum(axis=1), np.log2(k))
np.testing.assert_allclose(U, U.T)

Let $y_i,y_j$ be two map positions. A smooth curve converts their distance into an embedding strength $v_{ij}$. Positive shape parameters $a,b$ are fitted from min_dist and spread:

There is no global normalization of these strengths. For a fixed input edge strength $u$ and proposed map strength $v$, the binary cross-entropy term is

A strong input edge puts more weight on the attractive part, $-\log v$. A weak input edge puts more weight on the repulsive part, $-\log(1-v)$. If a pair’s strength were free to vary independently in this full objective, its optimum would be $v=u$. Real map positions couple all the pairs.

The usual implementation samples positive edges and negative vertices instead of enumerating all pairs. Its effective loss differs from the full cross-entropy expression; it is not simply an unbiased shortcut to that same objective. Damrich and Hamprecht (2021) analyze the resulting weighting. The formula here explains the attraction/repulsion idea, rather than claiming that every fitted $v$ should equal its input $u$.

import numpy as np
import umap
from sklearn.datasets import make_blobs

centers = np.zeros((3, 8))
centers[:, :2] = [[-6, 0], [0, 6], [6, 0]]
X, source_group = make_blobs(
    n_samples=[60, 90, 120], centers=centers,
    cluster_std=[0.3, 0.8, 1.3], random_state=17,
)
# These synthetic features already have comparable units.
embeddings, graphs = {}, {}
for neighbors in [5, 30]:
    for min_dist in [0.05, 0.6]:
        model = umap.UMAP(
            n_neighbors=neighbors, min_dist=min_dist,
            n_components=2, metric="euclidean", random_state=42, n_jobs=1,
        ).fit(X)  # No labels supplied to this unsupervised fit.
        embeddings[neighbors, min_dist] = model.embedding_.copy()
        if neighbors in graphs:
            np.testing.assert_allclose(graphs[neighbors], model.graph_.toarray())
        graphs[neighbors] = model.graph_.toarray()

These figures were generated with umap-learn 0.5.12. The dense graph copies are only a check for this small dataset; retain the sparse representation when working with larger data.

from sklearn.model_selection import train_test_split

# X is the synthetic feature matrix from the preceding example.
X_train, X_new = train_test_split(X, test_size=0.2, random_state=7)
mapper = umap.UMAP(n_neighbors=15, min_dist=0.1,
                   random_state=42, transform_seed=42, n_jobs=1).fit(X_train)
Z_train = mapper.embedding_.copy()
Z_new = mapper.transform(X_new)
np.testing.assert_array_equal(mapper.embedding_, Z_train)
print(Z_train.shape, Z_new.shape)  # (216, 2) (54, 2).

The code leaves labels out of UMAP.fit. UMAP also supports supervised fitting when labels are supplied. In a classifier pipeline that forwards y to intermediate estimators, be explicit about whether that supervised representation is intended.

#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <utility>
#include <vector>
using Matrix = std::vector<std::vector<double>>;
struct Graph { Matrix directed, combined; std::vector<double> rho, sigma; };

Graph build_graph(const std::vector<double>& x, std::size_t k) {
    const std::size_t n = x.size();
    if (k < 3 || k > n) throw std::invalid_argument("need 3 <= k <= row count");
    for (double v : x) if (!std::isfinite(v))
        throw std::invalid_argument("finite scalar inputs required");
    Graph graph{Matrix(n, std::vector<double>(n)), Matrix(n, std::vector<double>(n)),
                std::vector<double>(n), std::vector<double>(n)};
    for (std::size_t i = 0; i < n; ++i) {
        std::vector<std::pair<double, std::size_t>> neighbors;
        for (std::size_t j = 0; j < n; ++j) if (i != j) {
            double distance = std::abs(x[i]-x[j]);
            if (!std::isfinite(distance)) throw std::overflow_error("distance overflow");
            if (distance == 0) throw std::invalid_argument("distinct inputs required");
            neighbors.emplace_back(distance, j);
        }
        std::sort(neighbors.begin(), neighbors.end());
        neighbors.resize(k-1);  // Self occupies the omitted stored-list entry.
        double rho = neighbors[0].first;
        if (neighbors[1].first == rho)
            throw std::invalid_argument("this example requires a unique nearest neighbor");
        double spread = neighbors.back().first-rho;
        // Solve in relative units so the bisection is insensitive to scalar units.
        auto mass = [&](double s) {
            double sum = 0;
            for (auto [distance, j] : neighbors) {
                (void)j;
                sum += std::exp(-((distance-rho)/spread)/s);
            }
            return sum;
        };
        double target = std::log2(static_cast<double>(k)), low = 0, high = 1;
        while (mass(high) < target) high *= 2;
        for (int step = 0; step < 80; ++step) {
            double middle = (low+high)/2;
            if (mass(middle) < target) low = middle; else high = middle;
        }
        double s = (low+high)/2, sigma = spread*s;
        if (!(sigma > 0) || !std::isfinite(sigma))
            throw std::overflow_error("unrepresentable bandwidth");
        graph.rho[i] = rho; graph.sigma[i] = sigma;
        for (auto [distance, j] : neighbors)
            graph.directed[i][j] = std::exp(-((distance-rho)/spread)/s);
    }
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            double a = graph.directed[i][j], b = graph.directed[j][i];
            graph.combined[i][j] = a+b-a*b;
        }
    return graph;
}

Call build_graph({0, 1, 3, 6, 10}, 4). The directed A row is approximately [0, 1, 0.654045, 0.345955, 0], and the combined A–C strength is about 0.827023. This dense, exact-neighbor implementation makes the steps inspectable; the library uses more general neighbor handling and sparse graphs.