Fermat distances: sampling density, curvature, and geodesic recovery

What distance is being estimated?

For observations $X_1,\ldots,X_n\in\mathbb R^3$, the unrooted Fermat distance is

$$D_p(i,j)=\min_{\pi:i\leadsto j} \sum_{(u,v)\in\pi}\lVert X_u-X_v\rVert_2^p,\qquad p\geq1.$$

The minimum ranges over paths through any observed points, including the direct edge. At $p=1$, the triangle inequality gives $D_1(i,j)=\lVert X_i-X_j\rVert_2$: adding intermediate observations cannot shorten the ambient chord. For $p>1$, several short steps can cost less than a long step. More intermediate points can therefore change both the preferred route and its cost. The units are length to power $p$; this article takes no $p$-th root.

A surface geodesic measures length along the underlying surface. We compare Fermat distances to that intrinsic reference after fitting one multiplicative scale. These are finite-sample experiments, not a proof of asymptotic convergence or an estimate of its rate.

A compact example using the existing area sampler

Consider the saddle $(u,v,a(u^2-v^2))$ over $u^2+v^2\leq1$. For a quadratic form matrix $A$, the surface-area element is $\sqrt{1+\lVert2A(u,v)^T\rVert_2^2}$. Uniform draws in the parameter disk are therefore not uniform on the surface. Use the existing dgraphs sampler:

library(dgraphs)
A <- diag(c(2,-2))
geometry <- synthetic.quadform(2L,3L,forms=list(A))
sampling <- synthetic.sampling.quadform.lab('disk',extent=1,mode='area')
sample <- sample.synthetic.geometry(geometry,sampling,n=300L,seed=62117L)
X <- sample$predictors
endpoints <- 1:12
D150 <- fermat.distances(points=X[1:150,],p=1.25,backend='implicit',
                        sources=endpoints,targets=endpoints,rooted=FALSE)
D300 <- fermat.distances(points=X,p=1.25,backend='implicit',
                        sources=endpoints,targets=endpoints,rooted=FALSE)
stopifnot(all(D300 <= D150 + 1e-10))

Only the requested endpoint distances are returned, but paths may use every observation in the supplied matrix. This is an exact restriction of the queries, not a landmark restriction of the paths. The check illustrates that unrooted Fermat distances cannot increase when the available path vertices expand. These small teaching observations are separate from the frozen benchmark below.

For the ordinary neighborhood baseline, construct a union-symmetric graph: an edge exists if either endpoint selects the other among its $k$ nearest neighbors. Its weights remain Euclidean lengths, without a Fermat power.

graph <- create.sknn.graph(X,k=10,neighbor.method='exact',
                          connect.components=FALSE)
edges <- graph.edges(graph)
ig <- igraph::graph_from_data_frame(edges[,1:2],directed=FALSE,
                                  vertices=data.frame(name=seq_len(nrow(X))))
Dknn <- igraph::distances(ig,v=endpoints,to=endpoints,weights=edges$length)
igraph::components(ig)$no

If the graph is disconnected, some entries can be infinite. This study records that limitation rather than adding repair edges. This is an ordinary Euclidean kNN graph, distinct from create.fermat.sknn.graph(), which selects neighbors using Fermat distances.

The seven-surface benchmark

The quadratic surfaces are

$$S_a^+(u,v)=(u,v,a(u^2+v^2)),\qquad S_a^-(u,v)=(u,v,a(u^2-v^2)),\qquad u^2+v^2\leq1,$$

with $a\in\{1,2,4\}$. Principal curvature magnitudes at the origin are $2a$. Increasing $a$ also changes the surface area and curvature away from the origin; it does not uniformly multiply curvature everywhere. Both patches have area

$$A(a)=\frac{\pi}{6a^2}\left[(1+4a^2)^{3/2}-1\right].$$

Figure 1. The six quadratic patches. A common vertical range within each family makes the increased bending and height visible.

The seventh surface is $(t\cos t,h,t\sin t)$, where $1.5\pi\leq t\leq4.5\pi$ and $0\leq h\leq20$. Its arc coordinate is

$$s(t)=\frac12\left[t\sqrt{1+t^2}+\operatorname{asinh}(t)\right].$$

Uniform draws in the rectangle with coordinates $(s,h)$ give area-uniform samples. Its exact intrinsic distance is $g_{ij}=\sqrt{(s_i-s_j)^2+(h_i-h_j)^2}$.

#> The installed guide uses the static surface views in Figure 1. Interactive views are available when rebuilding with DGRAPHS_STATIC_VIGNETTES=false and the optional visualization packages.

Compare distance structure on a separate pair set

