The Atlas-Learn Approach to the Manifold Hypothesis

R-bloggers 2026-05-20

[This article was first published on R Works, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The manifold hypothesis, the idea that real-world high-dimensional data concentrates near a low-dimensional curved subspace, is foundational to modern machine learning. Many popular manifold learning methods such as UMAP, t-SNE, Isomap, and diffusion maps do achieve dimensionality reduction by embedding data into a flat Euclidean space $\mathbb{R}^D$ , but they do not attempt to directly learn the underlying manifold. In contrast, the 2025 paper by Robinett et al., Atlas-based manifold representations for interpretable Riemannian machine learning, offers a proof of concept for directly tackling the manifold hypothesis based on fundamental ideas from differential geometry. It provides an algorithm for learning a low dimensional manifold from point data by constructing an atlas of charts. The paper is also notable for the design of an efficient data structure for working with the learned atlas and for the extensive supplementary materials that include a GitHub Repository containing several practical Python algorithms for doing calculations on manifolds, and an extraordinary amount of implementation detail.

Reading through Robinett et al., however, requires a fairly deep background in the theory of differential geometry. This post is an attempt to provide an on-ramp to Robinett et al. by discussing the relatively simple example of the two dimensional sphere, $S^2$ embedded in $\mathbb{R^3}$ . It implements the Atlas-Learn data structures and algorithms in R, uses them to learn $S^2$ and then goes on to validate the Atlas-Learn algorithm for the sphere via three independent methods: 1) use numerical integration along the manifold to trace a great circle on the sphere, 3) recover the radius of curvature of the sphere from the atlas, and 4) verify the Gauss-Bonnet Theorem for the sphere.

The R code was mostly worked out by Claude Sonnet 4.3 in the context of participating in the Posit beta test for its AI Assistant. I found the integration of the AI engine into the RStudio IDE an effective means of communicating with Claude and managing the project workflow.

Atlas-Learn: Theory and Algorithm

This section provides some minimal theoretical background for understanding the Atlas-Learn algorithm. A smooth manifold $\mathcal{M}$ of intrinsic dimension $d$ embedded in $\mathbb{R}^D$ can be described by an atlas — a finite collection of charts $\{(\varphi_j, \Omega_j)\}_{j=1}^{k}$ such that the open sets $U_j = \varphi_j^{-1}(\Omega_j)$ cover $\mathcal{M}$ and each chart map $\varphi_j : U_j \to \mathbb{R}^d$ is a smooth bijection onto its image.

Normally, the definition of a smooth manifold also requires that any two charts be smoothly compatible, where two charts $(\mathcal{U}_1, \phi_1)$ and $(\mathcal{U}_2, \phi_2)$ are said to be iff $\phi_1(\mathcal{U}_1 \cap \mathcal{U}_2)$ and $\phi_2(\mathcal{U}_1 \cap \mathcal{U}_2)$ are both open in $\mathbb{R}^m$ and the transition map $\phi_{21} = \phi_2 \circ \phi_1^{-1} : \phi_1(\mathcal{U}_1 \cap \mathcal{U}_2) \to \phi_2(\mathcal{U}_1 \cap \mathcal{U}_2)$ is a diffeomorphism (e.g. see [2]). Robinett et al. relax the smoothly compatible requirement and define transition maps $\psi_{ij}$ separately from coordinate chart images. They then approximate a differentiable atlas by ensuring that the discrepancy between coordinate charts and transition maps $||{\phi(\xi) - \phi(\psi(xi)))}||_\mathbb{R^D}$ goes to 0 as the number of charts and the number of points sampled goes to infinity.

In the Atlas-Learn algorithm the manifold is a surface ( $d = 2$ ) embedded in $\mathbb{R}^3$ , and both the covering sets and the chart maps are learned from a finite point cloud $\{x_i\}_{i=1}^{N} \subset \mathcal{M}$ . The algorithm proceeds to construct an atlas in four basic steps.

The Atlas-Learn algorithm proceeds in four steps for each chart:

The point cloud comprising the data, the sphere in our case, is partitioned into k-medoids.
Local PCA is used to find the tangent plane and the normal plane for each point.
Quadratic regression is performed to find the curvature coefficients, K
The minimum ellipsoidal region enclosing the chart is estimated.

