Embedding Fermat distances on a paraboloid, saddle and helix • grip

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.

This example fits three-dimensional metric MDS directly to dense complete-graph Fermat distance matrices. No kNN graph is required. We then deliberately restrict paths to symmetric Fermat kNN graphs with k from 3 through 10 and compare their embeddings with the dense-target fits. The companion dgraphs vignette explains the distance estimator and supporting path graphs.

Sample surfaces and supply a distance matrix

We use the paraboloid $(u,v,u^2+v^2)$ and saddle $(u,v,u^2-v^2)$ over the unit disk. The existing dgraphs sampler draws uniformly with respect to surface area, using the area element of the graph surface; we do not replace it with uniform sampling of the two parameter coordinates.

library(grip)
s <- dgraphs::sample.synthetic.geometry(
  dgraphs::synthetic.quadform(2,3,forms=list(diag(c(1,-1)))),
  dgraphs::synthetic.sampling.quadform.lab('disk',mode='area'),
  n=150,seed=4101)
X <- s$predictors
D <- dgraphs::fermat.distances(points=X,p=2,backend='implicit',rooted=FALSE)
fit <- metric.mds(distance.matrix=D,dim=3,backend='sgd',init='random',
  seed=3051,max.iter=300,diagnostics=FALSE,
  sgd.control=list(final.rate=1e-5))
Y <- fit$coords

For observations $X_i$ , the supplied target is $D_p(i,j)=\min_{\pi:i\leadsto j}\sum_{(a,b)\in\pi}\|X_a-X_b\|^p$ . The minimum allows every observation as an intermediate point. The exponent is applied to individual ambient steps, and the returned distance is unrooted. The MDS coordinates $Y_i\in\mathbb R^3$ seek to minimize an unweighted sum of squared residuals $(\|Y_i-Y_j\|-D_p(i,j))^2$ over distinct pairs. The fixed 300-epoch SGD budget is an optimization budget, not a convergence certificate.

A matched comparison

The comparison uses powers 1, 2, 3 and 4. Every fit uses the same seeded random starting configuration for retained samples, with the initial global scale fitted separately to its target matrix, uses uniform pair weights and the same optimization controls. For a disconnected kNN graph only its largest component is embedded. All within-component graph shortest-path distances are supplied as a dense matrix; this isolates graph restriction, without adding a landmark approximation to this small example.

We keep three questions separate. Distance error compares the input distances with reference surface geodesics. Graph restriction error compares graph shortest paths with complete-graph Fermat distances. Embedding error compares Euclidean distances in the fitted coordinates with their supplied targets. For vectors of estimated distances $d$ and reference distances $g$ , we report $E(d,g;a)=\sqrt{\frac{\sum_r(ad_r-g_r)^2}{\sum_r g_r^2}}.$ Here $r$ indexes valid comparison pairs. For distance-versus-geodesic comparisons, $a=\sum_r d_rg_r/\sum_r d_r^2$ fits one global scale because different powers have different units. For embedding and graph restriction error, $a=1$ : the targets and their estimates are in the same units. Smaller errors are better; zero means agreement on those pairs. Reference geodesics use existing continuous solvers on 1,000 sampled pairs using best_of_six for both surfaces, with paraboloid_clairaut on the same paraboloid pairs. Each of the six refinement searches receives 64 rounds; no global-optimality certificate is implied.

Table 1. Dense-target fits. The two geodesic columns allow a global scale change; the embedding-versus-input column does not. A good MDS fit can faithfully reproduce an inaccurate geodesic estimator, and these columns distinguish those cases.

Figure 1. Scaled agreement between embedded Euclidean distances and the intrinsic reference. The point at k=0 is the complete-target baseline; points at k=3–10 are restricted graphs. Connecting segments help follow powers and do not represent fits at untested k values. Disconnected results use fewer pairs:

The table summarizes ranges across disconnected k values. Complete coverage is 150 samples. Small within-component errors do not establish whole-sample recovery.

Overlay embeddings and their supporting paths

The reference surface is translucent gray, original observations are gray points, and fitted coordinates are colored by their first latent coordinate. Thin blue lines show the selected shortest-path edges for the complete graph, or Fermat neighbor edges for the restricted graph. Views are aligned to the observations using least-squares similarity Procrustes alignment: translation, rotation or reflection, and one uniform scale, matching each embedded point to its original observation. Plot titles report the root mean squared Euclidean alignment residual per retained sample, in target-coordinate units (zero is exact overlap). Alignment changes only display coordinates; numerical errors above use the original fits.

Figure 2. Dense p=2 Fermat embedding on the paraboloid. Rotate the interactive view to inspect displacements from the gray surface.

Figure 3. Dense p=2 Fermat embedding on the saddle, with its supporting paths. An intrinsic metric need not have an exact three-dimensional Euclidean realization.

Figure 4. The same saddle samples after restriction to a symmetric k=3 Fermat graph. Its edges summarize potentially multistep Fermat paths; they are not the individual powered ambient edges shown in Figure 3. If disconnected, only the largest component appears. Edge legends report any display thinning beyond 1,800 edges; underlying graph data remain complete.

Why the helix should unfold

The helix $(\cos t,\sin t,0.2t)$ has arc-length distance $\sqrt{1.04}|t_i-t_j|$ . Therefore its intrinsic distance matrix has an exact one-dimensional realization. Fitting those distances in three coordinates should produce a nearly straight line, not reconstruct the coil. Three-dimensional metric MDS compares straight distances, not travel constrained to a curve.

Figure 5. Dense p=2 Fermat embedding over the original helix. Agreement with the coil is not the success criterion for intrinsic-distance reconstruction.

Figure 6. Exact arc-length target baseline. Its coordinates should approach a line; remaining transverse spread measures finite optimization error. The surface examples do not have corresponding all-pairs exact-reference baselines: their numerical reference checks cover selected pairs only.

Larger datasets

The direct dense API above is suitable for small comparisons. For larger datasets, dgraphs can compute only selected source rows with sources = landmarks. landmark.mds.constraints() converts those distances into explicit weighted pair fitting targets for sparse metric.mds(). This is a separate approximation in the embedding objective, not a restriction of the paths defining each Fermat distance.

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.