For calibration pairs $C$, evaluation pairs $T$, estimates $D_q$, and reference lengths $g_q$, fit the scale and evaluate the normalized error:

$$\widehat b=\frac{\sum_{q\in C}D_qg_q}{\sum_{q\in C}D_q^2},\qquad E_T=\sqrt{\frac{\sum_{q\in T}(\widehat bD_q-g_q)^2} {\sum_{q\in T}g_q^2}}.$$

Zero means proportional agreement with the reference. Plots show $100E_T$ as a percentage. Raw errors with $\widehat b=1$ are saved too; they answer a different question because powered lengths have different units.

Unreachable graph pairs are excluded from their finite-distance score, and coverage is reported. Scale fitting likewise uses only finite calibration pairs. Paired differences use ambient scores on the same reachable evaluation pairs. A small error with incomplete coverage does not demonstrate recovery of the entire surface.

Table 1. Coverage over all surfaces, repetitions and sample sizes. “Adaptive” denotes the fixed logarithmic rule above. A disconnected component may contain no evaluated endpoints, so connectivity and pair coverage are reported separately.

Intermediate powers and denser samples

Figure 2. Paraboloid errors against best-of-six references. Thin lines are individual repetitions; bold lines are means. Shading gives pointwise 95% Student-$t$ intervals with four degrees of freedom. The horizontal axis is logarithmic. All intervals here and below are conditional on the fixed endpoints and references, not simultaneous over the tested grid, and omit solver uncertainty.

At $a=1$, $p=1.25$ already improves on ambient chords at $n=150$: mean error is 3.73% versus 8.51%. For $p=2$, error falls from 11.72% to 5.67% as the sample increases to 2,400. Thus both the choice of power and the availability of intermediate observations matter. On these three paraboloids, $p=1.25$ is the best tested complete power throughout the sample-size grid.

Figure 3. The same comparison on saddles. On the strongest saddle at $n=2400$, $p=1.5$ has 2.55% error versus 28.32% for ambient and 11.99% for $p=1.25$. A power useful on a mild patch need not be sufficient on a more strongly bent patch.

Figure 4. Swiss-roll errors against exact intrinsic lengths. Ambient error is 55.84%; at $n=2400$, $p=3$ achieves 7.70%, whereas $p=1.25$ still has 51.99%. The best tested power changes from $p=4$ at $n=150$ to $p=3$ at $n=2400$.

Table 2. Descriptive minima over the tested powers at $n=2400$. Errors are means over five repetitions; wins count paired repetitions with lower error than ambient. These powers were identified using the displayed evaluation scores, so this table is not an independent assessment of a tuned selection procedure. A separate calibration-only power search at $n=300$ and transfer to the other sample sizes is reported below, after the original fixed-grid comparisons.

How much do ambient chords underestimate surface lengths?

Figure 5. Distributions of $R_q=g_q/c_q$ on the 1,000 evaluation pairs per surface, where $c_q$ is the ambient chord. One denotes equality. The axis has $\log_2$ spacing but is labeled in ratio units; two means twice the chord length. Violin widths show density on the log-ratio scale, with medians and quartiles inside. Shared endpoints make these descriptive pair distributions.

The rolled surface exposes very large shortcuts. Quadform coefficients change the distribution too, but need not increase every pair’s distortion: curvature varies spatially, and area-uniform sampling puts more mass in steeper regions.

Ordinary kNN distances are a strong baseline

Figure 6. Ordinary Euclidean kNN errors on paraboloids (top), saddles (middle) and the Swiss roll (bottom). Thin lines are repetitions; bold means and error bars use conditional 95% $t$ intervals. Intervals can extend below zero, although errors cannot. The dashed ambient baseline uses all evaluation pairs; consult Table 1 before comparing it with a restricted finite-pair score.

At $n=2400$, $k=20$ is more accurate than every tested complete Fermat power on each of the seven examples. Its errors on the mild paraboloid, mild saddle and Swiss roll are 0.62%, 0.49% and 0.35%. Larger neighborhoods are not automatically better: on the roll, $k=40$ has 13.03% error, consistent with shortcut edges. Connectivity alone is not a test of geodesic accuracy.

These findings concern ordinary geodesic recovery under uniform sampling. They do not settle whether density-sensitive distances are useful on nonuniform biological data, where the desired geometry itself may depend on density.

Reference accuracy can change the apparent error floor

Figure 7. The paraboloid comparison using Clairaut lengths. Observations, pairs and estimated distances are identical to Figure 2; only the reference and its fitted calibration scale change.

