Balancing Geometry and Density: Path Distances on High-Dimensional Data

New geometric and computational analyses of power-weighted shortest-path distances (PWSPDs) are presented. By illuminating the way these metrics balance density and geometry in the underlying data, we clarify their key parameters and discuss how they may be chosen in practice. Comparisons are made with related data-driven metrics, which illustrate the broader role of density in kernel-based unsupervised and semi-supervised machine learning. Computationally, we relate PWSPDs on complete weighted graphs to their analogues on weighted nearest neighbor graphs, providing high probability guarantees on their equivalence that are near-optimal. Connections with percolation theory are developed to establish estimates on the bias and variance of PWSPDs in the finite sample setting. The theoretical results are bolstered by illustrative experiments, demonstrating the versatility of PWSPDs for a wide range of data settings. Throughout the paper, our results require only that the underlying data is sampled from a low-dimensional manifold, and depend crucially on the intrinsic dimension of this manifold, rather than its ambient dimension.


Introduction
The analysis of high-dimensional data is a challenge in modern statistical and machine learning. In order to defeat the curse of dimensionality [38,34,10], distance metrics that efficiently and accurately capture intrinsically low-dimensional latent structure in high-dimensional data are required. Indeed, this need to capture low-dimensional linear and nonlinear structure in data has led to the development of a range of data-dependent distances and related dimension reduction methods, which have been widely employed in applications [44,57,8,26,21,58]. Understanding how these metrics trade off fundamental properties in the data (e.g. local versus global structure, geometry versus density) when making pointwise comparisons is an important challenge in their use, and may be understood as a form of model selection in unsupervised and semi-supervised machine learning problems. of the classical geodesic distance) with passing through high-density regions. We consider a finite data set X = {x i } n i=1 ⊂ R D , which we usually assume to be intrinsically low-dimensional, in the sense that there exists a compact d-dimensional Riemannian data manifold M ⊂ R D and a probability density function f (x) supported on M such that {x i } n i=1 i.i.d.
∼ f (x). Definition 1.1. For p ∈ [1, ∞) and for x, y ∈ X , the (discrete) p-weighted shortest path distance (PWSPD) from x to y is: where π is a path of points in X with x i 1 = x and x i T = y and · is the Euclidean norm.
Early uses of density-based distances for interpolation [54] led to the formulation of PWSPD in the context of unsupervised and semi-supervised learning and applications [30,60,17,53,18,13,47,46,42,64,16]. It will occasionally be useful to think of p p (·, ·) as the path distance in the complete graph on X with edge weights x i − x j p , which we shall denote G p X . When p = 1, 1 (x, y) = x − y , i.e. the Euclidean distance. As p increases, the largest elements in the set of path edge lengths x i j − x i j+1 T −1 j=1 begin to dominate the optimization (1), so that paths through higher density regions (with shorter edge lengths) are promoted. When p → ∞, p converges (up to rescaling by the number of edges achieving maximal length) to the longest-leg path distance ∞ (x, y) = min [42] and is thus driven by the density function f . Outside these extremes, p balances taking a "short" path and taking one through regions of high density. Note that p can be defined for p < 1, but it does not satisfy the triangle inequality and is thus not a metric ( p p however is a metric for all p > 0). This case was studied in [2], where it is shown to have counterintuitive properties that should preclude its use in machine learning and data analysis.
Note L 1 is simply the geodesic distance on M. However for p > 1 and a nonuniform density, the optimal path γ is generally not the geodesic distance on M: L p favors paths which travel along high-density regions, and detours off the classical L 1 geodesics are thus acceptable. The parameter p controls how large of a detour is optimal; for large p, optimal paths may become highly nonlocal and different from classical geodesic paths.
It is known [39,33] that when f is continuous and positive, for p > 1 and all x, y ∈ M, lim n→∞ n p−1 pd p (x, y) = C p,d L p (x, y) for an absolute constant C p,d depending only on p and d, i.e. that the discrete PWSPD computed on an i.i.d. sample from f (appropriately rescaled) is a consistent estimator for the continuum PWSPD.
In particular, (3) is established by [33] for C 1 , isometrically embedded manifolds and by [39] for smooth, compact manifolds without boundary and for p defined using geodesic distance. We thus define the normalized (discrete) path metric p (x, y) := n p−1 pd p (x, y) .
The n p−1 pd normalization factor accounts for the fact that for p > 1, p converges uniformly to 0 as n → ∞ [46]. Note that the 1/p exponent in (1) and (3) is necessary to obtain a metric that is homogeneous. Moreover, as p → ∞, L p is constant on regions of constant density, but L p p is not. Indeed, consider a uniform distribution on [0, 1] d , which has density f = 1 [0,1] d . Then for all x, y ∈ [0, 1] d and for all p, L p p (x, y) = x − y . On the other hand, for all x, y ∈ [0, 1] d , L p (x, y) = x − y 1/p → 1 as p → ∞, i.e. all points are equidistant in the limit p → ∞. Thus the 1/p exponent in (1) and (3) is necessary to obtain an entirely density-based metric for large p.
In practice, it is more efficient to compute PWSPDs in a sparse graph instead of a complete graph. It is thus natural to define PWSPDs with respect to a subgraph H of G p X .
Definition 1.3. Let H be any subgraph of G p X . For x, y ∈ X, let P H (x, y) be the set of paths connecting x and y in H. For p ∈ [1, ∞) and for x, y ∈ X , the (discrete) p-weighted shortest path distance (PWSPD) with respect to H from x to y is: In order to compute all-pairs PWSPDs in a complete graph with n nodes (i.e. p (x i , x j ) for all x i , x j ∈ X ), a direct application of Dijkstra's algorithm has complexity O(n 3 ). Let G p,k X denote the kNN graph, constructed from G p X by retaining only edges {x, y} if x is amongst the k nearest neighbors of y in X (we say: "x is a kNN of y" for short) or vice versa. In some cases the PWSPDs with respect to G p,k X are known to coincide with those computed in G p X [33,20]. If so, we say the kNN graph is a 1-spanner of G p X . This provides a significant computational advantage, since kNN graphs are much sparser, and reduces the complexity of computing all-pairs PWSPD to O(kn 2 ) [40].