Step 1: Partitioning via k-medoids

The point cloud is partitioned into $k$ clusters using the $k$ -medoids algorithm (PAM). Unlike $k$ -means, PAM selects actual data points as cluster centers (medoids), which makes the partition robust to outliers and avoids projection artefacts. Each point $x_i$ receives a chart label $j(i) \in \{1,\ldots,k\}$ , and the points belonging to chart $j$ together with their centroid are

$\mathcal{X}_j = \{ x_i : j(i) = j \}, \qquad \bar{x}_j = \frac{1}{|\mathcal{X}_j|} \sum_{x \in \mathcal{X}_j} x.$

Step 2: Local PCA and tangent-plane estimation

For each cluster $j$ , the centered data matrix $\widetilde{X}_j = \mathcal{X}_j - \bar{x}_j \in \mathbb{R}^{|\mathcal{X}_j| \times 3}$ is decomposed via the thin SVD:

$\widetilde{X}_j = U \Sigma V^\top, \qquad V = [v_1 \mid v_2 \mid v_3].$

The first two right singular vectors span the local tangent plane:

$L_j = [v_1 \mid v_2] \in \mathbb{R}^{3 \times 2},$

while the third singular vector $m_j = v_3 \in \mathbb{R}^3$ estimates the local surface normal (the direction of least variance). Each centered point is then decomposed into tangent and normal components:

$\tau_i = \widetilde{x}_i^\top L_j \in \mathbb{R}^2, \qquad \nu_i = \widetilde{x}_i^\top m_j \in \mathbb{R}.$

Step 3: Quadratic chart map

On a smooth surface the normal offset $\nu$ is a smooth function of the tangent coordinates $\tau$ . Atlas-Learn approximates this by a degree-2 polynomial (capturing local curvature):

$\nu \approx K_j \, \phi(\tau), \qquad \phi(\tau) = \bigl[\tau_1^2,\; \tau_1\tau_2,\; \tau_2^2\bigr]^\top,$

where $K_j \in \mathbb{R}^{1 \times 3}$ is estimated by ordinary least squares with a small ridge penalty $\varepsilon = 10^{-10}$ :

$K_j = \boldsymbol{\nu}^\top \Phi \bigl(\Phi^\top \Phi + \varepsilon I_3\bigr)^{-1}, \qquad \Phi = \bigl[\phi(\tau_1) \;\cdots\; \phi(\tau_{|\mathcal{X}_j|})\bigr]^\top \in \mathbb{R}^{|\mathcal{X}_j| \times 3}.$

The resulting inverse chart map $\varphi_j^{-1} : \mathbb{R}^2 \to \mathbb{R}^3$ reconstructs an ambient point from local coordinates $\xi \in \mathbb{R}^2$ :

$\varphi_j^{-1}(\xi) \;=\; \bar{x}_j + L_j\,\xi + m_j\,K_j\,\phi(\xi).$

Its Jacobian $J_j(\xi) = \partial\varphi_j^{-1}/\partial\xi \in \mathbb{R}^{3 \times 2}$ , required for geodesic integration, is:

$J_j(\xi) \;=\; L_j \;+\; m_j\,K_j\,\frac{\partial\phi}{\partial\xi}(\xi), \qquad \frac{\partial\phi}{\partial\xi}(\xi) = \begin{pmatrix} 2\xi_1 & 0 \\ \xi_2 & \xi_1 \\ 0 & 2\xi_2 \end{pmatrix}.$

Step 4: Ellipsoidal chart domains

Each chart is assigned an ellipsoidal domain $\Omega_j \subset \mathbb{R}^2$ defined by

$\Omega_j = \bigl\{\, \xi \in \mathbb{R}^2 : \xi^\top A_j\,\xi \leq 1 \,\bigr\},$

where $A_j$ is a rescaled inverse covariance of the projected points:

$A_j = \frac{\operatorname{Cov}(\tau)^{-1}} {\rho \cdot \max_i \bigl\{\tau_i^\top\operatorname{Cov}(\tau)^{-1}\tau_i\bigr\}}, \qquad \rho = \texttt{ellipsoid\_scale}.$