For $a=4$, the $p=1.25$ error at $n=2400$ changes from 4.82% against best-of-six to 2.05% against Clairaut. The $k=20$ error changes from 4.24% to 0.60%. Some apparent estimator error is therefore reference-method disagreement. The original best-of-six results remain part of the evidence.

A separate saddle diagnostic increased search rounds from 64 to 256 on 50 regularly spaced pair indices per saddle, selected without inspecting errors. The largest relative length reduction over those 150 pairs was 0.137%. This supports stability for those pairs and budgets, but does not certify global minimality or remove reference uncertainty elsewhere.

Paired improvements over ambient distances

Figure 8. Scaled error minus ambient error at $n=2400$, in percentage points. Negative values favor Fermat. Black points are individual repetitions; blue means and intervals show conditional 95% $t$ intervals. Paired comparisons hold the endpoint set and added sample pool fixed within each repetition.

Selecting a power at 300 observations and transferring it

This separate follow-up tunes each of the six quadforms at $n=300$, then keeps that power fixed at $n\in\{150,300,600,1200,2400\}$. It reuses the original area-uniform samples, fixed endpoints, five repetitions, and numerical reference lengths. No surface samples or reference geodesics were recomputed. The Swiss roll is outside this tuning experiment; its original results above are unchanged.

For each surface, the 500 calibration pairs were randomly assigned to five fixed folds of 100 pairs. The same folds were used for every power, repetition, and reference method. For repetition $r$, power $p$, and validation fold $F_j$, the scale was fitted on the other 400 pairs:

$$\widehat c_{r,p,-j}= \frac{\sum_{q\in\mathcal C\setminus F_j}D_{r,p,q}g_q} {\sum_{q\in\mathcal C\setminus F_j}D_{r,p,q}^{2}}, \qquad E_{r,p}^{\mathrm{CV}}= \sqrt{\frac{\sum_{j=1}^{5}\sum_{q\in F_j} (\widehat c_{r,p,-j}D_{r,p,q}-g_q)^2} {\sum_{q\in\mathcal C}g_q^2}}.$$

Here $\mathcal C$ is the calibration-pair set, $D_{r,p,q}$ the unrooted complete-graph Fermat distance, and $g_q$ the reference geodesic length. Selection minimizes $\overline E_p^{\mathrm{CV}}=\frac15\sum_{r=1}^{5} E_{r,p}^{\mathrm{CV}}$. Pooling the out-of-fold squared residuals before taking the square root gives each pair its contribution to the same normalized error; we do not average five differently normalized fold scores.

The coarse grid has spacing $0.025$ on $[1,1.5]$ for all paraboloids and the coefficient-one saddle, and on $[1,2]$ for the two stronger saddles. The best coarse power is refined over its $\pm0.025$ neighborhood in steps of $0.005$, clipped to the search interval. The reported choice is the best power found at this resolution, not a continuous optimum. The ambient baseline $p=1$ participated in selection. Exact ties would favor the smaller power; none of the nine selected powers lies on a search boundary.

All nine choices were saved before evaluation: six using best-of-six references and three independent paraboloid selections using Clairaut references. The follow-up completed 301 surface/reference/power combinations, each with five repetitions and five folds (7,525 fold scores), and 1,275 transfer comparison records. Every transfer record uses all 500 calibration and 1,000 evaluation pairs. No distance condition failed or was excluded. The number of records includes baseline comparisons and scoring each paraboloid selection against both reference methods; it is not a count of distinct distance searches.

Figure 9. Full tuning curves at $n=300$ using best-of-six references. Thin lines show the five repetitions; blue curves show their mean and pale bands pointwise conditional 95% $t$ intervals. Red diamonds and dotted lines mark the selected power. Panels use different horizontal ranges because the two stronger saddles were searched through $p=2$; vertical ranges also differ. All coarse and refined grid points are included. Smaller error is better.

Figure 10. Neighborhoods within $0.030$ of the selected power. The horizontal coordinate is $p-p_{\mathrm{chosen}}$; the vertical coordinate is $100(\overline E_p^{\mathrm{CV}}/\overline E_{p_{\mathrm{chosen}}}^{\mathrm{CV}}-1)$, the percentage increase relative to the minimum mean error. The dashed line marks a 1% relative increase. Bands use paired differences between each candidate and the selected power within each repetition, with pointwise conditional 95% $t$ intervals divided by the observed minimum mean error. They are not confidence intervals for the optimum’s location. Dots are evaluated powers and lines only connect them.