Summary of Contributions
This article develops new analyses, computational insights, and applications of PWSPDs, which may be summarized in three major contributions. First, we establish that when p d is not too large, PWSPDs locally are density-rescaled Euclidean distances. We give precise error bounds that improve over known bounds [39] and are tight enough to prove the local equivalence of Gaussian kernels constructed with PWSPD and density-rescaled Euclidean distances. We also develop related theory which clarifies the role of density in machine learning kernels more broadly. A range of machine learning kernels that normalize in order to mitigate or leverage differences in underlying density are considered and compared to PWSPD. Relatedly, we analyze how PWSPDs become increasingly influenced by the underlying density as p → ∞. We also illustrate the role of density and benefits of PWSPDs on illustrative data sets.
Second, we improve and extend known bounds on k [33,46,20] guaranteeing that the kNN graph is a 1-spanner of G p X . Specifically, we show that for any 1 < p < ∞, the kNN graph is a 1-spanner of G p X with probability exceeding 1 − 1/n if k ≥ C p,d,f,M · log(n), for an explicit constant C p,d,f,M that depends on the density power p, intrinsic dimension d, underlying density f , and the geometry of the manifold M, but is crucially independent of n. These results are proved both in the case that the manifold is isometrically embedded and in the case that the edge lengths are in terms of intrinsic geodesic distance on the manifold. Our results provide an essential computational tool for the practical use of PWSPDs, and their key dependencies are verified numerically with extensive large-scale experiments. Third, we bound the convergence rate of PWSPD to its continuum limit using a percolation theory framework, thereby quantifying the [39,33] asymptotic convergence result (4). Specifically, we develop bias and variance estimates by relating results on Euclidean first passage percolation (FPP) to the PWSPD setting. Surprisingly, these results suggest that the variance of PWSPD is essentially independent of p, and depends on the intrinsic dimension d in complex ways. Numerical experiments verify our theoretical analyses and suggest several conjectures related to Euclidean FPP that are of independent interest.

Notation
We shall use the notation in Table 1 consistently, though certain specialized notation will be introduced as required. We assume throughout that the data X is drawn from a compact Riemannian data manifold (M, g), with additional assumptions imposed on M as needed; we do not rigorously consider the more general case that X is drawn from a distribution supported near M. If M ⊂ R D , we assume that it is isometrically embedded in R D , i.e. g is the unique metric induced by restricting the Euclidean metric on R D to M, unless otherwise stated. If an event holds with probability 1−c/n, where n = |X | and c is independent of n, we say it holds with high probability (w.h.p.).

Local Analysis: Density and Kernels
Density-driven methods are commonly used for unsupervised and semi-supervised learning [19,27,21,51,13,7,52]. Despite this popularity, the role of density is not completely clear in this context. Indeed, some methods seek to leverage variations in density while others mitigate it. In this section, we explore the role that density plays in popular machine learning kernels, including those used in self-tuning spectral clustering and diffusion maps. We compare with the effect of density in p -based kernels, and illustrate the primary advantages and disadvantages on toy data sets.