Setting $\rho > 1$ (default $\rho = 1.1$ ) inflates each domain slightly beyond the convex hull of its own cluster, so that neighboring charts overlap and transitions are always possible. On $S^2$ specifically, because the sphere is isotropic and the $k$ -medoids partition tends to produce roughly equal-area, near-circular patches, the learned ellipsoids $\Omega_j$ are close to circles ( $A_j \approx \lambda I_2$ for some scalar $\lambda$ ). Each chart is assigned a domain $\Omega_j = \{\xi \in \mathbb{R}^2 : \xi^\top A_j\xi \leq 1\}$ where $A_j$ is a rescaled inverse covariance of the projected tangent-plane coordinates. Setting the scale factor $\rho > 1$ (default $\rho = 1.1$ ) inflates domains slightly so that neighboring charts overlap and transitions are always possible.

The Atlas-Learn Implementation

Here are the required R packages.

Required Packages

library(tidyverse)library(cluster)    # pam() for k-medoids partitioninglibrary(RANN)       # nn2() for fast k-nearest-neighbor querieslibrary(plotly)     # interactive 3D visualizationlibrary(purrr)      # map() / imap() for list operations

This block of code contains all of the functions for the Atlas Learn implementation.

Show the code

#| label: atlas-functions# ===========================================================================# PART 1: Quadratic feature helpers# ---------------------------------------------------------------------------# These functions implement the d=2 specialization of the general quadratic# feature map. For general d, phi(xi) would have choose(d+1, 2) components.# For d=2 it has exactly 3: [xi1^2, xi1*xi2, xi2^2].# ===========================================================================# Maps xi in R^2 to the three quadratic monomials used to model surface curvature.# General d would give choose(d+1,2) monomials; for d=2 this is exactly 3.quad_features <- function(xi) {  c(xi[1]^2, xi[1] * xi[2], xi[2]^2)}# Jacobian of quad_features w.r.t. xi: a 3x2 matrix (d=2 specialization).# Row i is the gradient of the i-th monomial:#   d/dxi [xi1^2]     = [2*xi1,   0   ]#   d/dxi [xi1*xi2]   = [xi2,     xi1 ]#   d/dxi [xi2^2]     = [0,       2*xi2]quad_jacobian <- function(xi) {  matrix(    c(2 * xi[1],   xi[2],         0,              0,   xi[1],  2 * xi[2]),    nrow = 3, ncol = 2, byrow = TRUE  )}# ===========================================================================# PART 2: Core chart operations# ---------------------------------------------------------------------------# Each chart stores:#   mean  : R^3 — chart center (centroid of the cluster)#   L     : R^{3x2} — orthonormal tangent basis (columns = v1, v2 from SVD)#   M     : R^{3x1} — unit surface normal (v3 from SVD; scalar normal because D-d=1)#   K     : R^{1x3} — quadratic curvature coefficients (1 row because D-d=1)#   ell_A : R^{2x2} — ellipsoidal domain matrix (2x2 because d=2)# ===========================================================================# Evaluate the inverse chart map: xi in R^2  ->  point in R^3.# Formula: x = mean + L*xi + M * (K * phi(xi))#   - mean + L*xi : linear move along the tangent plane#   - M*(K*phi)   : quadratic normal correction encoding surface curvaturechart_eval <- function(chart, xi) {  q <- quad_features(xi)          # 3-vector: [xi1^2, xi1*xi2, xi2^2]  as.vector(    chart$mean +    chart$L %*% xi +              # 3x2 * 2x1 = 3x1 tangent contribution    chart$M * as.numeric(chart$K %*% q)   # 3x1 * scalar normal correction  )}# Jacobian of the inverse chart map at xi: a 3x2 matrix.# J(xi) = L + M * (K * dQ(xi))#   dQ(xi) is the 3x2 Jacobian of the monomial map (quad_jacobian).#   K * dQ  is 1x3 * 3x2 = 1x2; then M * (K*dQ) is 3x1 * 1x2 = 3x2.# This is the key object for geodesic integration: it maps tangent vectors# in R^2 chart coordinates to ambient R^3 velocity vectors.chart_jacobian <- function(chart, xi) {  dQ <- quad_jacobian(xi)                   # 3x2 monomial Jacobian  chart$L + chart$M %*% (chart$K %*% dQ)   # 3x2 total Jacobian J(xi)}# Project an ambient R^3 point onto this chart's tangent-plane coordinates (R^2).# The linear projection xi = L^T * (x - mean) is exact for points on the tangent# plane; for points on the actual surface it is a first-order approximation.chart_project <- function(chart, x) {  as.vector(t(chart$L) %*% (x - chart$mean))}# Test whether tangent coords xi lie within the chart's ellipsoidal domain.# The domain is {xi in R^2 : xi^T * A * xi <= 1}, a 2D Mahalanobis ball.# This is O(d^2) = O(4) for d=2.chart_in_domain <- function(chart, xi) {  as.numeric(t(xi) %*% chart$ell_A %*% xi) <= 1}# ===========================================================================# PART 3: atlas_learn() — the main learning function# ---------------------------------------------------------------------------# Input : X (N x 3 matrix of surface points), k (number of charts)# Output: an S3 object of class "atlas" containing k atlas_chart objects## The algorithm runs four steps per chart:#   (a) k-medoids partitioning#   (b) local PCA to find the tangent plane and normal#   (c) quadratic regression for the curvature coefficients K#   (d) ellipsoidal domain construction## d=2, D=3 specializations that are hard-coded here:#   - SVD takes nv=3 (the full 3x3 right-singular-vector matrix)#   - L = V[,1:2] is 3x2; M = V[,3] is 3x1 (one normal vector, not a matrix)#   - nu (normal offset) is a scalar per point, not a vector#   - K is fitted as a 1x3 row vector (one row per normal direction)#   - ell_A is 2x2 (domain is a 2D ellipse)# ===========================================================================atlas_learn <- function(X, k, ellipsoid_scale = 1.1) {  # X              : N x 3 matrix of surface points  # k              : number of charts  # ellipsoid_scale: inflate domains by this factor so adjacent charts overlap  message("Fitting k-medoids (k = ", k, ") ...")  # PAM (Partitioning Around Medoids) is preferred over k-means for surfaces:  # medoids are actual data points, making the partition robust to outliers.  km <- cluster::pam(X, k, diss = FALSE)  charts <- seq_len(k) |>    purrr::map(\(j) {      # --- (a) Extract the cluster and center it ----------------------------      idx <- which(km$clustering == j)      Xj  <- X[idx, , drop = FALSE]   # N_j x 3 matrix of cluster points      m   <- colMeans(Xj)             # chart center (3-vector)      Xc  <- sweep(Xj, 2, m)         # centered cluster: N_j x 3      # --- (b) Local PCA: tangent plane and normal -------------------------      # SVD of the centered cluster reveals local geometry:      #   - first two right singular vectors (v1, v2) span the tangent plane      #   - third right singular vector (v3) is the surface normal      # d=2, D=3: we always take nv=3 because D=3 (fully determined).      sv  <- svd(Xc, nu = 0, nv = 3)      # L: 3x2 tangent basis (orthonormal columns). General form: R^{D x d}.      L   <- sv$v[, 1:2, drop = FALSE]      # M: 3x1 unit normal. d=2, D=3 specialization: D-d=1, so there is      # exactly ONE normal direction and M is a column vector, not a matrix.      M   <- sv$v[, 3, drop = FALSE]      # Project each centered point into tangent / normal coordinates.      # tau: N_j x 2 — tangent coordinates (2D because d=2)      # nu:  N_j x 1 — normal offsets (scalar because D-d=1)      tau <- Xc %*% L    # N_j x 2      nu  <- Xc %*% M    # N_j x 1 (scalar per point; would be N_j x (D-d) in general)      # --- (c) Quadratic curvature regression ------------------------------      # Fit: nu ~ K * phi(tau)  where phi(tau) = [tau1^2, tau1*tau2, tau2^2]      # K is 1x3 here (one output because D-d=1; would be (D-d)x3 in general).      #      # Design matrix Q: N_j x 3, each row is phi(tau_i)      Q   <- t(apply(tau, 1, quad_features))   # N_j x 3      # Ridge-regularized least squares: K = nu^T * Q * (Q^T*Q + esp*I)^{-1}      # The ridge penalty eps=1e-10 prevents singularity when clusters are      # nearly collinear in tangent coordinates.      K   <- t(nu) %*% Q %*% solve(crossprod(Q) + 1e-10 * diag(3))  # 1 x 3      # --- (d) Ellipsoidal domain ------------------------------------------      # Domain: {xi in R^2 : xi^T A xi <= 1}, a Mahalanobis ball.      # A_raw = Cov(tau)^{-1}: inverse covariance of the tangent coordinates.      # This adapts the domain shape to the data spread in each direction.      A_raw <- solve(cov(tau))      # Scale A so the outermost point lands on the boundary, then inflate      # by ellipsoid_scale (default 1.1) to ensure overlap with neighbors.      qvals <- apply(tau, 1, \(xi) as.numeric(t(xi) %*% A_raw %*% xi))      ell_A <- A_raw / (ellipsoid_scale * max(qvals))   # 2x2 positive-definite      # Pack all chart data into an S3 object      structure(        list(mean = m, L = L, M = M, K = K, ell_A = ell_A, idx = idx,             n_points = length(idx)),        class = "atlas_chart"      )    })  # Return the full atlas as an S3 object  structure(    list(      charts        = charts,      k             = k,      clustering    = km$clustering,      n_points      = nrow(X),      ambient_dim   = ncol(X),   # D = 3      intrinsic_dim = 2L         # d = 2 (hard-coded for surfaces in R^3)    ),    class = "atlas"  )}# ===========================================================================# PART 4: S3 print / summary methods# ===========================================================================print.atlas_chart <- function(x, ...) {  cat(sprintf(    "<atlas_chart>  %d points | mean [%s] | cond(ell_A) = %.1f\n",    x$n_points,    paste(round(x$mean, 3), collapse = ", "),    kappa(x$ell_A)  ))  invisible(x)}# Returns a per-chart summary tibblesummary.atlas <- function(object, ...) {  purrr::imap_dfr(object$charts, \(ch, i)    tibble::tibble(      chart      = i,      n_points   = ch$n_points,      mean_norm  = round(sqrt(sum(ch$mean^2)), 4),      cond_ell_A = round(kappa(ch$ell_A), 1)    )  )}print.atlas <- function(x, ...) {  cat(sprintf(    "<atlas>  k = %d | ambient R^%d | intrinsic R^%d | %d points\n\n",    x$k, x$ambient_dim, x$intrinsic_dim, x$n_points  ))  print(summary(x))  invisible(x)}# ===========================================================================# PART 5: Chart lookup# ===========================================================================# Return the index of the nearest chart whose domain contains x.# Falls back to the globally nearest chart center if none qualify.# Checking only the 6 nearest charts (by Euclidean center distance) is# sufficient in practice and avoids an O(k) domain test at every step.find_chart <- function(atlas, x) {  dists      <- map_dbl(atlas$charts, \(ch) sum((x - ch$mean)^2))  candidates <- order(dists)[seq_len(min(6L, atlas$k))]  in_domain  <- purrr::keep(    candidates,    \(i) chart_in_domain(atlas$charts[[i]],                         chart_project(atlas$charts[[i]], x))  )  if (length(in_domain) > 0L) in_domain[[1L]] else which.min(dists)}# ===========================================================================# PART 6: Quasi-Euclidean retraction# ---------------------------------------------------------------------------# Advances one step from (chart_idx, xi) in the ambient direction tau_r3 (R^3).## The pseudoinverse J^+ = (J^T J)^{-1} J^T is the key d=2, D=3 quantity:#   - J is 3x2 (D x d)#   - J^T J is 2x2 (d x d) — the pullback metric g; invertible for full-rank J#   - J^+ is 2x3 (d x D) — pulls ambient vectors back to chart coordinates## In general D > d, J^+ is the minimum-norm left inverse of J.# ===========================================================================retract_step <- function(atlas, chart_idx, xi, tau_r3) {  chart  <- atlas$charts[[chart_idx]]  J      <- chart_jacobian(chart, xi)        # 3x2 (D x d for d=2, D=3)  Jp     <- solve(crossprod(J)) %*% t(J)    # 2x3 pseudoinverse: (J^T J)^{-1} J^T  xi_new <- xi + as.vector(Jp %*% tau_r3)   # advance in chart coordinates  if (chart_in_domain(chart, xi_new)) {    # Still inside the same chart: evaluate and return    return(list(      chart_idx = chart_idx,      xi        = xi_new,      x         = chart_eval(chart, xi_new)    ))  }  # Step crossed a chart boundary: find a new host chart and re-project  x_cand <- chart_eval(chart, xi_new)       # approximate ambient position  ci_new <- find_chart(atlas, x_cand)       # nearest chart that contains x_cand  xi2    <- chart_project(atlas$charts[[ci_new]], x_cand)   # re-project  list(    chart_idx = ci_new,    xi        = xi2,    x         = chart_eval(atlas$charts[[ci_new]], xi2)  )}# ===========================================================================# PART 7: Geodesic path integration# ---------------------------------------------------------------------------# Traces a geodesic from x_start in direction direction_r3 (ambient R^3).# Returns a tibble with columns x, y, z, step, chart_idx.## At each step the ambient velocity vector is re-projected into the current# chart's tangent plane and pushed back to ambient.  This keeps the direction# of motion consistent across chart boundaries — without it, the fixed tau# vector drifts after every chart transition and the path wanders off the# intended geodesic.# ===========================================================================atlas_geodesic <- function(atlas, x_start, direction_r3,                           n_steps = 100L, step_size = 0.02) {  ci  <- find_chart(atlas, x_start)  ch  <- atlas$charts[[ci]]  xi  <- chart_project(ch, x_start)  # Project the initial direction onto the starting chart's tangent plane,  # then normalize to unit speed.  This ensures tau_r3 is always in T_x M.  J      <- chart_jacobian(ch, xi)  Jp     <- solve(crossprod(J)) %*% t(J)          # 2x3 pseudoinverse  tau_r3 <- as.vector(J %*% (Jp %*% direction_r3))  # project onto tangent plane  tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size  # unit-speed scaling  steps     <- vector("list", n_steps + 1L)  chart_ids <- integer(n_steps + 1L)  steps[[1L]]     <- chart_eval(ch, xi)  chart_ids[[1L]] <- ci  for (i in seq_len(n_steps)) {    res    <- retract_step(atlas, ci, xi, tau_r3)    ci     <- res$chart_idx    xi     <- res$xi    steps[[i + 1L]]     <- res$x    chart_ids[[i + 1L]] <- ci    # Re-project tau_r3 into the current chart's tangent plane and normalize.    # This is the identity-transport approximation: it holds the direction    # constant in the current chart's frame rather than applying the    # Christoffel correction that true parallel transport would require.    J_new  <- chart_jacobian(atlas$charts[[ci]], xi)    Jp_new <- solve(crossprod(J_new)) %*% t(J_new)    tau_r3 <- as.vector(J_new %*% (Jp_new %*% tau_r3))    tau_r3 <- tau_r3 / sqrt(sum(tau_r3^2)) * step_size    J      <- J_new  }  do.call(rbind, steps) |>    as_tibble(.name_repair = "minimal") |>    set_names(c("x", "y", "z")) |>    mutate(step = row_number(), chart_idx = chart_ids)}# ===========================================================================# PART 8: Sphere-specific helpers# ---------------------------------------------------------------------------# These are ground-truth functions for S^2 used only for validation;# they play no role in the Atlas-Learn algorithm itself.# ===========================================================================# Great-circle geodesic distance between two unit vectors (in radians).# d(p, q) = arccos(p . q), clamped to [-1, 1] to avoid NaN from floating point.sphere_dist <- function(p, q) {  acos(pmax(pmin(sum(p * q), 1.0), -1.0))}# Unit tangent vector at p pointing toward q along the great circle.# Computed by projecting q onto the tangent plane at p and normalizing.geodesic_direction_sphere <- function(p, q) {  v   <- q - sum(p * q) * p   # remove radial component  nrm <- sqrt(sum(v^2))  if (nrm < 1e-12) return(NULL)  v / nrm}# Run the Atlas geodesic from p toward q for the exact number of steps# needed to cover the true great-circle distance, then measure endpoint error.atlas_endpoint_error <- function(atlas, p, q, step_size = 0.02) {  d   <- sphere_dist(p, q)  dir <- geodesic_direction_sphere(p, q)  if (d < 1e-6 || is.null(dir)) {    return(tibble(true_dist = d, endpoint_error = 0, n_steps = 0L))  }  n_steps  <- max(1L, round(d / step_size))  path     <- atlas_geodesic(atlas, p, dir,                              n_steps = n_steps, step_size = step_size)  endpoint <- as.numeric(path[nrow(path), c("x", "y", "z")])  # Re-project endpoint onto S^2 before measuring angle error.  # This separates geodesic-direction error from on-manifold drift.  endpoint <- endpoint / sqrt(sum(endpoint^2))  tibble(    true_dist      = d,    endpoint_error = sphere_dist(endpoint, q),    n_steps        = as.integer(n_steps)  )}

