Fermat distances and their supporting graphs on surfaces • dgraphs

This executable example requires development versions of dgraphs and grip with the direct distance and supporting-graph APIs, plus igraph and plotly. Install these dependencies to render the calculated examples.

Fermat distances allow paths through every observation. A symmetric nearest-neighbor graph is an optional restriction, not part of the complete-graph estimator. This example asks how that restriction changes distances on an area-uniform paraboloid, an area-uniform saddle, and an arc-length-uniform helix. The companion grip vignette focuses on embedding the same distances into three dimensions.

Generate the observations

The quadratic surfaces are $(u,v,u^2+v^2)$ (index zero) and $(u,v,u^2-v^2)$ (index one), over the unit disk. Index counts negative eigenvalues of the quadratic form. We reuse the existing dgraphs area sampler:

library(dgraphs)
geometry <- synthetic.quadform(2, 3, forms=list(diag(c(1,-1))))
sampling <- synthetic.sampling.quadform.lab('disk',extent=1,mode='area')
sample <- sample.synthetic.geometry(geometry,sampling,n=150,seed=4101)
X <- sample$predictors

The sampler weights parameter-space proposals by $\sqrt{1+\|2Au\|^2}$ , the area element for $(u,u^TAu)$ , where $u$ is a two-dimensional parameter and $A$ is its quadratic-form matrix. This is uniform on surface area, not on the parameter disk. The helix is $(\cos t,\sin t,0.2t)$ for $0\leq t\leq4\pi$ ; uniform $t$ is uniform arc length because its speed is the constant $\sqrt{1+0.2^2}$ .

Distances and the path graph

For observation indices $i,j$ and exponent $p\geq1$ , we compute $D_p(i,j)=\min_{\pi:i\leadsto j}\sum_{(a,b)\in\pi}\|X_a-X_b\|^p.$ The minimum ranges over paths through the observed points, including direct edges. The distance is unrooted, so its units are length to power $p$ .

result <- fermat.distances(points=X,p=2,backend='implicit',rooted=FALSE,
                          return.graph=TRUE)
D <- result$distances
graph <- result$graph
head(graph.edges(graph))
result$metadata

The returned graph is the union of one selected shortest-path tree per source. Its edges carry the powered ambient step lengths. With every source requested, its shortest-path distances preserve the full matrix. With only landmark sources, only landmark-to-sample preservation is guaranteed. Ties may admit other valid unions, and the union need not be sparse: at $p=1$ , direct Euclidean edges usually remain necessary.

A Fermat symmetric kNN graph has different edges: they summarize Fermat paths to the $k$ nearest neighbors, then take the union of directed neighbor selections.

g <- create.fermat.sknn.graph(X,k=5,p=2,rooted=FALSE)

Compare with surface geodesics

The shared example uses powers 1–4 and k values 3–10. Numerical references on 1,000 reproducibly selected pairs use best_of_six in the existing quadform_geodesics() interface on both surfaces. The paraboloid also uses paraboloid_clairaut on the same 1,000 pairs. Best-of-six selects the shortest usable result from six refinement searches, each with 64 rounds. These are numerical references, not certified globally shortest paths. Helix references are exact arc lengths $\sqrt{1.04}|t_i-t_j|$ .

We compare estimates $d_r$ with reference distances $g_r$ on valid pairs $r$ : $a=\frac{\sum_r d_rg_r}{\sum_r d_r^2},\qquad E=\sqrt{\frac{\sum_r(ad_r-g_r)^2}{\sum_r g_r^2}}.$ The fitted global scale $a$ makes powers comparable; $E=0$ means proportional agreement. It does not establish equality in physical units. Spearman correlation measures rank agreement. For graph restriction error we compare graph distances with $D_p$ using the same formula with $a=1$ , on all pairs in the largest component. Disconnected graphs are not repaired, and retained sample counts are reported. If a retained component contains no reference pair, its geodesic error is undefined and omitted from the corresponding curve.

Table 1. Complete-graph Fermat distances compared with the surface or curve reference. Each row uses the same reference pairs within its example; these are finite-sample comparisons, not estimates of an asymptotic convergence rate.

Figure 1. Distance change caused by retaining only symmetric Fermat kNN edges. Zero means that the restricted graph preserves the dense Fermat distances on the retained component. The corresponding component counts and retained sample counts are available below; a small error on a small component would not demonstrate whole-sample preservation.

The table summarizes ranges among disconnected k values within each example and power. Complete coverage is 150 samples; restricting a disconnected graph can leave far fewer samples even when its within-component distance error is small.

Inspect supporting edges in three dimensions