Role of Density in Graph Laplacian Kernels
A large family of algorithms [8,9,56,48,61] view data points as the nodes of a graph, and define the corresponding edge weights via a kernel function. In general, by kernel we mean a function K : R D × R D → R that captures a notion of similarity between elements of R D . More precisely, we suppose that K is of the form K(x i , x j ) = h(δ(x i , x j )) for some metric δ on R D and smooth, positive, rapidly decaying (hence integrable) function h : R → R. Our technical results will pertain exclusively to the Gaussian kernel K(x i , x j ) = exp(−δ(x i , x j ) 2 / 2 ) for some metric δ and scaling parameter > 0, albeit more general kernels have been considered in the literature [4,23,11]. Given X ⊂ R D , one first defines a weight matrix W ∈ R n×n by W ij = K(x i , x j ) for some kernel K, and diagonal degree matrix Deg ∈ R n×n by Deg ii = n j=1 W ij . A graph Laplacian L is then defined using W, Deg. Then, the K lowest frequency eigenvectors of L, denoted φ 1 , . . . , φ K , define a K-dimensional spectral embedding of the data by Commonly, a standard clustering algorithm such as K-means is then applied to the spectral embedding. This procedure is known as spectral clustering (SC). In unnormalized SC, L = Deg − W , while in normalized SC density-based stretch of Euclidean distance with respect to f W, Deg, L weight, degree, and Laplacian matrices associated to a graph δ(·, ·) arbitrary metric B δ (x, ) {y | δ(x, y) ≤ }, ball of radius > 0 centered at x with respect to δ B(x, ) Euclidean ball of radius > 0 centered at x, dimension determined by context  either the random walk Laplacian L RW = Deg −1 L or the symmetric normalized Laplacian L SYM = Deg −1/2 LDeg −1/2 is used. Many modifications of this general framework have been considered. Although SC is better able to handle irregularly shaped clusters than many traditional algorithms [5,55], it is often unstable in the presence of low degree points and sensitive to the choice of scaling parameter when using the Gaussian kernel [61]. These shortcomings motivated [63] to apply SC with the self-tuning kernel , where σ i,k is Euclidean distance of x i to its k th NN. To clarify how the data density influences this kernel, consider how σ i,k relates to the kNN density estimator at x i : It is known [43] that if k = k(n) is such that k(n) → ∞ while k(n)/n → 0, then f n (x i ) is a consistent estimator of f (x i ), as long as f is continuous and positive. Furthermore, if f is uniformly continuous and k(n)/ log n → ∞ while k(n)/n → 0, then sup i |f n (x i ) − f (x i )| → 0 with probability 1 [25]. Although these results assume the density f is supported in R d , the density estimator (5) is consistent in the general case when f is supported on a d-dimensional Riemannian manifold M ⊆ R D for log n k(n) n [28]. For such k(n), σ i,k → n,d f (x i ) − 1 d for some constant n,d depending on n, d. Thus, for n large the kernel for self-tuning spectral clustering is approximately: Relative to a standard SC kernel, (6) weakens connections in high density regions and strengthens connections in low density regions. Diffusion maps [22,21] is a more general framework which reduces to SC for certain parameter choices. More specifically, [21] considered the family of kernels parametrized by a ∈ [0, 1], which determines the degree of density normalization. Since deg( is a kernel density estimator of the density f (x i ) [12] and, up to higher order terms, Note that f has an effect on the kernel similar to the self-tuning kernel (6): connections in high density regions are weakened, and connections in low density regions are strengthened. Let L a, RW denote the discrete random walk Laplacian using the weights W ij given in (7). The discrete operator −L a, RW / 2 converges to the continuum Kolmogorov operator L ψ = ∆ψ + (2 − 2a)∇ψ · ∇f f as n → ∞, → 0 + for Laplacian operator ∆ and gradient ∇, both taken with respect to the Riemannian metric inherited from the ambient space [8,21,12]. When a = 0, we recover standard spectral clustering; there is no density renormalization in the kernel but the limiting operator is density dependent. When a = 1, −L 1, RW / 2 → ∆; in this case the discrete operator is density dependent but the limiting operator is purely geometric, since the density term is eliminated. We note that Laplacians and diffusion maps with various metrics and norms have been considered in a range of settings [62,15,59,41].

Local Characterization of PWSPD-Based Kernels
While the kernels discussed in Section 2.1 compensate for discrepancies in density, PWSPD-based kernels strengthen connections through high-density regions and weaken connections through lowdensity regions. To illustrate more clearly the role of density in PWSPD-based kernels, we first show that locally the continuum PWSPD L p p is well-approximated by the density-based stretch of , as long as f does not vary too rapidly and M does not curve too quickly. This is quantified in Lemma 2.2, which is then used to prove Theorem 2.3, which bounds the local deviation of L p from D 1/p f,Euc . Finally, Corollary 2.4 establishes that Gaussian kernels constructed with L p and D 1/p f,Euc are locally similar. Throughout this section we assume M ∈ S(d, κ, 0 ) as defined below. [4,45]. Note that this is also equivalent to a positive lower bound on the reach [29] of M (e.g. Proposition 6.1 in [49] and Proposition A.1 in [1]); see Definition 3.8.
Let B L p p (x, ) and B D (x, ) denote, respectively, the (closed) L p p and geodesic balls centered at x of radius . Let f max = max y {f (y) : y ∈ M}, f min = min y {f (y) : y ∈ M} be the global density maximum and minimum. Define the following local quantities: Let ρ x, = f max (x, )/f min (x, ), which characterizes the local discrepancy in density in a ball of radius O( ) around the point x.
The following Lemma establishes that L p p and D f,Euc are locally equivalent, and that discrepancies depend on (ρ x, ) p−1 d and the curvature constant κ. We note similar estimates appear in [2] for the special case p = 0. The proof appears in Appendix A.
Note that corresponding bounds in terms of geodesic distance follow easily from the definition of . Lemma 2.2 thus establishes that the metrics L p p and D f,Euc are locally equivalent when (i) ρ x, is close to 1, (ii) p−1 d is not too large, and (iii) κ is not too large. However, when p−1 d 1, L p p balls may become highly nonlocal in terms of geodesics.
The following Theorem establishes the local equivalence of L p and D 1/p f,Euc (and thus kernels constructed using these metrics). Assuming the density does not vary too quickly, Lemma 2.2 can be used to show that locally the difference between D 1/p f,Euc and L p is small. Variations in density are controlled by requiring that f is F -Lipschitz with respect to geodesic distance, i.e.
|f (x)−f (y)| ≤ F D(x, y). This Lipschitz assumption allows us to establish a higher-order equivalence compared to existing results (e.g. Corollary 9 in [39]), which we leverage to obtain the local kernel equivalence stated in Corollary 2.4. The following analysis also establishes explicit dependencies of the equivalence on d, p, F, κ.
Then for all x, y ∈ M such that D(x, y) ≤ 0 and x − y ≤ , Proof. We first show that ρ x, is close to 1. Let (since these sets are compact, these points must exist). Then by the Lipschitz condition: Let γ 2 (t) be a path achieving L p p (x, y 1 ). Note that . We thus obtain Letting C x, = F (ρ p−1 d x, + 1)/f min (x, ), Taylor expanding around = 0 and (10) give ρ p−1 pd x, Rewriting the above yields: We thus obtain Note the coefficient C 1 increases exponentially in p; thus the equivalence between L p and D 1/p f,Euc is weaker for large p. We also emphasize that in a Euclidean ball of radius , the metric L p scales like 1 p ; Theorem 2.3 thus guarantees that the relative error of approximating When L p is locally well-approximated by D 1/p f,Euc , the kernels constructed from these two metrics are also locally similar. The following Corollary leverages the error term in Theorem 2.3 to make this precise for Gaussian kernels. It is a direct consequence of Theorem 2.3 and Taylor expanding the Gaussian kernel, and its proof is given in is the Gaussian kernel with metric δ(·, ·) and scaling parameter > 0. Note h 1 Corollary 2.4. Under the assumptions and notation of Theorem 2.3, When p−1 is not too large relative to d, a kernel constructed with L p is locally well-approximated by a kernel constructed with D 1/p f,Euc . Thus, in a Euclidean ball of radius , we may think of the Gaussian L p kernel as: Density plays a different role in this kernel compared with those of Section 2.1. This kernel strengthens connections in high density regions and weakens them in low density regions. We note that the 1 p -power in Definition 1.2 has a large impact, in that L p -based and L p p -based kernels have very different properties. More specifically, h 1 (L p p / ) is a local kernel as defined in [12], so it is sufficient to analyze the kernel locally. However h 1 (L p / 1/p ) is a non-local kernel, so that non-trivial connections between distant points are possible. The analysis in this Section thus establishes the global equivalence of h 1 (L p p / ) and h 1 (D f,Euc / ) (when p is not too large relative to d) but only the local equivalence of h 1

The Role of p: Examples
This subsection illustrates the useful properties of PWSPDs and the role of p on three synthetic data sets in R 2 : (1) Two Rings data, consisting of two non-convex clusters that are well-separated by a low-density region; (2) Long Bottleneck data, consisting of two isotropic clusters each with a density gap connected by a long, thin bottleneck; (3) Short Bottleneck data, where two elongated clusters are connected by a short bottleneck. The data sets are shown in Figures 1, 2, and 3, respectively. We also show the PWSPD spectral embedding (denoted PWSPD SE) for various p, computed from a symmetric normalized Laplacian constructed with PWSPD. The scaling parameter for each data set is chosen as the 15 th percentile of pairwise PWSPD distances. Different aspects of the data are emphasized in the low-dimensional PWSPD embedding as p varies. Indeed, in Figure 1, we see the PWSPD embedding separates the rings for large p but not for small p. In Figure 2, we see separation across the bottleneck for p small, while for p large there is separation with respect to the density gradients that appear in the two bells of the dumbbell. Interestingly, separation with respect to both density and geometry is observed for p = 2 (see Figure  2g). In Figure 3, the clusters are both elongated and lack robust density separation, but the PWSPD embedding well-separates the two clusters for moderate p. In general, p close to 1 emphasizes the geometry of the data, large p emphasizes the density structure of the data, and moderate p defines a metric balancing these two considerations. (h)K, PWSPD SC (binarized) Figure 1: Two Rings dataset. Because the underlying cluster structure is density-driven, the PWSPD SE separates the clusters for large p (Figure 1d). While taking small in Euclidean spectral clustering can allow for good clustering accuracy (see Figure 1e), the range is narrow and does not permit accurate estimation of K via the eigengap (see Figure 1g). On the other hand, PWSPD consistently clusters well and correctly captures K = 2 for a wide range of ( , p) pairs (see Figures 1g, 1h). Generally, PWSPD allows for fully unsupervised clustering as long as p is sufficiently large and not too small.

Comparison with Euclidean Spectral Clustering
To evaluate how p impacts the clusterability of the PWSPD spectral embedding, we consider experiments in which we run spectral clustering under various graph constructions. We run K-means for a range of parameters on the spectral embedding where φ k is the k th lowest frequency eigenvector of the Laplacian. We construct the symmetric normalized Laplacian using PWSPD (denoted PWSPD SC) and also using Euclidean distances (denoted SC) and the Laplacian with diffusion maps normalization a = 1 (denoted SC+DMN). We vary in the SC and SC+DMN methods, and both and p in the PWSPD SC method. Results for selt-tuning SC, in which the kNN used to compute the local scaling parameter varies, are in Appendix D. To allow for comparisons  . When varying p, the PWSPD SE separates by geometry (Figure 2b) for p near 1, before separating by density for p 1 ( Figure  2l). Given the correct choice of and a priori knowledge of K, any of the three natural clusterings can be learned by Euclidean SC (Figure 2e, 2j). However, in the Euclidean SC case, correct estimation of K fails to coincide with parameters that give good clustering results (Figure 2o). On the other hand, PWSPD SC is able to correctly estimate each of K = 2, 3, 4 for some choice of ( , p) parameters in the same region that such parameters yield high clustering accuracy (Figures 2c, 2d   (h)K, PWSPD SC (binarized) Figure 3: Short Bottleneck dataset. Because the underlying cluster structure is not driven entirely by geometry or density, the PWSPD SE separates the clusters for moderate p (see Figure 3c). We note PWSPD is able to correctly learn K and cluster accurately for somewhat large and p between 2 and 3 (see Figures 3f, 3h), while Euclidean SC cannot simultaneously learn K and cluster accurately (see Figures 3e, 3g).
across figures, is varied across the percentiles of the pairwise distances in the underlying data, up to the 25 th percentile. We measure two outputs of the clustering experiments: (i) The overall accuracy (OA), namely the proportion of data points correctly labeled after alignment when K is known a priori. For K = 2, similar results were observed when thresholding φ 2 at 0 instead of running K-means; see Appendix D.
(ii) The eigengap estimate of the number of latent clusters:K = arg max k≥2 λ k+1 − λ k , where 0 = λ 1 ≤ λ 2 ≤ · · · ≤ λ n are the eigenvalues of the corresponding graph Laplacian. We note that experiments estimating K by considering the ratio of consecutive eigenvalues were also performed, with similar results. In the case of PWSPD SC, we plot heatmaps of where K is correctly estimated, with yellow corresponding to success (K = K) and blue corresponding to failure (K = K).
The results in terms of OA andK as a function of and p are in Figures 1, 2, 3. We see that when density separates the data clearly, as in the Two Rings data, PWSPD SC with large p gives accurate clustering results, while small p may fail. In this dataset, very small allows for the data to be correctly clustered with SC and SC+DMN when K is known a priori. However, the regime of is so small that the eigenvalues become unhelpful for estimating the number of latent clusters. Unlike Euclidean spectral clustering, PWSPD SC correctly estimatesK = 2 for a range of parameters, and achieves near-perfect clustering results for those parameters as well. Indeed, as shown by Figures 1f,  1h, PWSPD SC with p large is able to do fully unsupervised clustering on the Two Rings data.
In the case of the Long Bottleneck dataset, there are three reasonable latent clusterings, depending on whether geometry, density, or both matter (see Figure 2a, 2k, 2f). PWSPD is able to balance between the geometry and density-driven cluster structure in the data. Indeed, all of the cluster configurations shown in Figure 2a, 2k, 2f are learnable without supervision for some choice of parameters ( , p). To capture the density cluster structure (K = 3), p should be taken large, as suggested in Figure 2m, 2n. To capture the geometry cluster structure (K = 2), p should be taken small and large, as suggested by Figures 2c, 2d. Interestingly, both cluster and geometry (K = 4) can be captured by choosing p moderate, as in Figure 2h, 2i. For Euclidean SC, varying is insufficient to capture the rich structure of this data.
In the case of the Short Bottleneck, taking large allows for the Euclidean methods to correctly estimate the number of clusters. But, in this regime, the methods do not cluster accurately. On the other hand, taking p between 2 and 3 and large allows PWSPD to correctly estimate K and also cluster accurately.
Overall, this suggests that varying p in PWSPD SC has a different impact than varying the scaling parameter , and can allow for richer cluster structures to be learned when compared to SC with Euclidean distances. In addition, PWSPDs generally allow for the underlying cluster structures to be learned in a fully unsupervised manner, while Euclidean methods may struggle to simultaneously cluster well and estimate K accurately.

Spanners for PWSPD
Let H ⊂ G p X denote a subgraph and recall the definition of H p (·, ·) given in Definition 1.3.
Clearly p (x, y) ≤ H p (x, y) always, as any path in H is a path in G p X . Hence if H is a 1-spanner we have equality: H p (x, y) = p (x, y). Define the kNN graph, G p,k X , by retaining only edges {x, y} if x is a kNN of y or vice versa. For appropriate k, p and M it is known that G p,k X is a 1-spanner of G p X w.h.p. Specifically, [33] shows this when M is an open connected set with C 1 boundary, 1 < p < ∞ and k = O(c p,d log(n)) for a constant c p,d depending on p, d. One can deduce c p,d ≥ 2 d+1 3 d d d/2 , while the dependence on p is more obscure. A different approach is used in [20] to show this for arbitrary smooth, closed, isometrically embedded M, 2 ≤ p < ∞ and k = O(2 d log(n)), where O hides constants depending on the geometry of M. In both cases f must be continuous and bounded away from zero. Under these assumptions, we prove G p,k X is a 1-spanner w.h.p., for any smooth, closed, isometrically embedded M with mild restrictions on its curvature. Our results hold generally for 1 < p < ∞ and enjoy improved dependence of k on d and explicit dependence of k on p and the geometry of M compared to [33,20]. We also consider an intrinsic version of PWSPD, where D(·, ·) is assumed known, which is not typically the case in data science. However this situation can occur when X is presented as a subset of R D , but one wishes to analyze X with an exotic metric (i.e. not · ). For example, if each x i ∈ X is an image, a Wasserstein metric may be more appropriate than · . As this case closely mirrors the statement and proof of Theorem 3.9 we leave it to Appendix E. Before proceeding we introduce some further terminology:  , 0)) for p = 1.01, 2, 100. We see that for smaller p, the set becomes quite small, converging to a line segment as p → 1 + . For p = 2, the p-elongated set is a circle. As p increases, D1,p((− 1 2 , 0), ( 1 2 , 0)) converges to a set resembling a vertically-oriented American football.

