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

In this paper we analyze power-weighted shortest path distances (PWSPDs) and develop their applications to problems in machine learning. These metrics compute the shortest path between two points in the data, accounting for the underlying density of the points along the path. Paths through low-density regions are penalized, so that the optimal path must balance being “short” (in the sense of the classical geodesic distance) with passing through high-density regions. We consider a finite data set $\mathcal{X}=\{x_{i}\}_{i=1}^{n}\subset\mathbb{R}^{D}$, which we usually assume to be intrinsically low-dimensional, in the sense that there exists a compact $d$-dimensional Riemannian data manifold $\mathcal{M}\subset\mathbb{R}^{D}$ and a probability density function $f(x)$ supported on $\mathcal{M}$ such that $\{x_{i}\}_{i=1}^{n}\mathbin{\overset{i.i.d.}{\kern 0.0pt\leavevmode\resizebox{17.63316pt}{3.66875pt}{$\sim$}}}f(x)$.

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 $\ell_{p}^{p}(\cdot,\cdot)$ as the path distance in the complete graph on $\mathcal{X}$ with edge weights $\|x_{i}-x_{j}\|^{p}$, which we shall denote $\mathcal{G}_{\mathcal{X}}^{p}$. When $p=1$, $\ell_{1}(x,y)=\|x-y\|$, i.e. the Euclidean distance. As $p$ increases, the largest elements in the set of path edge lengths $\left\{\left\|x_{i_{j}}-x_{i_{j+1}}\right\|\right\}_{j=1}^{T-1}$ begin to dominate the optimization (1), so that paths through higher density regions (with shorter edge lengths) are promoted. When $p\rightarrow\infty$, $\ell_{p}$ converges (up to rescaling by the number of edges achieving maximal length) to the longest-leg path distance $\displaystyle\ell_{\infty}(x,y)=\min_{\pi=\{x_{i_{j}}\}_{j=1}^{T}}\max_{j=1,\dots,T-1}\|x_{i_{j}}-x_{i_{j+1}}\|$ [42] and is thus driven by the density function $f$. Outside these extremes, $\ell_{p}$ balances taking a “short” path and taking one through regions of high density. Note that $\ell_{p}$ can be defined for $p<1$, but it does not satisfy the triangle inequality and is thus not a metric ($\ell_{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.

While (1) is defined for finite data, it admits a corresponding continuum formulation.

Note $\mathcal{L}_{1}$ is simply the geodesic distance on $\mathcal{M}$. However for $p>1$ and a nonuniform density, the optimal path $\gamma$ is generally not the geodesic distance on $\mathcal{M}$: $\mathcal{L}_{p}$ favors paths which travel along high-density regions, and detours off the classical $\mathcal{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\in\mathcal{M}$,

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 $\ell_{p}$ defined using geodesic distance. We thus define the normalized (discrete) path metric

The $n^{\frac{p-1}{pd}}$ normalization factor accounts for the fact that for $p>1$, $\ell_{p}$ converges uniformly to 0 as $n\rightarrow\infty$ [46]. Note that the $1/p$ exponent in (1) and (3) is necessary to obtain a metric that is homogeneous. Moreover, as $p\rightarrow\infty$, $\mathcal{L}_{p}$ is constant on regions of constant density, but $\mathcal{L}_{p}^{p}$ is not. Indeed, consider a uniform distribution on $[0,1]^{d}$, which has density $f=\mathbbm{1}_{[0,1]^{d}}$. Then for all $x,y\in[0,1]^{d}$ and for all $p$, $\mathcal{L}_{p}^{p}(x,y)=\|x-y\|$. On the other hand, for all $x,y\in[0,1]^{d}$, $\mathcal{L}_{p}(x,y)=\|x-y\|^{1/p}\rightarrow 1$ as $p\rightarrow\infty$, i.e. all points are equidistant in the limit $p\rightarrow\infty$. 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 $\mathcal{H}$ of $\mathcal{G}_{\mathcal{X}}^{p}$.

Clearly $\ell_{p}^{\mathcal{G}_{\mathcal{X}}^{p}}(\cdot,\cdot)=\ell_{p}(\cdot,\cdot)$. In order to compute all-pairs PWSPDs in a complete graph with $n$ nodes (i.e. $\ell_{p}(x_{i},x_{j})$ for all $x_{i},x_{j}\in\mathcal{X}$), a direct application of Dijkstra’s algorithm has complexity $O(n^{3})$. Let $\mathcal{G}^{p,k}_{\mathcal{X}}$ denote the $k$NN graph, constructed from $\mathcal{G}^{p}_{\mathcal{X}}$ by retaining only edges $\{x,y\}$ if $x$ is amongst the $k$ nearest neighbors of $y$ in $\mathcal{X}$ (we say: “$x$ is a $k$NN of $y$” for short) or vice versa. In some cases the PWSPDs with respect to $\mathcal{G}^{p,k}_{\mathcal{X}}$ are known to coincide with those computed in $\mathcal{G}_{\mathcal{X}}^{p}$ [33, 20]. If so, we say the $k$NN graph is a $1$-spanner of $\mathcal{G}_{\mathcal{X}}^{p}$. This provides a significant computational advantage, since $k$NN graphs are much sparser, and reduces the complexity of computing all-pairs PWSPD to $O(kn^{2})$ [40].

This article develops new analyses, computational insights, and applications of PWSPDs, which may be summarized in three major contributions. First, we establish that when $\frac{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\rightarrow\infty$. 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 $k$NN graph is a $1$-spanner of $\mathcal{G}_{\mathcal{X}}^{p}$. Specifically, we show that for any $1<p<\infty$, the $k$NN graph is a $1$-spanner of $\mathcal{G}_{\mathcal{X}}^{p}$ with probability exceeding $1-1/n$ if $k\geq C_{p,d,f,\mathcal{M}}\cdot\log(n)$, for an explicit constant $C_{p,d,f,\mathcal{M}}$ that depends on the density power $p$, intrinsic dimension $d$, underlying density $f$, and the geometry of the manifold $\mathcal{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.

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

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 $\mathcal{K}:\mathbb{R}^{D}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$ that captures a notion of similarity between elements of $\mathbb{R}^{D}$. More precisely, we suppose that $\mathcal{K}$ is of the form $\mathcal{K}(x_{i},x_{j})=h(\delta(x_{i},x_{j}))$ for some metric $\delta$ on $\mathbb{R}^{D}$ and smooth, positive, rapidly decaying (hence integrable) function $h:\mathbb{R}\rightarrow\mathbb{R}$. Our technical results will pertain exclusively to the Gaussian kernel $\mathcal{K}(x_{i},x_{j})=\exp(-\delta(x_{i},x_{j})^{2}/\epsilon^{2})$ for some metric $\delta$ and scaling parameter $\epsilon>0$, albeit more general kernels have been considered in the literature [4, 23, 11]. Given $\mathcal{X}\subset\mathbb{R}^{D}$, one first defines a weight matrix $W\in\mathbb{R}^{n\times n}$ by $W_{ij}=\mathcal{K}(x_{i},x_{j})$ for some kernel $\mathcal{K}$, and diagonal degree matrix $\text{Deg}\in\mathbb{R}^{n\times n}$ by $\text{Deg}_{ii}=\sum_{j=1}^{n}W_{ij}$. A graph Laplacian $L$ is then defined using $W,\text{Deg}$. Then, the $K$ lowest frequency eigenvectors of $L$, denoted $\phi_{1},\ldots,\phi_{K}$, define a $K$-dimensional spectral embedding of the data by $x_{i}\mapsto(\phi_{1}(x_{i}),\phi_{2}(x_{i}),\dots,\phi_{K}(x_{i}))$, where $\phi_{j}(x_{i})=(\phi_{j})_{i}$. 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=\text{Deg}-W$, while in normalized SC either the random walk Laplacian $L_{\text{RW}}=\text{Deg}^{-1}L$ or the symmetric normalized Laplacian $L_{\text{SYM}}=\text{Deg}^{-1/2}L\text{Deg}^{-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 $\epsilon$ when using the Gaussian kernel [61]. These shortcomings motivated [63] to apply SC with the self-tuning kernel $W_{ij}=\exp\left(-\frac{\|x_{i}-x_{j}\|^{2}}{\sigma_{i,k}\sigma_{j,k}}\right),$ where $\sigma_{i,k}$ is Euclidean distance of $x_{i}$ to its $k^{\text{th}}$ NN. To clarify how the data density influences this kernel, consider how $\sigma_{i,k}$ relates to the $k$NN density estimator at $x_{i}$:

It is known [43] that if $k=k(n)$ is such that $k(n)\rightarrow\infty$ while $k(n)/n\rightarrow 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\rightarrow\infty$ while $k(n)/n\rightarrow 0$, then $\sup_{i}|f_{n}(x_{i})-f(x_{i})|\rightarrow 0$ with probability 1 [25]. Although these results assume the density $f$ is supported in $\mathbb{R}^{d}$, the density estimator (5) is consistent in the general case when $f$ is supported on a $d$-dimensional Riemannian manifold $\mathcal{M}\subseteq\mathbb{R}^{D}$ for $\log n\ll k(n)\ll n$ [28]. For such $k(n)$, $\sigma_{i,k}\rightarrow\epsilon_{n,d}f(x_{i})^{-\frac{1}{d}}$ for some constant $\epsilon_{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\in[0,1]$, which determines the degree of density normalization. Since $\text{deg}(x_{i})\propto f(x_{i})+O(\epsilon^{2})$, $\text{deg}(x_{i})$ 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_{\text{RW}}^{a,\epsilon}$ denote the discrete random walk Laplacian using the weights $W_{ij}$ given in (7). The discrete operator $-L_{\text{RW}}^{a,\epsilon}/\epsilon^{2}$ converges to the continuum Kolmogorov operator $\mathscr{L}\psi=\Delta\psi+(2-2a)\nabla\psi\cdot\frac{\nabla f}{f}$ as $n\rightarrow\infty,\epsilon\rightarrow 0^{+}$ for Laplacian operator $\Delta$ and gradient $\nabla$, 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_{\text{RW}}^{1,\epsilon}/\epsilon^{2}\rightarrow\Delta$; 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].

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 low-density regions. To illustrate more clearly the role of density in PWSPD-based kernels, we first show that locally the continuum PWSPD $\mathcal{L}_{p}^{p}$ is well-approximated by the density-based stretch of Euclidean distance $\displaystyle\mathscr{D}_{f,\text{Euc}}(x,y)=\frac{\|x-y\|}{\left(f(x)f(y)\right)^{\frac{p-1}{2d}}},$ as long as $f$ does not vary too rapidly and $\mathcal{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 $\mathcal{L}_{p}$ from $\mathscr{D}_{f,\text{Euc}}^{1/p}$. Finally, Corollary 2.4 establishes that Gaussian kernels constructed with $\mathcal{L}_{p}$ and $\mathscr{D}_{f,\text{Euc}}^{1/p}$ are locally similar. Throughout this section we assume $\mathcal{M}\in S(d,\kappa,\epsilon_{0})$ as defined below.

The condition $\mathscr{D}(x,y)\leq\|x-y\|(1+\kappa\|x-y\|^{2})$ for all $x,y\in\mathcal{M}$ such that $\mathscr{D}(x,y)\leq\epsilon_{0}$ is equivalent to an upper bound on the second fundamental form: $\|II_{x}\|\leq\kappa$ for all $x\in\mathcal{M}$ [4, 45]. Note that this is also equivalent to a positive lower bound on the reach [29] of $\mathcal{M}$ (e.g. Proposition 6.1 in [49] and Proposition A.1 in [1]); see Definition 3.8.

Let $B_{\mathcal{L}_{p}^{p}}(x,\epsilon)$ and $B_{\mathscr{D}}(x,\epsilon)$ denote, respectively, the (closed) $\mathcal{L}_{p}^{p}$ and geodesic balls centered at $x$ of radius $\epsilon$. Let $f_{\text{max}}=\max_{y}\{f(y):y\in\mathcal{M}\}$, $f_{\text{min}}=\min_{y}\{f(y):y\in\mathcal{M}\}$ be the global density maximum and minimum. Define the following local quantities:

Let $\rho_{x,\epsilon}=f_{\text{max}}(x,\epsilon)/f_{\text{min}}(x,\epsilon)$, which characterizes the local discrepancy in density in a ball of radius $O(\epsilon)$ around the point $x$.

The following Lemma establishes that $\mathcal{L}_{p}^{p}$ and $\mathscr{D}_{f,\text{Euc}}$ are locally equivalent, and that discrepancies depend on $(\rho_{x,\epsilon})^{\frac{p-1}{d}}$ and the curvature constant $\kappa$. 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 $\mathcal{L}_{p}$: $f_{\max}(x,\epsilon)^{-\frac{p-1}{d}}\mathscr{D}(x,y)\leq\mathcal{L}^{p}_{p}(x,y)\leq f_{\min}(x,\epsilon)^{-\frac{p-1}{d}}\mathscr{D}(x,y)$. Lemma 2.2 thus establishes that the metrics $\mathcal{L}_{p}^{p}$ and $\mathscr{D}_{f,\text{Euc}}$ are locally equivalent when (i) $\rho_{x,\epsilon}$ is close to 1, (ii) $\frac{p-1}{d}$ is not too large, and (iii) $\kappa$ is not too large. However, when $\frac{p-1}{d}\gg 1$, $\mathcal{L}_{p}^{p}$ balls may become highly nonlocal in terms of geodesics.

The following Theorem establishes the local equivalence of $\mathcal{L}_{p}$ and $\mathscr{D}_{f,\text{Euc}}^{1/p}$ (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 $\mathscr{D}_{f,\text{Euc}}^{1/p}$ and $\mathcal{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)|\leq F\mathscr{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,\kappa$.

Note the coefficient $C_{1}$ increases exponentially in $p$; thus the equivalence between $\mathcal{L}_{p}$ and $\mathscr{D}_{f,\text{Euc}}^{1/p}$ is weaker for large $p$. We also emphasize that in a Euclidean ball of radius $\epsilon$, the metric $\mathcal{L}_{p}$ scales like $\epsilon^{\frac{1}{p}}$; Theorem 2.3 thus guarantees that the relative error of approximating $\mathcal{L}_{p}$ with $\mathscr{D}_{f,\text{Euc}}^{1/p}$ is $O(\epsilon)$.

When $\mathcal{L}_{p}$ is locally well-approximated by $\mathscr{D}_{f,\text{Euc}}^{1/p}$, 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 Appendix C. Let $h_{a}(x)=\exp(-x^{2a})$ so that $h_{1}\left(\frac{\delta(\cdot,\cdot)}{\epsilon}\right)$ is the Gaussian kernel with metric $\delta(\cdot,\cdot)$ and scaling parameter $\epsilon>0$. Note $h_{1}\left(\frac{\mathcal{L}_{p}}{\epsilon^{1/p}}\right)=h_{\frac{1}{p}}\left(\frac{\mathcal{L}_{p}^{p}}{\epsilon}\right)$.

When $p-1$ is not too large relative to $d$, a kernel constructed with $\mathcal{L}_{p}$ is locally well-approximated by a kernel constructed with $\mathscr{D}_{f,\text{Euc}}^{1/p}$. Thus, in a Euclidean ball of radius $\epsilon$, we may think of the Gaussian $\mathcal{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 $\frac{1}{p}$-power in Definition 1.2 has a large impact, in that $\mathcal{L}_{p}$-based and $\mathcal{L}_{p}^{p}$-based kernels have very different properties. More specifically, $h_{1}(\mathcal{L}_{p}^{p}/\epsilon)$ is a local kernel as defined in [12], so it is sufficient to analyze the kernel locally. However $h_{1}(\mathcal{L}_{p}/\epsilon^{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}(\mathcal{L}_{p}^{p}/\epsilon)$ and $h_{1}(\mathscr{D}_{f,\text{Euc}}/\epsilon)$ (when $p$ is not too large relative to $d$) but only the local equivalence of $h_{\frac{1}{p}}(\mathcal{L}_{p}^{p}/\epsilon)$ and $h_{\frac{1}{p}}(\mathscr{D}_{f,\text{Euc}}/\epsilon)$.

This subsection illustrates the useful properties of PWSPDs and the role of $p$ on three synthetic data sets in $\mathbb{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 $\epsilon$ for each data set is chosen as the $15^{\text{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 2(g)). 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.