Construct the Atlas

This section of code constructs the sphere point cloud by sampling 2,000 points uniformly from $S^2$ using the standard parametrization $(\theta, \phi)$ with $\cos\phi$ drawn uniformly on $[-1,1]$ , which avoids the pole-crowding that arises from drawing $\phi$ uniformly.

Show the code

#| code-summary: "Show the code"#| label: generate-sphereset.seed(4471)N     <- 2000Ltheta <- runif(N, 0, 2 * pi)phi   <- acos(runif(N, -1, 1))   # uniform area measure on S^2sphere_pts <- tibble(  x = sin(phi) * cos(theta),  y = sin(phi) * sin(theta),  z = cos(phi))X_sphere <- as.matrix(sphere_pts)

This block of code fits an atlas with $k = 25$ charts to the sampled sphere. The fit is cached to atlas_s2_cache.rds to avoid unnecessarily repeating the this time consuming process during development sessions.

Show the code

#| label: learn-atlasatlas_cache <- "atlas_s2_cache.rds"if (file.exists(atlas_cache)) {  message("Loading cached atlas from ", atlas_cache)  atlas_s2 <- readRDS(atlas_cache)} else {  message("Cache not found — fitting atlas (k = 25) ...")  atlas_s2 <- atlas_learn(X_sphere, k = 25L)  saveRDS(atlas_s2, atlas_cache)  message("Atlas saved to ", atlas_cache)}# Attach chart labels to the point-cloud tibble for visualizationsphere_pts <- sphere_pts |>  mutate(chart = factor(atlas_s2$clustering))