Nearest Neighbors and PWSPD Spanners
A key proof ingredient is the following definition, which generalizes the role of spheres in the proof of Theorem 1.3 in [20].

Definition 3.4.
For any x, y ∈ R d and α ∈ (0, 1], the p-elongated set associated to x, y is Visualizations of D 1,p (x, y) ⊂ R 2 are shown in Figure 4. D 1,p (x, y) is the set of points z such that the two-hop path, x → z → y, is p -shorter than the one-hop path, x → y. Hence: Lemma 3.5. If there exists z ∈ D 1,p (x, y) ∩ X then the edge {x, y} is not critical.
We defer the proof of the following technical Lemma to Appendix B. For α = 1, [33] makes a similar claim but crucially does not quantify the dependence of the radius of this ball on p. Before proceeding, we introduce two regularity assumptions: Proof. In light of Lemma 3.3 we prove that, with probability at least 1 − 1/n, G p,k X contains every critical edge of G p X . Equivalently, we show every edge of G p X not contained in G p,k X is not critical. For any c, > 0, P max x,y∈X p (x, y) ≤ ≥ 1 − c/n for n sufficiently large [46]. So, let n be sufficiently large so that P p (x, y) ≤ min 0 , ζ d 1 for all x, y ∈ X ≥ 1 − 1 2n . Pick any x, y ∈ X which are not kNNs and let r := x − y . If r > min 0 , ζ  y is not a kNN of x, x − x i j ≤ x − y = r for j = 1, . . . , k. Thus, x i j ∈ B(x, r) and so by Lemma B.2 we bound for fixed j Because the x i j are all independently drawn: A routine calculations reveals that for k ≥ By Lemma 3.5 we conclude the edge {x, y} is not critical with probability exceeding 1 − 1 n 3 . There are fewer than n(n − 1)/2 such non-kNN pairs x, y ∈ X . These edges {x, y} are precisely those contained in G p X but not in G p,k X . By the union bound and (14) we conclude that none of these are critical with probability greater than 1 − n(n−1) for all x, y ∈ X , which holds with probability exceeding 1 − 1 2n . Thus, all critical edges are contained in G X p,k with probability exceeding 1 − 1 2n + 1 2n = 1 − 1 n . Unpacking ε M,p,f yields the claimed lower bound on k.
In (11), the explicit dependence of k on κ 0 , p, and d are shown. The 4κ 2 0 factor corresponds to the geometry of M. The numerical constant 4, which is not tight, stems from accounting for the reach of M. If M is convex (i.e. ζ = ∞) then it can be replaced with 3. The second factor in (11) is controlled by the probability distribution while the third corresponds to p and d. For p = 2 and ignoring geometric and density factors we attain k = O(2 d log(n)) as in [20]. For large p we get k ≈ O 4 quantitative lower bound n ≥ C max d ζ pd p−1 for a constant C depending on the geometry of M. So, when M is high-dimensional, has small reach, or when p is close to 1, n may need to be quite large for k as in (11) to yield a 1-spanner.

