Compute Complete or Graph-Restricted Sample Fermat Distances

Minimize sums of powered edge lengths. Coordinate input defaults to the complete Euclidean graph, implementing the sample Fermat distance on the observed points. A symmetric kNN restriction or a supplied graph is an explicit alternative, with no automatic guarantee of the same distances.

Usage

fermat.distances(
  graph = NULL,
  p = 2,
  rooted = FALSE,
  vertices = NULL,
  stage = "final",
  points = NULL,
  graph.type = NULL,
  k = NULL,
  backend = "explicit",
  sources = NULL,
  targets = NULL,
  max.workspace.bytes = 512 * 1024^2,
  return.graph = FALSE
)

Arguments

graph: A dgraph with finite nonnegative edge lengths, interpreted as base lengths before powering. Supply exactly one of graph and points.
p: Finite numeric scalar at least one; the edge-length exponent.
rooted: Logical scalar. If TRUE, return the 1/p power of the minimum sum. Defaults to FALSE; rooting does not change minimizing paths.
vertices: Optional unique 1-based vertex indices in output order. Paths may use all vertices. NULL selects all; an empty vector selects none.
stage: Stored graph stage, default "final". Coordinate input only permits "final".
points: Finite numeric matrix or data frame, observations in rows.
graph.type: Coordinate construction: "complete" (default) or "sknn" (union-symmetric Euclidean kNN, without pruning or repair). Cannot accompany graph.
k: Required for graph.type = "sknn", with 1 <= k < n. Not accepted for complete or supplied graphs.
backend: Complete-graph coordinate backend: "explicit" (the default reference) or "implicit" (evaluate edges on demand without storing adjacency). Explicitly supplied only for complete coordinate input.
sources: Optional unique source indices for rectangular queries with backend = "implicit". For pivot-to-all distances, supply the pivot ids. Cannot accompany vertices; every point remains a possible transit vertex.
targets: Optional unique target indices with sources; defaults to all observations. Row and column order follow sources and targets.
max.workspace.bytes: Positive finite workspace allowance for the implicit backend (default 512 MiB). The estimate includes coordinates, output and linear native scratch, not total R process memory or R copies.
return.graph: Logical scalar; return distances together with the union of selected shortest-path edges and coverage metadata. Defaults to FALSE. Complete coordinate input uses the implicit implementation for this option, even when backend = "explicit" is left at its default.

Value

Numeric square distance matrix, or a source-by-target matrix when sources is supplied, in requested order. Unreachable pairs are Inf. Zero-length edges are retained; duplicate points have zero distance. Point row names or graph adjacency names label the matrix. With return.graph = TRUE, return a list containing distances (this matrix), graph (a dgraph with powered base edge lengths), and metadata. Metadata records sources, targets, coverage (all_pairs or requested_pairs), p, rooted, edge.weight.type, and the backend used. Graph distances are unrooted even when the returned matrix is rooted.

Details

For base lengths \(l_e\), the unrooted result is \(\min_\pi \sum_{e\in\pi} l_e^p\). At p = 1, complete-graph coordinate distances are Euclidean; restricted distances are ordinary graph distances. Restricting edges cannot decrease the result when retained weights agree. Taking k = n - 1 recovers the complete graph.

Groisman, Jonckheere and Sapienza define the unrooted complete-graph estimator and prove consistency after sample-size normalization under their sampling assumptions. Their Proposition 2.12 establishes a high-probability equality for a sufficiently large nearest-neighbor restriction under its stated assumptions. It does not certify arbitrary fixed k, pruned graphs, or arbitrary supplied lengths. This function applies neither sample-size normalization nor an estimated limit constant.

Both coordinate constructions evaluate all Euclidean pairs. Explicit complete graphs store quadratic adjacency. The implicit backend uses dense Dijkstra searches with linear scratch and no graph restriction: with h sources, n observations and d coordinates its worst-case time is O(hn^2d). A global numeric-range check takes O(n^2d) time. Pivot-to-all output uses O(hn) storage; full all-pairs output is still quadratic. Select pivots outside this function; no automatic pivot sampling or distance approximation is performed. Distances are computed in floating-point arithmetic. Overflow, positive powered lengths underflowing to zero, and a conservative upper bound on simple-path cost exceeding the numeric range cause errors; rescale base lengths before retrying. Multiplying lengths by a > 0 multiplies unrooted distances by a^p (rooted distances by a).

With return.graph = TRUE, one deterministic shortest-path tree is selected per source (ties need not include every minimizing path). Only paths to the requested targets contribute edges. The returned graph retains all input vertices, including unused isolated vertices. Its edge lengths are powered base lengths even when the returned distances are rooted. Thus graph geodesics reproduce the unrooted requested distances. With all sources and targets it preserves every pair; a subset guarantees only the requested pairs. The union need not be sparse. The implicit workspace estimate includes a conservative allowance for the predecessor-edge union when requested.

References

Groisman, P., Jonckheere, M. and Sapienza, F. (2022). Nonhomogeneous Euclidean first-passage percolation and distance learning. Bernoulli, 28(1), 255–276. doi:10.3150/21-BEJ1341 . See Definitions 2.1 and 2.11 and Proposition 2.12 in https://arxiv.org/html/1810.09398v2.

Examples

X <- matrix(c(0, 1, 2), ncol = 1)
fermat.distances(points = X, p = 2) # 0 -> 1 -> 2 costs 2, not 4
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    0    1
#> [3,]    2    1    0
fermat.distances(points = X, p = 2, rooted = TRUE)
#>          [,1] [,2]     [,3]
#> [1,] 0.000000    1 1.414214
#> [2,] 1.000000    0 1.000000
#> [3,] 1.414214    1 0.000000
fermat.distances(points = X, graph.type = "sknn", k = 1)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    0    1
#> [3,]    2    1    0
g <- dgraph(list(2L, c(1L, 3L), 2L), list(1, c(1, 1), 1))
fermat.distances(g, vertices = c(1, 3))
#>      [,1] [,2]
#> [1,]    0    2
#> [2,]    2    0