The atlas data structure

In this section we provide an overview of the structure of the Atlas. Theatlas_learn() function returns an S3 object of class "atlas" containing a nested list of 25 "atlas_chart" objects, one per cluster.

atlas_s2├── k              integer          number of charts (25)├── n_points       integer          total data points (2000)├── ambient_dim    integer          dimension of embedding space (3 = D)├── intrinsic_dim  integer          dimension of the manifold (2 = d)├── clustering     integer[2000]    chart label for every data point└── charts         list[25]         one atlas_chart per cluster    └── charts[[j]]        ├── mean     numeric[3]     centroid of cluster j in R³        ├── L        numeric[3×2]   orthonormal tangent basis (columns = v₁, v₂)        ├── M        numeric[3×1]   unit surface normal = v₃  [For d=2: single vector]        ├── K        numeric[1×3]   quadratic curvature coefficients  [for d=2: 1 row]        ├── ell_A    numeric[2×2]   ellipsoidal domain matrix  [d=2 specialisation: 2×2]        ├── idx      integer[n_j]   row indices into the N×3 data matrix        └── n_points integer        cluster size n_j

Here are the values for each field of the $k = 25$ charts.

Show the code

purrr::imap_dfr(atlas_s2$charts, function(ch, j) {  tibble(    Chart  = j,    `n pts` = ch$n_points,    `centre (x,y,z)` = paste(round(ch$mean, 2), collapse = ", "),    # For the unit sphere, M should align with the radial direction (M·r̂ ≈ 1)    `M · r̂`  = round(abs(sum(ch$M * ch$mean / sqrt(sum(ch$mean^2)))), 3),    # Quadratic curvature coefficients    k1 = round(ch$K[1], 3),    k2 = round(ch$K[2], 3),    k3 = round(ch$K[3], 3),    # Semi-axes of the 2D ellipsoidal domain (1/sqrt of eigenvalues of ell_A)    `ell axis a` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[2]), 3),    `ell axis b` = round(1 / sqrt(eigen(ch$ell_A, only.values=TRUE)$values[1]), 3)  )}) |> print(n = 25)

# A tibble: 25 × 9   Chart `n pts` `centre (x,y,z)`    `M · r̂`     k1     k2     k3 `ell axis a`   <int>   <int> <chr>                 <dbl>  <dbl>  <dbl>  <dbl>        <dbl> 1     1      82 0.91, 0.3, -0.16      0.999 -0.216 -0.078 -0.068        0.498 2     2      70 0.35, -0.9, -0.06     1     -0.136  0     -0.163        0.486 3     3      84 -0.64, 0.71, -0.13    1      0.182 -0.031  0.101        0.44  4     4      73 -0.04, 0.95, -0.12    1      0.162 -0.024  0.112        0.463 5     5      81 -0.75, 0.4, 0.46      1      0.115 -0.056  0.184        0.431 6     6      73 0.28, 0.6, -0.7       0.999 -0.203  0.016 -0.114        0.518 7     7      78 -0.93, 0.07, -0.22    1      0.169 -0.087  0.055        0.46  8     8      86 -0.08, -0.3, 0.91     1     -0.237  0.051  0.099        0.53  9     9      73 -0.18, -0.87, -0.38   1      0.111 -0.036  0.127        0.42110    10      85 0.87, 0.07, 0.4       1     -0.148 -0.002 -0.13         0.47111    11      75 -0.36, 0.21, 0.87     1      0.2   -0.034  0.054        0.55612    12      67 -0.79, -0.52, -0.21   1      0.178 -0.001  0.109        0.43113    13      84 0.04, 0.05, -0.96     1      0.182  0.012  0.081        0.44614    14      66 -0.39, -0.83, 0.31    1     -0.189 -0.006 -0.066        0.51115    15      84 0.28, -0.55, -0.74    1      0.132 -0.066  0.105        0.43116    16      87 -0.52, 0.43, -0.68    1     -0.251  0.009  0.026        0.51117    17     100 0.83, -0.44, -0.2     1      0.144  0.082  0.133        0.52618    18      68 0.07, -0.8, 0.55      1     -0.24   0.072 -0.088        0.47419    19      91 0.59, 0.71, 0.24      0.999  0.229 -0.002 -0.003        0.53720    20      93 -0.81, -0.26, 0.45    1      0.17   0.084  0.115        0.45521    21      95 0.41, 0.27, 0.83      1      0.155 -0.087  0.056        0.53222    22      66 -0.46, -0.32, -0.78   1     -0.179 -0.025 -0.072        0.52523    23      70 -0.12, 0.78, 0.54     1      0.183  0.045  0.1          0.49524    24      93 0.67, 0.07, -0.69     1      0.192 -0.001  0.116        0.45325    25      76 0.62, -0.45, 0.58     1     -0.215 -0.039  0.028        0.512# ℹ 1 more variable: `ell axis b` <dbl>

Key things to notice:

M · r̂ is uniformly $\approx 1$ , confirming the SVD-learned normal aligns with the sphere’s outward radial direction.
k1, k2, k3 encode local curvature. The Hessian $\begin{pmatrix}2k_1 & k_2 \\ k_2 & 2k_3\end{pmatrix}$ should have eigenvalues $\approx \pm 1$ for a unit sphere (developed further in the curvature proof section).
axes a/b semi-axes being nearly equal (near-circular domains) reflects the sphere’s isotropy.

The following plot shows how the charts cover the sphere. Each point is colored by the chart it belongs to. With $k = 25$ charts, the sphere is tiled into roughly equal-area patches.

Show the code

#| label: atlas-chart-plotn_charts <- length(unique(sphere_pts$chart))chart_colors <- colorRampPalette(RColorBrewer::brewer.pal(9, "Set1"))(n_charts)plot_ly(  data   = sphere_pts,  x = ~x, y = ~y, z = ~z,  type   = "scatter3d",  mode   = "markers",  color  = ~chart,  colors = chart_colors,  marker = list(size = 2.5, opacity = 0.6),  text   = ~paste0("Chart ", chart),  hoverinfo = "text") |>  layout(    title  = "S² coloured by Atlas-Learn chart membership (k = 25)",    scene  = list(      xaxis = list(title = "x"),      yaxis = list(title = "y"),      zaxis = list(title = "z"),      aspectmode = "cube"    ),    showlegend = FALSE  )

To leave a comment for the author, please follow the link and comment on their blog: R Works.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Continue reading: The Atlas-Learn Approach to the Manifold Hypothesis