Numerical Experiments
We verify the claimed dependence of k on n, p and d ensures that G p,k X is a 1-spanner of G p X numerically. To generate Figures 5a-5f we: (1) Fix p, d, M, and f , then generate a sequence of (n, k) pairs.
(3) Repeat step 2 twenty times and compute the proportion of successes.
As can be seen from Figure 5, there is a sharp transition between an "all failures" and an "all successes" regime. The transition line is roughly linear when viewed using semi-log-x axes, i.e. k ∝ log(n). Moreover the slope of the line-of-best-fit to this transition line decreases with increasing p (compare Figure 5a-5c) and depends on intrinsic, not extrinsic dimension (compare Figure 5b and 5d), as predicted by Theorem 3.9. Intriguingly, there is little difference between Figure 5b (uniform distribution) and Figure 5e (Gaussian distribution), suggesting that perhaps the assumption f min > 0 in Theorem 3.9 is unnecessary. Finally, we observe that the constant of proportionality (i.e. C such that k = C log n) predicted by Theorem 3.9 appears pessimistic. For Figure 5a-5c, Theorem 3.9 predicts C = 484.03, 128 and 21.76 respectively (taking κ 0 = 1 due to the flat domain), while empirically the slope of the line-of-best-fit is 43.43, 25.33 and 0.29 respectively.
In Figure 5f, we consider an intrinsically 4-dimensional set corrupted with Gaussian noise (standard deviation 0.1) in the fifth dimension. Interestingly, the scaling with k is more efficient than as shown in Figure 5a for the intrinsically 5-dimensional data. This suggests that measures which concentrate near low-dimensional sets benefit from that low-dimensionality, even if they are not supported exactly on it.
We also consider relaxing the success condition (2.iii). We define H to be a (t, ω)-spanner if H p (x, y) ≤ t p (x, y) for ω ∈ (0, 1] proportion of the edges, so that Theorem 3.9 pertains to (1,1)spanners. Figures 5g and 5h show the minimal ω (averaged across simulations) for which G p,k X is a (1.1, ω)-spanner and a (1.01, ω)-spanner respectively; the red lines trace out the requirements for G p,k X to be a (1.1, 1)-spanner and (1.01, 1)-spanner respectively. Comparing with Figure 5b, we see that the required scaling for G p,k X to be a (1 + , 1)-spanner is similar to the required scaling to be a (1, 1)-spanner, at least for > 0 small. However, the required scaling for (1 + , ω)-spanners (ω < 1) is quite different and much less restrictive, even for ω very close to 1; for example the requirement for G p,k X to be a (1.01, 0.95)-spanner appears sublinear in the log 2 (n) versus k plot (see Figure 5h). If this notion of approximation is acceptable, our empirical results suggest one can enjoy much greater sparsity. Finally, in Figure 5i we compute the minimal t ≥ 1 such that G p,k X is a (t, 1)-spanner of G p X ; again the overall transition patterns for (t, 1)-spanners are similar to the (1, 1)-spanner case in Figure 5b when t is close to 1. Overall we see that greater sparsity is permissible in these relaxed cases, and analyzing such notions rigorously is a topic of ongoing research.   Figures 5a-5f show the proportion of randomly generated data sets for which G p,k X is a 1-spanner of G p X . The red line is the line of best fit through the cells representing the first value of k, for each value of n, for which all trials were successful, i.e. it is the line ensuring G p,k X is a 1-spanner. The slopes of Figure 5a-5f are, respectively, 43.43, 25.33, 0.29, 14.18, 18.79, and 40.0. Figures 5g and 5h show the minimal ω (averaged across simulations) for which G p,k X is a (1.1, ω)-spanner and a (1.01, ω)-spanner respectively; the red lines trace out the requirements for G p,k X to be a (1.1, 1)-spanner and (1.01, 1)-spanner respectively. Figure 5i shows the minimal t ≥ 1 such that G p,k X is a (t, 1)-spanner of G p X , and the red line traces out the (1, 1)-spanner requirement.