Table 3. Best powers using best-of-six references, a descriptive near-minimum range, and mean held-out errors after refitting scale on all calibration pairs. The range contains the smallest and largest tested powers with mean cross-validation error no more than 1% above its minimum; it is not a confidence interval and does not classify untested powers. The 1% threshold is relative to the minimum error, not one percentage point of distance error. The minima span two to eight qualifying grid points. Their breadth varies with the surface, so the third decimal place records search resolution rather than certainty about a uniquely preferred power.

Figure 11. Transfer of the power selected at $n=300$ to the five sample sizes, scored against best-of-six references. Scale is refitted on the 500 calibration pairs at each size and repetition; the 1,000 evaluation pairs are unchanged. Ambient chords, fixed powers $1.25$ and $2$, and ordinary Euclidean symmetric-kNN shortest paths with $k=20$ provide comparisons. Thin lines show repetitions; thicker lines and points show means, with pointwise conditional 95% $t$ intervals. The sample-size axis is logarithmic and panel vertical ranges differ. All displayed neighborhood graphs are connected.

The selected powers improve on ambient chords on all six surfaces at every sample size, in all five repetitions. Their mean errors decrease as the sample size increases. At $n=2400$, errors range from 1.00% on the mild saddle to 4.78% on the steepest paraboloid against best-of-six references. This supports useful transfer, but does not mean that selection at $n=300$ remains optimal at every size.

Figure 12. Selected-power error minus fixed-$p=1.25$ error on the same evaluation pairs and sample pools, in percentage points. Negative values favor the selected power. Thin lines show the five paired differences; blue means and pointwise conditional 95% $t$ intervals summarize them.

For example, on the intermediate saddle the selected $p=1.335$ has 3.95% error at $n=300$, versus 5.91% for $p=1.25$. At $n=2400$ the ordering reverses: 2.00% versus 1.58%. On the strongest saddle, transferred $p=1.680$ reaches 3.25%, whereas the original grid’s $p=1.5$ reaches 2.55%. Selection at one sample size therefore identifies a useful finite-sample choice, not a universal power. Ordinary Euclidean $k=20$ graph distances remain more accurate at $n=2400$ than these selected complete-graph Fermat distances on all six examples.

Reproduce the evidence without slowing article builds

The installed snapshot can be inspected directly:

path <- system.file('extdata','fermat-density-curvature','benchmark.rds',
                    package='dgraphs',mustWork=TRUE)
benchmark <- readRDS(path)
nrow(benchmark$metrics)       # 2,100 condition-level records
nrow(benchmark$references)    # 10,500 main reference pairs
benchmark$power_tuning$selection
nrow(benchmark$power_tuning$transfer_metrics)  # 1,275 follow-up records
benchmark$provenance

The companion surface-views.rds contains only the seven saved 300-point coordinate sets used in the interactive display; surface-views-provenance.R records their source checksums and sampling repetition.

The benchmark snapshot contains condition-level errors and coverage, summaries, reference comparisons, the saddle sensitivity diagnostic, the power-selection fold scores and transfer results, and provenance. All data are synthetic. The companion provenance.R records source-file checksums and the research revision from which the snapshot and figures were imported. The full research experiment is named fermat_density_curvature_01_oct_2026.

To reproduce the expensive calculations, copy the installed scripts to a writable directory:

src <- system.file('examples','fermat-density-curvature',package='dgraphs',
                   mustWork=TRUE)
work <- file.path(tempdir(),'fermat-recovery-reproduction')
dir.create(work,showWarnings=FALSE)
file.copy(list.files(src,full.names=TRUE),work,overwrite=FALSE)

From that directory run the following commands. This is intentionally separate from rendering the article. It requires dgraphs, igraph and jsonlite; use one worker on platforms without forked multiprocessing.

Rscript --vanilla prepare.R
Rscript --vanilla references.R 1
Rscript --vanilla distances.R 1
Rscript --vanilla reference_sensitivity.R
Rscript --vanilla evaluate.R
Rscript --vanilla validate.R
Rscript --vanilla power_tuning.R 1
Rscript --vanilla power_transfer.R 1
Rscript --vanilla power_validate.R

Reference searches can be expensive. Completed chunks are reused; use a new working directory when changing the design. The scripts use an internal numerical reference interface, so retain the recorded package versions and source revision when interpreting differences between reruns. The packaged snapshot remains the fixed evidence for this article.

The follow-up tuning script refuses to overwrite a frozen selection. Its recorded checksum is checked again during evaluation. The additional power_summarize.R script regenerates follow-up figures using ggplot2; article builds read saved figures and do not repeat power selection.

The practical conclusion is to compare useful intermediate powers, sample sizes, and an ordinary neighborhood-graph baseline before attributing a result to Fermat geometry. Keep distance-estimation error separate from embedding distortion, and show reference uncertainty when it can change the comparison.