The following view uses dense-target SGD metric MDS from grip. The gray patch is the underlying surface; colored points are the embedding after one global least-squares similarity Procrustes alignment for display, and thin blue segments are supporting path edges. Gray points identify the original observations. Plot titles report the root mean squared Euclidean alignment residual per retained sample, in target-coordinate units; zero means exact overlap. Similarity alignment uses translation, rotation/reflection and uniform scale, not a deformation of the fit.

Figure 2. Saddle embedding and complete-graph Fermat supporting paths. Drag to rotate. The legend states how many edges are shown; a deterministic display subset is used above 1,800 edges. The saved graph contains every selected edge.

Intrinsic distances measure travel along a surface. Euclidean distances in a three-dimensional embedding measure straight chords, so geometric recovery of the original surface is not implied by a good distance fit. In particular, exact helix arc lengths embed on a straight line. See the grip companion for the full embedding comparison and additional overlays.

Chord lengths before and after scaling

At p=1, complete-graph Fermat distances equal ambient Euclidean chord lengths. A 3D MDS fit can therefore reproduce the sampled shape without reproducing its intrinsic distances. The earlier distance score fits a global scale; it is not the discrepancy of the original chord lengths.

For each of the same evaluated pairs, write $c_r$ for the chord and $g_r$ for the intrinsic reference. The following table reports $E(c,g;1)$ and $E(c,g;a)$ with $a=\sum c_rg_r/\sum c_r^2$ , together with the median and maximum ratio $g_r/c_r$ . Errors are dimensionless and zero means agreement. The ratio is one for equal lengths; 1.3 means the geodesic is 30% longer than its chord. These are the 1,000 evaluated pairs for seed 4101, not all sample pairs.

Scaling reduces chord underestimation on these patches, but does not make chords intrinsic distances. Larger discrepancies remain for individual endpoint pairs. On the helix, accurately preserving arc lengths unfolds the curve into a line; matching the original coil is a different objective.

Independent solver checks are attempts, not guaranteed successes. Candidate lengths and close agreement where both solvers succeed do not certify global minimality. The table below makes the successful comparison count explicit.

Paraboloid reference-method comparison

Best-of-six and Clairaut evaluate the same 1,000 pairs. Positive signed relative difference $(b_r-h_r)/h_r$ means the best-of-six candidate $b_r$ is longer than the Clairaut candidate $h_r$ . Neither is assumed globally certified.

The second table compares scaled Fermat-distance errors under each reference. It uses identical targets and endpoint pairs, isolating reference-method effects.

Follow-up: sampling density, intermediate powers and stronger curvature

A separate distance-only experiment varied total sample size from 150 to 2,400, with five repetitions, on paraboloids and saddles with coefficients 1, 2 and 4, and an area-uniform Swiss roll. It reused the existing quadform area sampler. For each surface, 64 fixed sampled endpoints supplied 500 calibration pairs and 1,000 different evaluation pairs. Additional path-support points varied across repetitions and were nested across sample sizes. All observations could be used by the complete-graph paths; this was not a landmark approximation.

Errors below are percentages: the square root of the ratio of summed squared distance residuals to summed squared reference distances. Each method’s single multiplicative scale was fitted on calibration pairs only. Values are means over five repetitions conditional on the fixed endpoints; they are new measurements, not replacements for the original 150-point fits above. No MDS error enters them.

Surface	Power p	Error at n=150 (%)	Error at n=2,400 (%)	Ambient error (%)
Paraboloid, coefficient 1	1.25	3.73	1.82	8.51
Paraboloid, coefficient 1	2	11.72	5.67	8.51
Saddle, coefficient 1	1.25	4.55	1.75	8.74
Saddle, coefficient 1	2	14.27	5.58	8.74
Saddle, coefficient 4	1.5	14.82	2.55	28.32
Swiss roll	3	42.33	7.70	55.84

Ambient error is constant because the endpoints are fixed. These results expose both a sample-size effect and a limitation of the original integer-only power grid: p=1.25 already beats ambient chords on the new mild-surface samples at n=150. Larger powers are needed to suppress shortcuts on the rolled surface. The displayed powers illustrate these effects; they are not universal tuning recommendations.

Ordinary Euclidean kNN graph distances were more accurate on these uniform examples at the larger sample size. For k=20 and n=2,400, errors were 0.62%, 0.49% and 0.35% on the mild paraboloid, mild saddle and Swiss roll. This separates recovering ordinary geodesics from deliberately imposing a density-weighted geometry.

References used best-of-six on quadforms and analytic distances on the roll, with Clairaut sensitivity checks on every paraboloid pair. For the steepest paraboloid, p=1.25 at n=2,400 had 4.82% error against best-of-six but 2.05% against Clairaut. Reference-method disagreement therefore matters for this harder surface. Increasing the quadform coefficient changes its curvature profile and area; it does not multiply curvature uniformly over the patch. The complete experiment, pair definitions, coverage and uncertainty are recorded in the research experiment fermat_density_curvature_01_oct_2026.