Global Analysis: Statistics on PWSPD and Percolation
We recall that after a suitable normalization, p is a consistent estimator for L p . Indeed, [39,33] prove that for any d ≥ 1, p > 1, there exists a constant C p,d independent of n such that lim n→∞˜ p (x, y) = C p,d L p (x, y). The important question then arises: how quickly does p converge? How large does n need to be to guarantee the error incurred by approximating L p with p is small? To answer this question we turn to results from Euclidean first passage percolation (FPP) [36,37,6,24]. For any discrete set X , we let p (x, y, X ) denote the PWSPD computed in the set X ∪ {x} ∪ {y}.

Overview of Euclidean First Passage Percolation
Euclidean FPP analyzes p p (0, z, H 1 ), where H 1 is a homogeneous, unit intensity Poisson point process (PPP) on R d . It is known that where µ = µ p,d is a constant depending only on p, d known as the time constant. The convergence of p p (0, z, H 1 ) is studied by decomposing the error into random and deterministic fluctuations, i.e.
In terms of mean squared error (MSE), one has the standard bias-variance decomposition: The following Proposition is well known in the Euclidean FPP literature.
constant C depending only on p, d.
Proof. By Theorem 2.1 in [37], Var [ p p (0, z, H 1 )] ≤ C z . By Theorem 2.1 in [3], Although Var [ p p (0, z, H 1 )] ≤ C z is the best bound which has been proved, the fluctuation rate is known to in fact depend on the dimension, i.e. Var [ p p (0, z, H 1 )] ∼ z 2χ for some exponent χ = χ(d) ≤ 1 2 . Strong evidence is provided in [24] that the bias can be bounded by the variance, so the exponent χ very likely controls the total convergence rate.
The following tail bound is also known [37].

Convergence Rates for PWSPD
We wish to utilize the results in Section 4.1 to obtain convergence rates for PWSPD. However, we are interested in PWSPD computed on a compact set with boundary M and the convergence rate of p rather than p p . To simplify the analysis, we restrict our attention to the following idealized model. We establish three things: (i) Euclidean FPP results apply away from ∂M ; (ii) the time constant µ equals the constant C p,d in (3); (iii) p has the same convergence rate as p p . To establish (i), we let H n denote a homogeneous PPP with rate λ = n, and let p (0, y, H n ) denote the length of the shortest path connecting 0 and y in H n . We also let X N = H n ∩ M and p (0, y, X N ) denote the PWSPD in X N ; note E[|X N |] = n. To apply percolation results to our setting, the statistical equivalence of p (0, y, X n ), p (0, y, X N ), and p (0, y, H n ) must be established. For n large, the equivalence of p (0, y, X n ) and p (0, y, X N ) is standard and we omit any analysis. The equivalence of p (0, y, X N ) and p (0, y, H n ) is less clear. In particular, how far away from ∂M do 0, y need to be to ensure these metrics are the same? The following Proposition is a direct consequence of Theorem 2.4 from [37], and essentially guarantees the equivalence of the metrics as long as 0 and y are at least distance O(n − 1 4d ) from ∂M .
Next we establish the equivalence of µ p,d (percolation time constant) and C p,d (PWSPD discreteto-continuum normalization constant).

Conclusion and Future Work
This article establishes local equivalence of PWSPD to a density-based stretch of Euclidean distance. We derive a near-optimal condition on k for the kNN graph to be a 1-spanner for PWSPD, quantifying and improving the dependence on p and d. Moreover, we leverage the theory of Euclidean FPP to establish statistical convergence rates for PWSPD to its continuum limit, and apply our spanner results to empirically support conjectures on the optimal dimension-dependent rates of convergence. Many directions remain for future work. Our statistical convergence rates for PWSPD in Section 4 are limited to uniform distributions. Preliminary numerical experiments indicate that these rates also hold for PWSPDs defined with varying density, but rigorous convergence rates for nonhomogeneous PPPs are lacking in the literature.
The analysis of Section 2 proved the local equivalence of PWSPDs with density-stretched Euclidean distances. These results and the convergence results of Section 4 are the first steps in a program of developing a discrete-to-continuum limit analysis for PWSPDs and PWSPD-based operators. A major problem is to develop conditions so that the discrete graph Laplacian (defined with p ) converges to a continuum second order differential operator as n → ∞. A related direction is the analysis of how data clusterability with PWSPDs depends on p for various random data models and in specific applications.
The numerical results of Section 3.2 confirm that k ∝ log(n) is required for the kNN graph to be a 1-spanner, as predicted by theory. Relaxing the notion of t-spanners to (t, ω)-spanners, as suggested in Section 3.2, is a topic of future research.
Finally, the results of this article require data to be generated from a distribution supported exactly on a low-dimensional manifold M. An arguably more realistic setting is the noisy one in which the data is distributed only approximately on M. Two potential models are of interest: (i) replacing M with B(M, τ ) = {x ∈ R D | dist(x, M) ≤ τ } (tube model) and (ii) considering a density that concentrates on M, rather than being supported on it (concentration model). PWSPDs may exhibit very different properties under these two noise models, for example under bounded uniform noise and Gaussian noise, especially for large p. For the concentration model one expects noisy PWSPDs to converge to manifold PWSPDs for p large, since the optimal PWSPD paths are density driven. Preliminary empirical results (Figure 5f) suggest that when the measure concentrates sufficiently near a low-dimensional set M, the number of nearest neighbors needed for a 1-spanner benefits from the intrinsic low-dimensional structure. For the tube model, although noisy PWSPDs will not converge to manifold PWSPDs, they will still scale according to the intrinsic manifold dimension for τ small. For both models, incorporating a denoising procedure such as local averaging [31] or diffusion [35] before computing PWSPDs is expected to be advantageous. Future research will investigate robust denoising procedures for PWSPD, computing PWSPDs after dimension reduction, and which type of noise distributions are most adversarial to PWSPD. By Definition 3.7, B(x M ,r 2 )∩M f ≥ f min vol (B(x M , r 2 ) ∩ M) ≥ f min κ −1 0 (r 2 ) d vol (B(0, 1)) and B(x,r)∩M f ≤ f max vol (B(x, r) ∩ M) ≤ f max κ 0 r d vol (B(0, 1)). Returning to (20): The result follows by noting . C Local Analysis: Proof of Corollary 2.4 Proof. First note that by Theorem 2.3

D.1 Spectral Clustering with Self-Tuning
Clustering results for Euclidean SC with self-tuning (ST) are shown in Figure 7. Similarly to Euclidean SC and Euclidean SC with the diffusion maps normalization, SC+ST can cluster well given K a priori but struggles to simultaneously learn K and cluster accurately.

D.2 Clustering with Fiedler Eigenvector
Clustering results on the Two Rings and Short Bottleneck data (both of which have K = 2) appear in Figure 8. The results are very similar to running K-means on the second eigenvector of L.

E PWSPD Spanners: Intrinsic Path Distance
In this section we assume that M is a compact Riemannian manifold with metric g, but we do not assume that M is isometrically embedded in R D . Let us first establish some notation. For precise definitions of terms in italics, we refer the reader to [50].
• For any v, w ∈ T x M, R x (v, w) denotes the Riemannian curvature while K x (v, w) denotes the sectional curvature. In this notation, R x (v, w) is a matrix while K x (v, w) is a scalar. These notions of curvature are related: We shall drop the "x" in the subscript when it is clear from context. Let K min be such that V min (r) = 1 − S max 6(d + 2) V max (r) = 1 − S min 6(d + 2) r 2 + O(r 4 ) vol(B(0, 1))r d .
• For any x ∈ M, exp x : T x M → M denotes the exponential map. By construction, D(x, exp x (v)) = g x (v, v). The exponential map is used to construct normal coordinates.
• The injectivity radius at x is denoted by inj(x), while inj(M) denotes the injectivity radius of M. Because M is closed, the injectivity radius is bounded away from 0.
Proposition E.1. In normal coordinates the metric has the following expansion: where δ kl = 1 if k = l and is zero otherwise. Moreover, for any v, w ∈ B(0, 1) and r ≤ inj(M): Proof. See the discussion in [50] above Theorem 5.5.8 and Exercises 5.9.26 and 5.9.27.
Combining (22) and (25) yields: For any x i ∈ X let r i,k denote the distance to the kNN of x i . Because r i,k → 0 + as n → ∞, we have that r i,k ≤ inj(M) for all x i , almost surely for n large enough.
The proof of Theorem E.2 proceeds in a similar manner to the proof of Theorem 3.9. However, care must be taken to account for the curvature. Note that for any quantity c := c(n) implicitly dependent on n, we shall say that c = o(1) if for any > 0 there exists an N such that for all n > N we have that c(n) < . Let G p M,X denote the complete graph on X with edge weights D(x i , x j ) p while G p,k M,X shall denote the kNN subgraph of G p M,X .
G p,k M,X is a 1-spanner of G p M,X with probability at least 1 − 1/n.
Proof. As in Theorem 3.9, we prove this by showing that every edge of G p M,X not contained in G p,k M,X is not critical. As before, w.h.p. M,p (x, y) → 0 uniformly in x, y [39]. So, let n be large enough so that P [ M,p (x, y) ≤ inj(M) for all x, y ∈ X ] ≥ 1 − 1 2n . Pick any x, y ∈ X which are not kNNs and let r := D(x, y). If r > inj(M), then p (x, y) < D(x, y) and thus the edge {x, y} is not critical. So, suppose without loss of generality in what follows that r ≤ inj(M). Let x i 1 , . . . , x i k denote the kNNs of x. Because y is not a kNN of x, D(x, x i j ) ≤ D(x, y) = r for j = 1, . . . , k. We show {x, y} is not critical by showing there exists an ∈ {1, . . . , k} such that D(x, x i ) p + D(x i , y) p ≤ D(x, y) p with probability at least 1 − 1/n.
For both cases, consider the p-elongated set defined in the tangent space: As in Theorem 3.9: . (by definition of V min and V max ) Using (23), (24) and r = o(1): .
If x i j ∈ exp (D 1−α,p ) then from (27) (resp. (29)): and so {x, y} is not critical. Thus, {x, y} is not critical with probability exceeding 1− 1 n 3 . These edges {x, y} are precisely those contained in G p M,X but not in G p,k M,X . There are fewer than n(n − 1)/2 such non k-NN pairs x, y ∈ X . By the union bound and (33) we conclude that none of these are critical with probability greater than 1 − n(n−1)  Table 3 shows confidence interval estimates for χ obtained by computing˜ p in a sparse graph.  Table 3: Confidence interval estimates of χ for uniform data for different density weightings (p) and different dimensions (d).