Power Weighted Shortest Paths for Clustering Euclidean Data

TL;DR

This paper introduces power weighted shortest path distances for clustering high-dimensional Euclidean data, demonstrating improved accuracy and providing an efficient algorithm for computation.

Contribution

It proposes a novel distance function for clustering that leverages power weights, with theoretical justification and a fast algorithm for practical use.

Findings

Higher clustering accuracy with power weighted shortest paths

Theoretical support for the effectiveness of the method

Efficient algorithm for computing the distances

Abstract

We study the use of power weighted shortest path distance functions for clustering high dimensional Euclidean data, under the assumption that the data is drawn from a collection of disjoint low dimensional manifolds. We argue, theoretically and experimentally, that this leads to higher clustering accuracy. We also present a fast algorithm for computing these distances.

Figures

Tables

Table 1. Table 1: Classification accuracy of spectral clustering. Note that A ( 1 ) superscript 𝐴 1 A^{(1)} represents using the Euclidean metric.

	$A^{(f,1)}$	$A^{\left(1\right)}$	$A^{\left(2\right)}$	$A^{\left(10\right)}$	$A^{\left(\infty \right)}$
3 Lines	$66.11\pm 0.94\%$	$66.35\pm 3.73\%$	$66.87\pm 3.37\%$	$95.38\pm 9.22\%$	$95.38\pm 9.1\%$
3 Moons	$85.90\pm 1.13\%$	$94.40\pm 1.48\%$	$94.40\pm 1.48\%$	$96.20\pm 1.76\%$	$94.35\pm 3.34\%$
3 Circles	$51.87\pm 0.00\%$	$51.93\pm 0.32\%$	$51.94\pm 0.36\%$	$71.22\pm 9.50\%$	$73.61\pm 10.47\%$
DrivFace	$78.88\%$	$71.62\%$	$71.62\%$	$74.71\%$	$85.38\%$
COIL-20	$63.24\%$	$75.28\%$	$78.61\%$	$77.45\%$	$60.92\%$
OptDigits	$77.73\%$	$91.49\%$	$91.54\%$	$88.39\%$	$83.17\%$
USPS	$48.65\%$	$65.05\%$	$65.02\%$	$76.20\%$	$77.92\%$
MNIST	-	$76.11\%$	$75.63\%$	$84.54\%$	$86.77\%$

Table 2. Table 2: Run time of spectral clustering, in seconds. Note that this includes the time rquired to construct the similarity matrix. A ( 1 ) superscript 𝐴 1 A^{(1)} represents using the Euclidean metric.

	$A^{(f,1)}$	$A^{\left(1\right)}$	$A^{\left(2\right)}$	$A^{\left(10\right)}$	$A^{\left(\infty \right)}$
3 Lines	$0.32$	$0.16$	$1.20$	$1.22$	$1.22$
3 Moons	$0.33$	$0.17$	$1.31$	$1.30$	$1.36$
3 Circles	$0.35$	$0.16$	$1.00$	$1.06$	$1.07$
DrivFace	$0.37$	$1.24$	$1.55$	$1.64$	$1.64$
COIL-20	$0.57$	$0.72$	$1.57$	$1.82$	$1.78$
OptDigits	$5.40$	$1.41$	$5.28$	$5.58$	$5.67$
USPS	$27.40$	$17.12$	$26.75$	$22.78$	$23.79$
MNIST	-	$2060.23$	$2031.38$	$1554.15$	$1613.41$

Equations

$$A_{ij}=\left\{\begin{array}[]{cc}1&\text{ if ${\bf x}_{i}$ among the $k$ nearest neighbors of ${\bf x}_{j}$ with respect to $d(\cdot,\cdot)$}\\ 0&\text{ otherwise }\end{array}\right.$$

$${\text{dist}}_a(\mathbf{x},\mathbf{y}):=\underset{\gamma }{\inf }{\int }_0^1\sqrt{g_a({\gamma }^{{}^{\prime }}(t),{\gamma }^{{}^{\prime }}(t))}dt$$

$$\text{diam}({\mathcal{M}}_a):=\underset{\mathbf{x},\mathbf{y}\in {\mathcal{M}}_a}{\sup }{\text{dist}}_a(\mathbf{x},\mathbf{y})$$

$$\text{dist}({\mathcal{M}}_a,{\mathcal{M}}_b)=\underset{\mathbf{x}\in {\mathcal{M}}_a,\mathbf{y}\in {\mathcal{M}}_b}{\min }\Vert \mathbf{x}-\mathbf{y}\Vert \ge \delta >0\text{ for all},1\le a<b\le \ell$$

$$B({\mathcal{M}}_a,\tau )=\left\{\mathbf{x}\in {\mathbb{R}}^D:\underset{\mathbf{y}\in {\mathcal{M}}_a}{\min }\Vert \mathbf{x}-\mathbf{y}{\Vert }_2\le \tau \right\}.$$

$$L^{(p)}(\gamma ):={\left(\sum _{j=0}^m\Vert {\mathbf{x}}_{i_{j+1}}-{\mathbf{x}}_{i_j}{\Vert }^p\right)}^{1/p}$$

$$d_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }):=\min \left\{L^{(p)}(\gamma ):\gamma \text{ a path from xα to xβ through }\mathcal{X}\right\}$$

$$L^{(\infty )}(\gamma )=\underset{j=0,\dots ,m}{\max }\Vert {\mathbf{x}}_{i_{j+1}}-{\mathbf{x}}_{i_j}\Vert$$

$$d_{\mathcal{X}}^{(\infty )}({\mathbf{x}}_i,{\mathbf{x}}_j)=\min \left\{L^{(\infty )}(\gamma ):\gamma \text{ a path from xα to xβ through X}\right\}$$

$$L^{(q)}(\gamma )\le L^{(p)}(\gamma )\le n^{(1/p-1/q)}{L}^{(q)}(\gamma )$$

$$L^{(\infty )}(\gamma )\le L^{(p)}(\gamma )\le n^{1/p}{L}^{(\infty )}(\gamma )$$

$$\Vert v{\Vert }_q\le \Vert v{\Vert }_p\le m^{(1/p-1/q)}\Vert v{\Vert }_q$$

$$\Vert v{\Vert }_{\infty }\le \Vert v{\Vert }_p\le m^{(1/p)}\Vert v{\Vert }_{\infty }$$

$$v_{\gamma }={\left[\Vert {\mathbf{x}}_1-{\mathbf{x}}_0\Vert ,\Vert {\mathbf{x}}_2-{\mathbf{x}}_1\Vert ,\dots ,\Vert {\mathbf{x}}_{m+1}-{\mathbf{x}}_m\Vert \right]}^{\top }$$

$$d_{\mathcal{X}}^{(q)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le d_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le n^{(1/p-1/q)}{d}_{\mathcal{X}}^{(q)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })$$

$$\underset{\gamma }{\min }{L}^{(q)}(\gamma )\le \underset{\gamma }{\min }{L}^{(p)}(\gamma )\le n^{(1/p-1/q)}\underset{\gamma }{\min }{L}^{(q)}(\gamma )$$

$$d_{\mathcal{X}}^{(\infty )}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le d_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le n^{(1/p)}{d}_{\mathcal{X}}^{(\infty )}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }).$$

$$d_{\mathcal{X}}^{(\infty )}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le \underset{p\to \infty }{\lim }{d}_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le \underset{p\to \infty }{\lim }{n}^{(1/p)}{d}_{\mathcal{X}}^{(\infty )}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })$$

$$L^{(1)}({\gamma }_{\alpha \to \beta })=\Vert {\mathbf{x}}_{\alpha }-{\mathbf{x}}_{\beta }\Vert$$

$${ϵ}_2:=\underset{\begin{array}{c}a,b\in [\ell ]\\ a\ne b\end{array}}{\min }\underset{\begin{array}{c}{\mathbf{x}}_{\alpha }\in {\mathcal{X}}_a\\ {\mathbf{x}}_{\beta }\in {\mathcal{X}}_b\end{array}}{\min }{d}_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })$$

$$L^{(p)}(\gamma ):={\left(\sum _{j=0}^m\Vert {\mathbf{x}}_{i_{j+1}}-{\mathbf{x}}_{i_j}{\Vert }^p\right)}^{1/p}\ge {\left(\Vert {\mathbf{x}}_{i_{j^{\ast }+1}}-{\mathbf{x}}_{i_{j^{\ast }}}{\Vert }^p\right)}^{1/p}\ge \delta .$$

$$\underset{{\mathbf{x}}_{\alpha }\in {\mathcal{X}}_a,{\mathbf{x}}_{\beta }\in {\mathcal{X}}_b}{\min }{d}_{\mathcal{X}}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\ge \delta$$

$${\text{dist}}_a^{(p)}(\mathbf{x},\mathbf{y})=\underset{\eta }{\inf }{\int }_0^1\frac{\sqrt{g_a({\eta }_t^{{}^{\prime }},{\eta }_t^{{}^{\prime }})}}{{\mu }_a({\eta }_t{)}^{(1-p)/d_a}}dt$$

$$d_{{\mathcal{M}}_a}^{(p)}(\mathbf{x},\mathbf{y}):=\underset{\gamma }{\min }{\left(\sum _{j=0}^m{\text{dist}}_a({\mathbf{x}}_{i_{j+1}},{\mathbf{x}}_{i_j}{)}^p\right)}^{1/p}$$

$$\mathbb{P}\left[\underset{\begin{array}{c}\mathbf{x},\mathbf{y}\in {\mathcal{M}}_a\\ {\text{dist}}_a(\mathbf{x},\mathbf{y})\ge r_a\end{array}}{\sup }\left|\frac{{\left(d_{{\mathcal{M}}_a}^{(p)}(\mathbf{x},\mathbf{y})\right)}^p}{n_a^{(1-p)/d_a}{\text{dist}}_a^{(p)}(\mathbf{x},\mathbf{y})}-C(d_a,p)\right|>ϵ\right]=e^{-{\theta }_0^{{}^{\prime }}{n}_a^{1/p(d_a+2p)}+\mathcal{O}(\log (n_a))}$$

$$\mathbb{P}\left[\underset{{\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }\in {\mathcal{X}}_a}{\max }{d}_{{\mathcal{M}}_a}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le C_a{n}_a^{(1-p)/p{d}_a}\right]\ge 1-\exp \left(-{\theta }_0^{{}^{\prime }}{n}_a^{1/p(d_a+2p)}+\mathcal{O}(\log (n_a)\right)$$

$$d_{{\mathcal{M}}_a}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le {\left({\text{dist}}_a({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }{)}^p\right)}^{1/p}<r_a=n_a^{(1-p)/p{d}_a}$$

$$\mathbb{P}\left[\underset{{\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }\in {\mathcal{X}}_a}{\max }{d}_{{\mathcal{M}}_a}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le C_a{n}_a^{(1-p)/p{d}_a}\right]=\mathbb{P}\left[\underset{\begin{array}{c}{\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }\in {\mathcal{X}}_a\\ {\text{dist}}_a({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\ge r_a\end{array}}{\max }{d}_{{\mathcal{M}}_a}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\le C_a{n}_a^{(1-p)/p{d}_a}\right]$$

$${\left(d_{{\mathcal{M}}_a}^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })\right)}^p\le (C(d_a,p)+ϵ)n_a^{(1-p)/d_a}{\text{dist}}_a^{(p)}({\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta })$$

$$\displaystyle\Rightarrow\$$

$${\text{dist}}_a^{(p)}(\mathbf{x},\mathbf{y})$$

Full text

Keywords: clustering, shortest path distance, manifold hypothesis, unsupervised learning.

1 Introduction

Clustering high dimensional data is an increasingly important problem in contemporary unsupervised machine learning. Here, we shall consider this problem in the case where our data is presented as a subset of a Euclidean space, $\mathcal{X}\subset\mathbb{R}^{D}$, although our results easily extend to more general metric spaces. Loosely speaking, by clustering we mean partitioning $\mathcal{X}$ into $\ell$ subsets, or clusters, $\mathcal{X}=\mathcal{X}_{1}\cup\ldots\cup\mathcal{X}_{\ell}$ such that data points in the same $\mathcal{X}_{a}$ are more “similar” than data points in different subsets. Clearly, the notion of similarity is context dependent. Although there exist algorithms that operate on the data directly, for example $k$-means, many modern algorithms proceed by first representing the data as a weighted graph $G=(V,E,A)$ with $V=\{1,\ldots,n\}$ and $A_{ij}$ representing the similarity between ${\bf x}_{i}$ and ${\bf x}_{j}$ and then using a graph clustering algorithm on $G$. Spectral clustering [25] is an archetypal example of such an approach. Constructing $A$ requires a choice of metric $d(\cdot,\cdot):\mathbb{R}^{D}\times\mathbb{R}^{D}\to\mathbb{R}$. Ideally, one should choose $d$ such that points in the same cluster are close with respect to $d(\cdot,\cdot)$, while points in different clusters remain distant. Thus, the choice of metric should reflect, in some way, our assumptions about the data $\mathcal{X}$ and the notion of similarity we would like the clusters to reflect.

A common assumption, frequently referred to as the manifold hypothesis (see, for example, [15]) posits that each $\mathcal{X}_{a}$ is sampled from a latent data manifold $\mathcal{M}_{a}$. Many types of data sets are known or suspected to satisfy this hypothesis, for example motion segmentation [12, 2], images of faces or objects taken from different angles or under different lighting [4, 18] or handwritten digits [27]. It is also usually assumed that the dimension of each $\mathcal{M}_{a}$ is much lower than the ambient dimension $D$. Although it can be shown that taking $d(\cdot,\cdot)$ to be the Euclidean distance can be successful [3] for such data, data-driven metrics have been increasingly favored [10, 7, 9, 22].

Once $d(\cdot,\cdot)$ has been chosen, $A$ can be constructed. A common choice [25, 30] is to use some variant of a Gaussian kernel, whereby $A_{ij}=\exp(-d^{2}({\bf x}_{i},{\bf x}_{j})/\sigma^{2})$ for a user defined parameter $\sigma$. However this is unsuitable for large data sets, as the resulting similarity matrix is dense and thus may be too large to store in memory. Moreover, operations which form core parts of many clustering algorithms, for example computing eigenvectors, are prohibitively expensive when $A$ is full. Hence in this case, a common choice that yields a sparse similarity matrix is to use a $k$ nearest neighbors ($k$-NN) graph constructed as:

[TABLE]

Here $k$ is a user specified parameter.

In this article we consider taking $d(\cdot,\cdot)$ to be a power weighted shortest path metric (henceforth: $p$-wspm and defined in §2.2) with an emphasis on cases where the data satisfies the manifold hypothesis and where the data set is so large that a $k$-NN similarity matrix is preferable to a full similarity matrix. The use of shortest path metrics in clustering data is not new (see the discussion in §2.4), but has typically been hindered by high computational cost. Indeed finding the pairwise distance between all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$ in the shortest path metric is equivalent to the all pairs shortest paths problem on a complete weighted graph, which requires $O(n^{3})$ operations using the Floyd-Warshall algorithm. We provide a way around this computational barrier, and also contribute to the theoretical analysis of $p$-wspm’s. Specifically, our contributions are:

We prove that $p$-wspm’s behave as expected for data satisfying the manifold hypothesis. That is, we show that the maximum distance between points in the same cluster is small with high probability, and tends to zero as the number of data points tends to infinity. On the other hand, the maximum distance between points in different clusters remains bounded away from zero. 2. 2.

We show how $p$-wspm’s can be thought of as interpolants between the Euclidean metric and the longest leg path distance (defined in §2.3), which we shall abbreviate to LLPD. 3. 3.

We introduce a novel modified version of Dijkstra’s algorithm that computes the $k$ nearest neighbors, with respect to any $p$-wspm or the LLPD, of any ${\bf x}_{\alpha}$ in $\mathcal{X}$ in $O(k^{2}\mathcal{T}_{Enn})$ time, where $\mathcal{T}_{Enn}$ is the cost of a Euclidean nearest-neighbor query. Hence one can construct a $p$-wspm k-NN graph in $O(nk^{2}\mathcal{T}_{Enn})$. As we typically have $k\ll n$, i.e. $k=O(\log(n))$ or even $k=O(1)$, this means that constructing a $p$-wspm k-NN graph requires only marginally more time than constructing a Euclidean $k$-NN graph (which requires $O(nk\mathcal{T}_{Enn})$). 4. 4.

We verify experimentally that using a $p$-wspm in lieu of the Euclidean metric results in an appreciable increase in clustering accuracy, at the cost of a small increase in run time, for a wide range of real and synthetic data sets.

After establishing notation and surveying the literature in §2, we prove our main results in §3 and §4. In §5 we present our algorithm for computing $k$ nearest neighbors in any $p$-wspm, while in §6 we report the results of our numerical experiments.

2 Definitions and Notation

Let us first fix some notation. Throughout this paper, $D$ shall denote the ambient dimension, while $\mathcal{X}$ will denote a fixed, finite sets of distinct points in $\mathbb{R}^{D}$. We shall denote the Euclidean (i.e. $\ell_{2}$) norm on $\mathbb{R}^{D}$ as $\|\cdot\|$. For any finite set $S$, by $|S|$ we shall mean its cardinality. For any positive integer $\ell$, by $[\ell]$ we mean the set $\{1,2,\ldots,\ell\}$. Finally, for two functions $f(n)$ and $g(n)$ by $f(n)=O(g(n))$ we shall mean that there exist constants $C$ and $n_{0}$ such that $f(n)\leq Cg(n)$ for all $n\geq n_{0}$. Similarly, by $f(n)=o(g(n))$ we shall mean that $f(n)/g(n)\to 0$ as $n\to\infty$. Occasionally we shall explicitly indicate the dependence on the variable $n$ (later, also $n_{a}$) by writing $O_{n}$ or $o_{n}$.

2.1 Data Model

We consider data sets $\mathcal{X}=\mathcal{X}_{1}\cup\mathcal{X}_{2}\cup\ldots\cup\mathcal{X}_{\ell}\subset\mathbb{R}^{D}$ consisting naturally of $\ell$ clusters, which are a priori unknown. Let $|\mathcal{X}_{a}|=n_{a}$. We posit that for each $\mathcal{X}_{a}$ there is a smooth, compact, embedded manifold $\mathcal{M}_{a}\hookrightarrow\mathbb{R}^{D}$ such that $\mathcal{X}_{a}\subset\mathcal{M}_{a}$. Let $g_{a}$ denote the restriction of the Euclidean metric to $\mathcal{M}_{a}$, then $(\mathcal{M}_{a},g_{a})$ is a compact Riemannian manifold. We shall further assume that $\mathcal{X}_{a}$ is sampled according to a probability density function $\mu_{a}$ supported on $\mathcal{M}_{a}$ and continuous with respect to the metric $g_{a}$. For any ${\bf x},{\bf y}\in\mathcal{M}_{a}$ let

[TABLE]

denote the metric induced by $g_{a}$, where the infimum is over all piecewise smooth curves $\gamma:[0,1]\to\mathcal{M}_{a}$ with $\gamma(0)={\bf x}$ and $\gamma(1)={\bf y}$. Define the diameter of $\mathcal{M}_{a}$ to be the supremum over all distances between points in $\mathcal{M}_{a}$:

[TABLE]

Since each $\mathcal{M}_{a}$ is compact this supremum is in fact a maximum and $\text{diam}(\mathcal{M}_{a})$ is finite. We assume that the data manifolds are fairly well separated, that is,

[TABLE]

Note that frequently (for example, in [3]), this model is extended to allow for noisy sampling, whereby for some $\tau>0$, $\mathcal{X}_{a}$ is sampled from the tube $B(\mathcal{M}_{a},\tau)$ defined as:

[TABLE]

but we leave this extension to future work.

2.2 Power Weighted Shortest Path Metrics

For any distinct pair ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$ and any path $\gamma={\bf x}_{\alpha}\to{\bf x}_{1}\to\ldots\to{\bf x}_{m}\to{\bf x}_{\beta}$, define the $p$-weighted length of $\gamma$ to be:

[TABLE]

where by convention we declare ${\bf x}_{i_{0}}={\bf x}_{\alpha}$ and ${\bf x}_{i_{m+1}}={\bf x}_{\beta}$. We define the $p$-weighted shortest path distance from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}$ to be the minimum length over all such paths:

[TABLE]

Note that $d^{(p)}_{\mathcal{X}}$ is a metric on the set $\mathcal{X}$ for $p\geq 1$ (see, for example, [19]). As several authors [20, 9, 1] have noted, the metric $d^{(p)}_{\mathcal{X}}$ is density-dependent, so that if ${\bf x}_{\alpha}$ and ${\bf x}_{\beta}$ are contained in a region of high density (i.e. a cluster) the path distance $d^{(p)}({\bf x}_{\alpha},{\bf x}_{\beta})$ will likely be shorter than the Euclidean distance $\|{\bf x}_{\alpha}-{\bf x}_{\beta}\|$ (as long as $p>1$).

2.3 Longest-Leg Path Distance

Another common path-based distance is the longest-leg path distance (LLPD), which we shall denote as $d^{(\infty)}_{\mathcal{X}}$ (the choice of this notation should become clear shortly). It is defined as the minimum, over all paths from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}$, of the maximum distance between consecutive points in the path (i.e. legs). Before formally defining $d^{(\infty)}_{\mathcal{X}}$, define, for any path $\gamma$ from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}$, the longest-leg length of $\gamma$ as:

[TABLE]

again we are using the convention that ${\bf x}_{i_{0}}={\bf x}_{\alpha}$ and ${\bf x}_{i_{m+1}}={\bf x}_{\beta}$. Now, in analogy with (3):

[TABLE]

$d^{\infty}_{\mathcal{X}}$ is also a metric, in fact an ultrametric [22], on $\mathcal{X}$.

2.4 Prior work

The idea of using $p$-wspm’s for clustering was proposed in [28], and further explored in [26]. Recently, several papers [16, 8, 22] have considered the use of LLPD for clustering and in particular [22] provides performance guarantees for spectral clustering with LLPD for a data model that is similar to ours. [9] studies $p$-wspm’s for $p\geq 2$ and proposes to use such metrics with density-based clustering algorithms, such as DBScan, although they do not provide any experimental results. The paper [7] proposes the use of $p$-wspm’s for semi-supervised learning and provides a fast Dijkstra-style algorithm for finding, for every ${\bf x}\in\mathcal{X}$, its nearest neighbor, with respect to a $p$-wspm, in some set of labeled data points $\mathcal{L}$. They consider a similar data model to ours, but do not provide any quantitative results on the behaviour of shortest path distances. More generally, shortest path distances are a core part of the ISOMAP dimension reduction algorithm [27], although we emphasize that here not all paths through $\mathcal{X}$ are considered—first a $k$-NN graph $G^{(k)}$ is computed from $\mathcal{X}$ and only paths in this graph are admissible.

The asympotic behaviour of power weighted shortest path distances are analyzed for Euclidean Poisson processes in [19] and for points sampled from an arbitrary probability distribution supported on a Riemannian manifold $\mathcal{M}$ in [20]. Note that in [20] the lengths of the legs of the path are measured using geodesic distance on $\mathcal{M}$, which is not computable in the data model we are considering as the $\mathcal{M}_{a}$ are unknown. Finally, in [9] the results of [20] are used to show that, for certain $k$ and $p\geq 2$, with high probability the Euclidean $k$-NN graph can be used to compute $p$-wspm distances in the case where the data is sampled from a single Riemannian manifold $\mathcal{M}$. We discuss this further in §5.3.

On the computational side, we are unaware of any prior mention of Algorithm 2 in the literature, although similar algorithms, which solve slightly different problems, are presented in [17], [23] and [7]. In [7] a variation of Dijkstra’s algorithm is presented which starts from a set of labeled data $\mathcal{L}$ and expands its search through the rest of the data $\mathcal{X}$. The algorithm is designed to find, for all ${\bf x}_{i}\in\mathcal{X}$, the nearest point in $\mathcal{L}$ to ${\bf x}_{i}$ with respect to a $p$-wspm and then terminate. Thus, it is not clear how one would extend this algorithm to finding $k$ nearest neighbors for $k>1$. The algorithm of [17] is formulated for any weighted graph $G=(V,E,A)$ (i.e. not just graphs obtained from data sets $\mathcal{X}\subset\mathbb{R}^{D}$) and as such is not well-adapted to the problem at hand. In particular, it has run time $O(k(n\log(n)+|E|))$. Because the distance graph obtained from $\mathcal{X}$ is implicitly complete, $|E|=O(n^{2})$ and this results in a run time proportional to $kn^{2}$, which is infeasible for large data sets. Finally, the algorithm presented in [23], although adapted to the situation of distance graphs of data sets, actually solves a slightly different problem. Specifically they consider finding the $k_{1}$ $p$-wspm nearest neighbors of each ${\bf x}\in\mathcal{X}$ in a $k_{2}$ Euclidean nearest neighbors graph of $\mathcal{X}$. As such, it is not clear whether the set of nearest neighbors produced by their algorithm are truly the $p$-wspm nearest neighbors in $\mathcal{X}$.

Let us also mention that our approach is “one at a time”, whereas the other three algorithms mentioned are “all at once”. That is, our algorithm takes as input ${\bf x}\in\mathcal{X}$ and outputs the $k$ $p$-wpsm nearest neighbors of ${\bf x}$. This can then be iterated to find the $p$-wspm nearest neighbors of all ${\bf x}\in\mathcal{X}$. In contrast, “all at once” algorithms directly return the sets of $k$ nearest neighbors for each ${\bf x}\in\mathcal{X}$. Thus it is possible our algorithm will have applications in other scenarios where the $p$-wspm nearest neighbors of only some small subset of points of $\mathcal{X}$ are required or in “online” scenarios where new data points are continuously received.

3 Relation between $p$-wspm’s for different values of $p$

Here we compare $p$-wspm’s for different values of the power weighting, $p$, in the interval $[1,\infty]$.

Theorem 3.1**.**

For any fixed ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$, and any fixed path $\gamma$ from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$:

If $1\leq p<q<\infty$ then:

[TABLE] 2. 2.

For any $1\leq p<\infty$:

[TABLE]

where $n:=|\mathcal{X}|$.

Proof.

We shall make use of the following well-known results from analysis:

for any fixed $m$, any $v\in\mathbb{R}^{m}$ and any $0<p<q$ we have that

[TABLE] 2. 2.

for any fixed $m$, any $v\in\mathbb{R}^{m}$ and any $0<p<\infty$:

[TABLE]

For any path $\gamma={\bf x}_{\alpha}\to{\bf x}_{1}\to\ldots\to{\bf x}_{m}\to{\bf x}_{\beta}$, observe that $L^{(p)}(\gamma)$ can be thought of as $\|v_{\gamma}\|_{p}$, where $v_{\gamma}\in\mathbb{R}^{m+1}$ is the vector:

[TABLE]

It follows from (5) that $L^{(q)}(\gamma)\leq L^{(p)}(\gamma)\leq(m+1)^{(1/p-1/q)}L^{(q)}(\gamma)$. Observe that $m+1$, the path length, is always less than or equal to $n$, the number of points in $\mathcal{X}$. Thus $(m+1)^{(1/p-1/q)}\leq n^{(1/p-1/q)}$ and part 1. of the theorem follows. In a similar fashion, we obtain part 2. from (6). ∎

Corollary 3.2**.**

For any fixed $\mathcal{X}$, any $1\leq p<q<\infty$ and all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$:

[TABLE]

where $n:=|\mathcal{X}|$.

Proof.

By Theorem 3.1: $L^{(q)}(\gamma)\leq L^{(p)}(\gamma)\leq n^{(1/p-1/q)}L^{(q)}(\gamma)$ for all paths $\gamma$ between ${\bf x}_{\alpha}$ and ${\bf x}_{\beta}$, hence the same relationship must hold for the minimum over all paths from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$:

[TABLE]

which, appealing to the definition of $d_{\mathcal{X}}^{(p)}({\bf x}_{\alpha},{\bf x}_{\beta})$ (see (3)) yields the corollary. ∎

Corollary 3.3**.**

For any fixed $\mathcal{X}$, we have that $\displaystyle\lim_{p\to\infty}d^{(p)}({\bf x}_{\alpha},{\bf x}_{\beta})=d^{(\infty)}({\bf x}_{\alpha},{\bf x}_{\beta})$ for all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$.

Proof.

Again by Theorem 3.1 we have that $L^{(\infty)}(\gamma)\leq L^{(p)}(\gamma)\leq n^{(1/p)}L^{(\infty)}(\gamma)$ for all paths $\gamma$ between ${\bf x}_{\alpha}$ and ${\bf x}_{\beta}$, hence by the same reasoning as in the proof of Corollary 3.2:

[TABLE]

Because this holds for all $p<\infty$, taking the limit we get:

[TABLE]

and the result follows from the fact that $\lim_{p\to\infty}n^{(1/p)}=1$. ∎

Theorem 3.4**.**

For all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$, $d^{(1)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})=\|{\bf x}_{\alpha}-{\bf x}_{\beta}\|$

Proof.

$d^{(1)}$ is defined as a minimum over all paths from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}$, and in particular the one hop path $\gamma_{\alpha\to\beta}={\bf x}_{\alpha}\to{\bf x}_{\beta}$ is such a path. We claim it is the shortest such path as for any other path $\gamma={\bf x}_{\alpha}\to{\bf x}_{i_{1}}\to\ldots\to{\bf x}_{i_{m}}\to{\bf x}_{\beta}$ by repeated applications of the triangle inequality:

[TABLE]

∎

4 $p$-wspm’s in the Multi-Manifold Setting

One of the most useful aspects of $p$-wspm’s, when applied to clustering problems, is that they tend to “squeeze” points in the same cluster together, while (hopefully) keeping points in different clusters separated. Here we make this more precise. Specifically we show that for any $1<p<\infty$ if the data comes from the model described in §2.1 then:

• •

$\displaystyle\min_{{\bf x}_{\alpha}\in\mathcal{X}_{a},{\bf x}_{\beta}\in\mathcal{X}_{b}}d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})\geq\delta>0$ (see Lemma 4.1). Recall that $\delta$ is the minimal separation between data manifolds.

• •

$\displaystyle\max_{a\in[k]}\max_{{\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}_{a}}d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})=O(n_{\min}^{(1-p)/pd_{\max}})$ with probability tending to $1$ as $n$ tends to $\infty$. (see Theorem 4.6).

In this section it is sometimes necessary to enlarge our definition of $p$-wspm to allow for paths between ${\bf x},{\bf y}\in\mathbb{R}^{D}$ that are not necessarily in $\mathcal{X}$ (and points that are not in $\mathcal{X}$ shall be denoted without a subscript). Thus $d^{(p)}_{\mathcal{X}}({\bf x},{\bf y})$ is technically defined as, using the notation of §2.2, $d^{(p)}_{\mathcal{X}\cup\{{\bf x},{\bf y}\}}({\bf x},{\bf y})$.

4.1 Paths between points in different clusters

Here we prove that $p$-wspm’s maintain a separation between points in different clusters.

Theorem 4.1**.**

Let $\epsilon_{2}$ denote the minimal distance between points in different clusters. That is:

[TABLE]

Then $\epsilon_{2}\geq\delta$ with $\delta$ as defined in (1).

Proof.

For any ${\bf x}_{\alpha}\in\mathcal{X}_{a}$ and ${\bf x}_{\beta}\in\mathcal{X}_{b}$ let $\gamma={\bf x}_{\alpha}\to{\bf x}_{i_{1}}\to,\ldots,\to{\bf x}_{i_{m}}\to{\bf x}_{\beta}$ be any path from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}$, where again we are using the convention that ${\bf x}_{i_{0}}:={\bf x}_{\alpha}$ and ${\bf x}_{i_{m+1}}={\bf x}_{\beta}$. If ${\bf x}_{\alpha}\in\mathcal{X}_{a}$ and ${\bf x}_{\beta}\in\mathcal{X}_{b}$ there must exist (at least one) $j^{*}\in[m]$ such that ${\bf x}_{i_{j^{*}}}\in\mathcal{X}_{a}$ while ${\bf x}_{i_{j^{*}+1}}\in\mathcal{X}_{b}$. By the assumptions on the generative model, $\mathcal{X}_{a}\subset\mathcal{M}_{a}$ and $\mathcal{X}_{b}\subset\mathcal{M}_{b}$ and so: $\|{\bf x}_{i_{j^{*}+1}}-{\bf x}_{i_{j^{*}}}\|^{p}\geq\left(\text{dist}(\mathcal{M}_{a},\mathcal{M}_{b})\right)^{p}=\delta^{p}$ thus:

[TABLE]

Because this holds for all such $\gamma$ we have $d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta}):=\min_{\gamma}\left\{L^{(p)}(\gamma)\right\}\geq\delta$ and because this holds for all such ${\bf x}_{\alpha}$ and ${\bf x}_{\beta}$:

[TABLE]

Finally, this holds for all $a\neq b$, yielding the lemma. ∎

4.2 Asymptotic Limits of power weighted shortest paths

For all $a\in[\ell]$, define $d^{(p)}_{\mathcal{X}_{a}}({\bf x}_{\alpha},{\bf x}_{\beta})$ as the minimum $p$-weighted length of paths from ${\bf x}_{\alpha}$ to ${\bf x}_{\beta}$ through $\mathcal{X}_{a}$ (i.e. we are excluding paths that may pass through points in $\mathcal{X}\setminus\mathcal{X}_{a}$). Because $\mathcal{X}_{a}\subset\mathcal{X}$, it follows that $d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})\leq d^{(p)}_{\mathcal{X}_{a}}({\bf x}_{\alpha},{\bf x}_{\beta})$111More generally the reader is invited to check that for any $\mathcal{Y}\subset\mathcal{X}$ we have that $d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})\leq d^{(p)}_{\mathcal{Y}}({\bf x}_{\alpha},{\bf x}_{\beta})$.. In this section we address the asymptotic behaviour of $d^{(p)}_{\mathcal{X}_{a}}({\bf x}_{\alpha},{\bf x}_{\beta})$. Here is where we make critical use of the main theorem of [20], which we state as Theorem 4.2. Recall that $\mu_{a}$ is the probability density function with respect to which $\mathcal{X}_{a}$ is sampled from $\mathcal{M}_{a}$, and that by assumption $\mu^{\min}_{a}:=\displaystyle\min_{x\in\mathcal{M}_{a}}\mu_{a}(x)>0$. Define the following power-weighted geodesic distance on $\mathcal{M}_{a}$:

[TABLE]

where here the infimum is over all piecewise smooth paths $\eta:[0,1]\to\mathcal{M}_{a}$ with $\eta(0)={\bf x}$ and $\eta(1)={\bf y}$. As in §2.1, for the Riemannian manifold $(\mathcal{M}_{a},g_{a})$ let $\text{dist}_{a}({\bf x},{\bf y})$ denotes the geodesic distance from ${\bf x}$ to ${\bf y}$ on $\mathcal{M}_{a}$ with respect to $g_{a}$.

In order to bound $d^{(p)}_{\mathcal{X}_{a}}({\bf x}_{\alpha},{\bf x}_{\beta})$ we study an auxiliary shortest path distance, $d^{(p)}_{\mathcal{M}_{a}}$. This distance will also be defined as a minimum over $p$-weighted path lengths, but instead of measuring the length of the legs using the Euclidean distance $\|\cdot\|$, we measure them with respect to the geodesic distance $\text{dist}_{a}$:

[TABLE]

where again the $\min$ is over all paths $\gamma$ from ${\bf x}$ to ${\bf y}$ through $\mathcal{X}_{a}$.

Theorem 4.2**.**

Let $\mathcal{M}_{a}$ be a compact Riemannian manifold, and assume that $\mathcal{X}_{a}$ is drawn from $\mathcal{M}_{a}$ with continuous probability distribution $\mu_{a}$ satisfying $\min_{x\in\mathcal{M}_{a}}\mu_{a}(x)>0$. Let $n_{a}:=|\mathcal{X}_{a}|$. For all $n_{a}$, let $r_{a}:=n_{a}^{(1-p)/pd_{a}}$. Then for any $1\leq p<\infty$ and any fixed $\epsilon>0$ there exists a constant $\theta^{{}^{\prime}}_{0}$ independent of $n_{a}$ such that:

[TABLE]

where $C(d_{a},p)$ is a constant depending only on $d_{a}$ and $p$, but not on $n_{a}$.

Remark 4.3**.**

This is a slightly modified version Theorem 1 in [20]. As stated in [20], the $\sup$ is over ${\bf x},{\bf y}\in\mathcal{M}$ satisfying $\text{dist}_{a}({\bf x},{\bf y})\geq b$ for a fixed constant $b$. However immediately below the statement of Theorem 1 the authors acknowledge that one can weaken this assumption to $\text{dist}_{a}({\bf x},{\bf y})\geq r_{a}$ as long as $n_{a}r_{a}^{d_{a}}/\log n_{a}\to\infty$, which is the case for our choice of $r_{a}$. Note that there is a slight notational discrepancy here. What is called $d^{(p)}_{\mathcal{X}}({\bf x}_{\alpha},{\bf x}_{\beta})$ in [20] is our $\left(d^{(p)}_{\mathcal{M}_{a}}({\bf x}_{\alpha},{\bf x}_{\beta})\right)^{p}$.

Corollary 4.4**.**

With assumptions as in Theorem 4.2, there exists a constant $\theta^{{}^{\prime}}_{0}$, independent of $n_{a}$ such that:

[TABLE]

where $C_{a}$ is a constant depending on $d_{a},p,\mu_{a}^{\min}$ and $\text{diam}(\mathcal{M}_{a})$ but not on $n_{a}$.

Proof.

First observe that because the because the one leg path $\gamma_{\alpha\to\beta}={\bf x}_{\alpha}\to{\bf x}_{\beta}$ is trivially a path through $\mathcal{X}_{a}$, for any ${\bf x}_{\alpha},{\bf x}_{\beta}$ satisfying $\text{dist}_{a}({\bf x}_{\alpha},{\bf x}_{\beta})<r_{a}$ we have that:

[TABLE]

Hence as long as $C_{a}\geq 1$ we get that:

[TABLE]

Because $\mathcal{X}_{a}\subset\mathcal{M}_{a}$ we may use Theorem 4.2 to bound this probability. Indeed, fix any small $\epsilon<1$. Then there exists a $\theta_{0}^{{}^{\prime}}$ such that with probability at least $1-\exp\left(-\theta_{0}^{{}^{\prime}}n_{a}^{1/p(d_{a}+2p)}+\mathcal{O}(\log(n_{a})\right)$ we have:

[TABLE]

for all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}_{a}$ satisfying $\text{dist}_{a}({\bf x}_{\alpha},{\bf x}_{\beta})\geq r_{a}$. We now upper-bound the bracketed quantity in (10). From the definition of $\text{dist}_{a}^{(p)}$ (see (7))

[TABLE]

Because $\mathcal{M}_{a}$ is compact and embedded, its diameter (see §2.1) is finite, and $\text{dist}_{a}({\bf x},{\bf y})\leq\text{diam}(\mathcal{M}_{a})$. So for all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}$ with $\text{dist}_{a}({\bf x}_{\alpha},{\bf x}_{\beta})\geq r_{a}$:

[TABLE]

Defining $C_{a}=\max\{\tilde{C}_{a},1\}$ we get that

[TABLE]

thus proving the corollary. ∎

Finally, it remains to compare the path distance with Euclidean legs, $d^{(p)}_{\mathcal{X}_{a}}$, to the path distance with geodesic legs, $d^{(p)}_{\mathcal{M}_{a}}$.

Lemma 4.5**.**

For any ${\bf x},{\bf y}\in\mathcal{M}_{a}$, and for all $a\in[k]$, $d^{(p)}_{\mathcal{X}_{a}}({\bf x},{\bf y})\leq d^{(p)}_{\mathcal{M}_{a}}({\bf x},{\bf y})$

Proof.

Observe that for any ${\bf x},{\bf y}\in\mathcal{M}_{a}$, $\|{\bf x}-{\bf y}\|\leq\text{dist}_{a}({\bf x},{\bf y})$. It follows that for any path $\gamma={\bf x}\to{\bf x}_{i_{1}}\to\ldots\to{\bf x}_{i_{m}}\to{\bf y}$ through $\mathcal{X}_{a}$:

[TABLE]

and so:

[TABLE]

whence the result follows. ∎

4.3 Paths Between Points in the Same Cluster

Let us now return to the full distance function $d^{(p)}_{\mathcal{X}}$.

Theorem 4.6**.**

Define $\epsilon_{1}$ to be the maximal distance between points in the same cluster:

[TABLE]

With assumptions as in §2.1, for any $1\leq p<\infty$:

[TABLE]

Proof.

First, for all ${\bf x}_{\alpha},{\bf x}_{\beta}\in\mathcal{X}_{a}$ observe that:

[TABLE]

where the first inequality is because $\mathcal{X}_{a}\subset\mathcal{X}$ and the second is Lemma 4.5. Now let $C:=\max_{a}C_{a}$. Clearly $C_{a}n_{a}^{(1-p)/pd_{a}}\leq Cn^{(1-p)/pd_{\max}}$ and similarly:

[TABLE]

combining these observations, (12) and Corollary 4.4:

[TABLE]

for all $a=1,\ldots,k$. By the union bound, and the definition of $\epsilon_{1}$:

[TABLE]

∎

4.4 Discussion

Theorem 4.6 reveals an interesting tradeoff, already present in the work of [19, 20], namely that by increasing $p$ we get a tighter upper bound on $\epsilon_{1}$, but it holds with lower probability. We find experimental evidence for this in §6. In [3] and [22] bounds analogous to those in Theorems 4.1 and 4.6, but for $\|\cdot\|$ and $d^{(\infty)}_{\mathcal{X}}$ respectively, are used to provide bounds on the performance of single-linkage heirarchical clustering and spectral clustering with a full similarity matrix. As the focus of this article is clustering with a $k$ nearest neighbors similarity matrix, we do not pursue this line of inquiry further here.

5 A Fast Algorithm for $p$-wspm Nearest Neighbors

In this section we start from a more general perspective. Let $G=(V,E,A)$ be a weighted graph with weighted adjacency matrix $A$. Occasionally we shall find it more convenient not to fix an ordering of the vertices, in which case $A(u,v)$ will represent the weight of the edge $\{u,v\}$. We assume all edge weights are positive. For any $v\in V$ we denote its set of neighbors by $\mathcal{N}(v)$. By $\gamma=u\to w_{1}\to\ldots\to w_{m}\to v$ we shall mean the path from $u$ to $v$ in $G$ through $w_{1},\ldots,w_{m}$. Here, this is only valid if $\{u,w_{1}\},\ldots,\{w_{i},w_{i+1}\},\ldots,\{w_{m},v\}$ are all edges in $G$. In analogy with §2.2 we maintain the convention that for such a path $\gamma$, $w_{0}=u$ and $w_{m+1}=v$. Define the length of $\gamma$ as the sum of all its edge weights: $L(\gamma):=\sum_{i=0}^{m}A(w_{i},w_{i+1})$ and similarly define the longest-leg length of $\gamma$ as: $L^{(\infty)}(\gamma)=\max_{i=0}^{m}A(w_{i},w_{i+1})$. For any $u,v\in V$ define the shortest path distance as:

[TABLE]

and analogously define the longest-leg path distance as:

[TABLE]

Definition 5.1**.**

Let $\mathcal{N}_{k,G}(v)$ denote the set of $k$ nearest neighbors of $v\in V$. That is, $\mathcal{N}_{k,G}(v)=\{w_{1},\ldots,w_{k}\}$ with $A(v,w_{1})\leq A(v,w_{2})\leq\ldots\leq A(v,w_{k})\leq A(v,w)$ for all $w\in V\setminus\{w_{1},\ldots,w_{k}\}$. By convention, we take $v\in\mathcal{N}_{k,G}(v)$

Definition 5.2**.**

For any graph $G$, define a directed $k$ nearest neighbors graph $G^{(k)}$ with directed edges $(u,v)$ whenever $v\in\mathcal{N}_{k,G}(u)$.

In practice we do not compute the entire edge set of $G^{(k)}$, but rather just compute the sets $\mathcal{N}_{k,G}(u)$ as it becomes necessary.

Definition 5.3**.**

Let $\mathcal{N}^{d_{G}}_{k,G}(v)$ denote the set of $k$ vertices which are closest to $v$ in the shortest-path distance $d_{G}$. That is, $\mathcal{N}^{d_{G}}_{k,G}(v)=\{w_{1},\ldots,w_{k}\}$ and $d_{G}(v,w_{1})\leq d_{G}(v,w_{2})\leq\ldots d_{G}(v,w_{k})\leq d_{G}(v,w)$ for all $w\in V\setminus\mathcal{N}^{d_{G}}_{k,G}(v)$. By convention, we take $v$ to be in $\mathcal{N}^{d_{G}}_{k,G}(v)$. Similarly, we define $\mathcal{N}^{(\infty)}_{k,G}(v)$ to be the $k$ vertices closest to $v$ in the metric $d^{(\infty)}_{G}$.

We have not specified how to break ties in the definition of $\mathcal{N}_{k,G}(v),\mathcal{N}^{d_{G}}_{k,G}(v)$ or $\mathcal{N}^{(\infty)}_{k,G}(v)$. For the results of this section to hold, any method will suffice, as long as we use the same method in all three cases. To simplify the exposition, we shall assume henceforth that all distances are distinct.

Remark 5.4**.**

Let us relate this to the discussion in previous sections. For any set of data points $\mathcal{X}=\{{\bf x}_{1},\ldots,{\bf x}_{n}\}\subset\mathbb{R}^{D}$ and any power weighting $1\leq p<\infty$ one can form a graph $G$ on $n$ vertices, where the vertex $v_{i}$ corresponds to the data point ${\bf x}_{i}$, and edge weights $A_{ij}=\|{\bf x}_{i}-{\bf x}_{j}\|^{p}$. Then:

[TABLE]

Moreover, if we denote by $\mathcal{N}^{(p)}_{k,\mathcal{X}}({\bf x}_{i})$ the $k$ nearest neighbors, with respect to $d^{(p)}_{\mathcal{X}}$, of ${\bf x}_{i}\in\mathcal{X}$, then:

[TABLE]

The analogous results also hold for $d_{\mathcal{X}}^{(\infty)}$, ie

[TABLE]

Before proceeding, let us briefly review how Dijkstra’s algorithm works. For a graph $G$ with non-negative edge weights, and a given source vertex, $s$, Dijkstra will return a list of pairs of the form $(u,d_{G}(s,u))$. The computational complexity of Dijkstra’s algorithm is $\mathcal{O}(|E|+n\log n)$ so when $G$ is a complete graph the complexity is $\mathcal{O}(n^{2})$. Intuitively, our proposed algorithm (Algorithm 2) sidesteps this complexity by allowing one to run Dijkstra’s algorithm on a much sparser graph, as long as one is only interested in determining the identity and path distance to the $k$ nearest neighbors of $s$, with respect to the path distance $d_{G}$ (see Lemma 5.6).

The following implementation of Dijkstra’s algorithm is as in [11]. The min-priority queue operations $\text{\tt decreaseKey},\text{\tt insert}$ and extractMin have their standard definitions (see, for example Chpt. 6 of [11]). For any vertex $s\in V$ and any subset $W\subset V$, we shall also use the shorthand $\text{\tt makeQueue}(W,s)$ to denote the process of initializing a min-priority queue with $\text{\tt key}[s]=0$ and $\text{\tt key}[v]=+\infty$ for all $v\in W\setminus s$.

Note once $u$ is popped in step $4$, $\text{\tt key}[u]=d_{G}(s,u)$. Our first key observation is the following:

Lemma 5.5**.**

Suppose that all weights are non-negative: $A_{ij}\geq 0$. If $u_{i}$ is the $i$-th vertex to be removed from $Q$ at step $11$, then $u_{i}$ is the $i$-th closest vertex to $s$.

Proof.

See, for example, the discussion in [11]. ∎

It follows that, if one is only interested in finding the $k$ nearest neighbors of $s$ in the path distance $d_{G}$, one need only iterate through the “while” loop $3\to 12$ $k$ times. There is a further inefficiency, which was also highlighted in [7]. The “for” loop 6–10 iterates over all neighbors of $u$. The graphs we are interested in are, implicitly, fully connected, hence this for loop iterates over all $n-1$ other vertices at each step. We fix this with the following observation:

Lemma 5.6**.**

For any graph $G$, let $G^{(k)}$ denote its $k$-Nearest-Neighbor graph (see Definition 5.2). Then:

[TABLE]

Note that in the directed graph $G^{(k)}$, we consider only paths that traverse each edge in the ‘correct’ direction.

Concretely: the path-nearest-neighbors in $G$ are the same as the path-nearest neighbors in $G^{(k)}$, hence one can find $\mathcal{N}^{d_{G}}_{k,G}(v)$ by running a Dijkstra-style algorithm on $G^{(k)}$, instead of $G$. As each vertex in $G^{(k)}$ has a small number of neighbors (precisely $k$), this alleviates the computational burden of the “for” loop 6–10 highlighted above.

Before proving this lemma, let us explain why it may seem counterintuitive. If $w\in\mathcal{N}^{d_{G}}_{k,G}(v)$ there is a path $\gamma$ from $v$ to $w$ that is short (at least shorter than the shortest paths to all $u\notin\mathcal{N}^{d_{G}}_{k,G}(v)$). In forming $G^{(k)}$ from $G$, one deletes a lot of edges. Thus it is not clear that $\gamma$ is still a path in $G^{(k)}$ (some of its edges may now be “missing”). Hence it would seem possible that $w$ is now far away from $v$ in the shortest-path distance in $G^{(k)}$. The lemma asserts that this cannot be the case.

Proof.

Since the sets $\mathcal{N}^{d_{G}}_{k,G}(v)$ and $\mathcal{N}^{d_{G^{(k)}}}_{k,G^{(k)}}(v)$ have the same cardinality (i.e. $k$), to prove equality it suffices to prove one containment. We shall show that $\mathcal{N}^{d_{G}}_{k,G}(v)\subset\mathcal{N}^{d_{G^{(k)}}}_{k,G^{(k)}}(v)$. Consider any $w\in\mathcal{N}^{d_{G}}_{k,G}(v)$. Let $\tilde{\gamma}=v\to u_{1}\to\ldots\to u_{m}\to w$ be a shortest path from $v$ to $w$. That is, L(\tilde{\gamma})=\min\{L(\gamma):\gamma\text{ a path from v$to$w}\}.

We claim that $\tilde{\gamma}$ is a path in $G^{(k)}$. If this not the case, then there is an edge $\{u_{i},u_{i+1}\}$ that is in $\tilde{\gamma}$ but $(u_{i},u_{i+1})$ is not an edge in $G^{(k)}$ (we again adopt the convention that $u_{0}:=v$ and $u_{m+1}:=w$). By the construction of $G^{(k)}$ this implies that there are $k$ vertices $x_{1},\ldots,x_{k}$ that are closer to $u_{i}$ than $u_{i+1}$. (Note that the sets $\{u_{0},\ldots,u_{i-1}\}$ and $\{x_{1},\ldots,x_{k}\}$ need not be disjoint). But then the paths $\gamma_{j}=v\to u_{1}\to\ldots\to u_{i}\to x_{j}$ in $G$ are all shorter than the path $v\to u_{1}\to\ldots\to u_{i+1}$ and hence shorter than $\tilde{\gamma}$, as all edge weights are assumed positive. It follows that $d_{G}(v,x_{j})<d_{G}(v,w)$ for $j=1,\ldots,k$, contradicting the assumption that $w\in\mathcal{N}^{d_{G}}_{k,G}(v)$.

*Now, we claim that * $w\in\mathcal{N}^{d_{G^{(k)}}}_{k,G^{(k)}}(v)$. If this were not the case, there would exists $k$ other vertices $w_{1},\ldots,w_{k}$ that are closer in the shortest-path distance $d_{G^{(k)}}$ to $v$ than $w$. That is, there would be paths $\gamma_{1},\ldots,\gamma_{k}$ from $v$ to $w_{1},\ldots,w_{k}$ respectively that are shorter than $\gamma$. But every path in $G^{(k)}$ is also a path in $G$, hence $w_{1},\ldots,w_{k}$ are also closer to $v$ than $w$ in the shortest-path distance $d_{G}$. This contradicts the assumption that $w\in\mathcal{N}^{d_{G}}_{k,G}(v)$.∎

There is a final, minor, inefficiency in Algorithm 1 that we can improve upon; $Q$ is initialized to contain all vertices $V$ when it is actually only necessary to initialize it to contain the neighbors of $s$. Combining these three insights we arrive at Algorithm 2. We call this algorithm Dijkstra-with-pruning as the key idea, expressed in Lemma 5.6, is to “prune” the neighborhood of each $v\in V$ down to the $k$ nearest neighbors of $v$. Note that we use DecreaseOrInsert as shorthand for the function that decreases $\text{\tt key}[v]$ to tempDist if $\text{tempDist}<\text{\tt key}[v]$ and $v\in Q$, inserts $v$ into $Q$ with priority $\text{\tt key}[v]=\text{tempDist}$ if $v\notin Q$ and does nothing if $v\in Q$ but $\text{tempDist}\geq\text{\tt key}[v]$. In fact, this is equivalent to inserting a copy of $v$ into $Q$ with priority $\text{\tt key}[v]=\text{tempDist}$, hence DecreaseOrInsert has the same computational complexity as insert (see also [23]). Note that in this implementation the size of $Q$ grows by one on every iteration of the inner for loop, 10–13.

Theorem 5.7**.**

For any $s$ and any $G$ with positive weights, Algorithm 2 is correct. That is, $S$ contains precisely the pairs $(v,d_{G}(s,v))$ for all $v\in\mathcal{N}^{d_{G}}_{G,k}(s)$.

Proof.

By only using $\mathcal{N}_{k,G}(u)$ in step 8, Algorithm 2 is essentially running Dijkstra’s algorithm on $G^{(k)}$. By Lemma 5.5, the first $k$ elements to be popped off the queue in line 9 are indeed the $k$ closest vertices to $s$ in the graph $G^{(k)}$. That is, $S$ contains $\left(v,d_{G}(s,v)\right)$ for all $v\in\mathcal{N}^{d_{G^{(k)}}}_{k,G^{(k)}}(s)$. By Lemma 5.6, $\mathcal{N}^{d_{G^{(k)}}}_{k,G^{(k)}}(s)=\mathcal{N}^{d_{G}}_{k,G}(s)$. ∎

5.1 Analysis of complexity

Let us determine the computational complexity of Algorithm 2. We shall remain agnostic for the moment about the precise implementation of the min-priority queue, and hence shall use the symbols $\mathcal{T}_{in},\mathcal{T}_{dk}$ and $\mathcal{T}_{em}$ to denote the computational complexity of insert, decreaseKey and extractMin respectively. As discussed above, the complexity of DecreaseOrInsert is also $\mathcal{T}_{in}$. Let $\mathcal{T}_{nn}$ denote the cost of a nearest neighbor query in $G$. Then the cost of a $k$ nearest neighbor query, i.e. the cost of determining $\mathcal{N}_{k,G}(u)$ as in line 8, is $k\mathcal{T}_{nn}$.

Initializing the queue in line 4 requires $k$ insertions, for a cost of $k\mathcal{T}_{in}$. Precisely $k$ extractMin operations are performed, for a total cost of $k\mathcal{T}_{em}$. DecreaseOrInsert is performed $k^{2}$ times, for a cost of $k^{2}\mathcal{T}_{in}$. Finally $k+1$ $k$-Nearest Neighbor queries are performed, for a cost of $(k+1)k\mathcal{T}_{nn}$. This gives a total cost of $k\mathcal{T}_{em}+(k+k^{2})\mathcal{T}_{in}+(k^{2}+k)\mathcal{T}_{nn}$. If the min priority queue is implemented using a Fibonacci heap, insert and decreaseKey both run in constant time (i.e. $\mathcal{T}_{in},\mathcal{T}_{dk}=O(1)$) while for extractMin $\mathcal{T}_{em}=O(\log(|Q|))$. Note that $|Q|$ never exceeds $k^{2}+k$ as at most one element is added to $Q$ during every pass through the inner for loop 9–12, which happens $k^{2}$ times. Hence $\mathcal{T}_{em}=O(\log(k))$ and we have a net cost of $O(k\log(k)+k^{2}\mathcal{T}_{nn}))$ where $\mathcal{T}_{nn}$ depends on the specifics of $G$.

Let us return to the case of primary interest in this paper; where $G$ is the complete graph on $n$ vertices and $A_{ij}=\|{\bf x}_{i}-{\bf x}_{j}\|^{p}$. By remark 5.4 to compute $\mathcal{N}^{(p)}_{k,\mathcal{X}}({\bf x}_{i})$ it will suffice to compute $\mathcal{N}_{k,G}^{d_{G}}(v_{i})$. Here, $\mathcal{T}_{nn}$ is equal to the cost of a Euclidean nearest neighbors query on $\mathcal{X}$, namely $\mathcal{T}_{Enn}$. Because $\mathcal{T}_{Enn}\gg\log(k)/k$ we get that for this case Algorithm 2 runs in $O(k^{2}\mathcal{T}_{Enn})$, as advertised in the introduction. For a totally general data set, $\mathcal{T}_{Enn}=O(Dn)$. However if $\mathcal{X}$ is intrinsically low-dimensional, which we are assuming, it is possible to speed this up. For example if $\mathcal{X}$ is stored in an efficient data structure such as a k-d tree [5] or a cover tree [6] then $\mathcal{T}_{Enn}=O(\log(n))$. Note that initializing a Cover tree requires $O(c^{d_{\max}}Dn\log(n))$ time, where $c$ is a fixed constant [6]. Hence finding $\mathcal{N}^{(p)}_{k,\mathcal{X}}({\bf x}_{i})$ for all ${\bf x}_{i}\in\mathcal{X}$ requires $O(k^{2}n\log(n)+c^{d_{\max}}Dn\log(n))$.

5.2 Extension to Longest-Leg Path Distance

A small modification to Algorithm 2 allows one to compute the $k$ nearest neighbors in the longest-leg-path distance, simply change the ‘$+$’ in line 10 to a ‘max’. This guarantees that tempDist represents the longest-leg length of the path $s\to\ldots\to u\to v$. For completeness, we present this algorithm below as Algorithm 3. The proof of correctness is analogous to Theorem 5.7, and we leave it to the interested reader.

Remark 5.8**.**

By the same arguments as in §5.1, for any $\mathcal{X}=\{{\bf x}_{1},\ldots,{\bf x}_{n}\}\subset\mathbb{R}^{D}$ one may find $\mathcal{N}^{(\infty)}_{k,\mathcal{X}}({\bf x}_{i})$ by finding $\mathcal{N}^{(\infty)}_{k,G}(v_{i})$ in the complete graph $G$ with edge weights $A_{ij}=\|{\bf x}_{i}-{\bf x}_{j}\|$. The complexity of Algorithm 3 is the same as Algorithm 2, hence one may find $\mathcal{N}^{(\infty)}_{k,\mathcal{X}}({\bf x}_{i})$ for all ${\bf x}_{i}\in\mathcal{X}$ in $O(k^{2}n\log(n)+c^{d_{\max}}Dn\log(n))$.

5.3 Comparison with Results of [9]

Recall the following definition:

Definition 5.9**.**

Let $G=(V,E,A)$ be a weighted graph. A sub-graph $H\subset G$ is called a $1$-spanner of $G$ if $H$ has the same vertex set as $G$ and, for all $u,v\in V$ we have that $d_{G}(u,v)=d_{H}(u,v)$.

When preparing this manuscript for publication, the authors became aware of the following result of Chu, Miller and Sheehy:

Theorem 5.10**.**

Let $\mathcal{M}\subset\mathbb{R}^{D}$ be a compact Riemannian manifold of dimension $d$, and let $\mathcal{X}=\{{\bf x}_{1},\ldots,{\bf x}_{n}\}$ be sampled from $\mathcal{M}$ according to a Lipschitz continuous probability distribution $\mu$ satisfying $\mu_{\min}>0$. Let $G$ be the complete graph on $V=\mathcal{X}$ with $A_{ij}=\|{\bf x}_{i}-{\bf x}_{j}\|^{p}$. If $p\geq 2$ and $k=O(2^{d}\log(n))$ then with probability $1-o(1)$ the Euclidean $k$ nearest neighbors graph $G^{(k)}$ is a $1$-spanner of $G$.

Proof.

See Theorem 1.6 and Corollary 6.1.1 of [9]. ∎

This theorem give another way to deduce Theorem 5.7, as the statement that $G^{(k)}$ is a $1$-spanner of $G$ implies Lemma 5.6. However Lemma 5.6 holds in more generality. In particular:

It is not conditional (ie holds with probability $1$). 2. 2.

It places no restriction on $k$. 3. 3.

It holds for any data model. Note that Theorem 5.10 does not hold for the Data Model of §2.1, where multiple manifolds are under consideration.

6 Numerical Experiments

In this section we verify that using a $p$-wspm in lieu of the Euclidean distance does indeed result in more accurate clustering results, at a modest increase in run time. Specifically, we consider eight datasets, and compare the accuracy of spectral clustering using $k$-NN graphs constructed using $p$-wspm’s for $p=2,10$ and $\infty$ and using the Euclidean metric. As a baseline, we also consider spectral clustering with a full similarity matrix based on the Euclidean metric. For notational reasons it is convenient to denote the Euclidean metric as $d^{(1)}_{\mathcal{X}}$ (which is correct by Theorem 3.4), in which case the four metrics under consideration are $d^{(1)}_{\mathcal{X}},d^{(2)}_{\mathcal{X}},d^{(10)}_{\mathcal{X}}$ and $d^{(\infty)}_{\mathcal{X}}$. All numerical experiments described in this section were implemented in MATLAB on a mid 2012 Mac Pro with 2 2.4 GHz 6-Core Intel Xeon processors and 32 GB of RAM. All code used is available at danielmckenzie.github.io.

6.1 The Data Sets

Three Lines.

We draw data uniformly from three horizontal line segments of length $5$ in the x-y plane, namely $y=0$, $y=1$ and $y=2$. We draw $500$ points from each line to create three clusters. We then embed the data into $\mathbb{R}^{50}$ by appending zeros to the coordinates, and add i.i.d. random Gaussian noise to each coordinate (with standard deviation $\sigma=0.14$).

Three Moons.

This data set is as described in [29] and elsewhere. It has three clusters, generated by sampling points uniformly at random from the upper semi-circle of radius $1$ centered at $(0,0)$, the lower semi-circle of radius $1.5$ centered at $(1.5,0.4)$ and the upper semi-circle of radius $1$ centered at $(3,0)$. As for the Three Lines data set, we draw $500$ data points from each semi-circle, embed the data into $\mathbb{R}^{50}$ by appending zeros, and then add Gaussian noise to each coordinate with standard deviation $\sigma=0.14$.

Three Circles.

Here we draw data points uniformly from three concentric circles, of radii $1,2.25$ and $3.5$. We draw $222$ points from the smallest circle, $500$ points from the middle circle and $778$ points from the largest circle (the numbers are chosen so that the total number of points is $1500$). As before, we embed this data into $\mathbb{R}^{50}$ and add i.i.d Gaussian noise to each component, this time with standard deviation of $\sigma=0.14$.

Two dimensional projections of these data sets are shown in Figure 2. We also considered five real data sets. We focused on image data sets that are suspected to satisfy the manifold hypothesis, namely images of faces and objects taken from different angles and handwritten digits. We obtained most of our datasets from the UCI Machine Learning Repository [14].

DrivFace

consists of $80\times 80$ greyscale images of four drivers, from a variety of angles. There are $606$ images in total, and the largest class contains $179$ images while the smallest class contains $90$ images.222available at https://archive.ics.uci.edu/ml/datasets/DrivFace and see also [13]

COIL-20

The Columbia Object Image Library (COIL) contains greyscale images of a variety of objects. There are $72$ images of each object, all from different angles. The COIL-20 dataset contains all $72$ images for $20$ objects, for a total of $1440$ images.333available at http://www.cs.columbia.edu/CAVE/software/softlib/coil-20.php and see also [24].

OptDigits

This data set consists of downsampled, $8\times 8$ greyscale images of handwritten digits $0-9$ and is available at . There are $150$ images of zero, and approximately $550$ images each of the remaining digits, for a total of $5620$ images.444available at https://archive.ics.uci.edu/ml/datasets/optical+recognition+of+handwritten+digits

USPS

This data set consists of $16\times 16$, greyscale images of the handwritten digits 0–9. There are $1100$ images per class for a total of $11\ 000$ images.555available at: https://cs.nyu.edu/~roweis/data.html

MNIST

This data set consists of $28\times 28$ greyscale images of the handwritten digits 0–9. We combined the test and training sets to get a total of $70\ 000$ images.666available at http://yann.lecun.com/exdb/mnist/

6.2 Preprocessing the Data

Zelnik-Manor and Perona [30] propose the following, locally scaled kernel for spectral clustering:

[TABLE]

where $\sigma_{i}:=d({\bf x}_{i},{\bf x}_{[r,i]})$ and ${\bf x}_{[r,i]}$ denotes the $r$-th closest point in $\mathcal{X}$ to ${\bf x}_{i}$. It is shown in [30] and elsewhere that this kernel tends to outperform the unscaled Gaussian kernel (and indeed most other choices of kernel function), hence we adopt this for our experiments. We construct the full Euclidean similarity matrix as:

[TABLE]

with $r=10$. We also construct weighted $k$-NN similarity matrices for all four metrics using the same kernel and described in detail as Algorithm 4. All the real data sets were initially vectorized, so the $80\times 80$ DrivFace data set becomes vectors in $\mathbb{R}^{6400}$ and so on. Note that this procedure for constructing weighted $k$-NN similarity matrices was inspired by [21].

Remark 6.1**.**

We certainly make no claim that the choice of parameters $r=10$ and $k=15$ is optimal, and indeed playing around with them can result in slightly better (or worse) results on certain data sets. However, we observed that changing the parameters had little qualitative effect on the results, and in particular on the ordering of the similarity matrices from least to most accurate (see Table 1). As a sanity check, we also experimented with an unweighted k-NN graph, whereby for $p=1,2,10$ and $\infty$ we define $A^{(\text{uw},p)}$ as:

[TABLE]

Again, the relative ordering of the results changed little, although the accuracy was several points lower for all four metrics. All code to reproduce the experiments is available on the second author’s website, and we invite the curious reader to experiment further for themselves.

6.3 Experimental Results

For each similarity matrix we perform normalized spectral clustering as described in Ng, Jordan and Weiss [25] using freely available code 777See: https://www.mathworks.com/matlabcentral/fileexchange/34412-fast-and-efficient-spectral-clustering. We compute $k$ nearest neighbors in the $p$-wspm’s using Algorithms 2 and 3, with the data points stored in a k-d tree. We calculate the accuracy by comparing the output of spectral clustering to the ground truth and recorded the running time. For the randomly generated data sets (ie Three Lines, Three Circles and Three Moons) we ran fifty independent trials and report the mean and standard deviation. For the deterministic data sets (ie all the others) we ran ten independent trials and report the mean. The results are displayed in Tables 1 and 2. We do not attempt clustering with a full similarity matrix, $A^{(f,1)}$ on MNIST as the resulting matrix is too large to hold in memory.

From these results, we may draw several broad conclusions. Observe that for smaller or low dimensional data sets, constructing the full similarity matrix $A^{(f,1)}$ is fastest. This is due to the overhead cost of constructing the k-d tree incurred by the nearest neighbors methods. This situation is reversed in higher dimensions or for larger data sets. The gap between the run-times for $A^{(1)}$ and $A^{(2)},A^{(10)}$ and for $A^{(\infty)}$ is attributable to the cost of running Algorithm 2 or 3. Although this gap is large, relatively, for smaller data sets it becomes less relevant for larger data sets. In fact, the entire process of spectral clustering is faster with a p-wspm for the MNIST data set. We attribute this to the fact that $A^{(2)},A^{(10)}$ and $A^{(\infty)}$ are more “block-diagonal” than $A^{(1)}$, hence finding their leading eigenvectors takes less time. This more than offsets the extra time required to construct them. With regards to accuracy, observe that for every data set clustering using the Euclidean metric (either the full version or the $k$-NN version) is less accurate then using a $p$-wspm. The catch here is that which $p$-wspm varies. As a general rule, $d^{(\infty)}_{\mathcal{X}}$ appears best for elongated data (ie the three lines) or when there is a large gap between the intrinsic dimension of the data and the ambient dimension (ie USPS or MNIST). On the other hand, $d^{(2)}_{\mathcal{X}}$ appears best for more globular data (ie the three moons) or when the gap between intrinsic and ambient dimension is less pronounced (ie OptDigits). In all cases, $d^{(10)}_{\mathcal{X}}$ seems to be a good compromise between these two extremes.

Finally, note that the standard deviation of the accuracy is much higher for $d^{(10)}_{\mathcal{X}}$ (and $d^{(\infty)}_{\mathcal{X}}$) than it is for $d^{(1)}_{\mathcal{X}}$ or $d^{(2)}_{\mathcal{X}}$. This is in agreement with Theorem 4.6, where it is shown that the bound on $\epsilon_{1}$ holds with probability inversely proportional to $p$.

6.4 Varying the Power Weighting

From the analysis of §4 it would appear that taking $p$ to be as large as possible always results in the best clustering results. However, this is true only in an asymptotic sense, and indeed the results contained in Table 1 indicate that for finite sample sizes this is not always the case. In Figure 3 we show the results of varying $p$ from $1$ to $20$ for the three lines data set, this time with $300$ points drawn from each cluster. We do this for three values of the ambient dimension, $D=10,50$ and $100$. As is clear, the optimal value of $p$ depends on the dimension888The observant reader will notice that we are only varying the ambient dimension, which according to the analysis of §4 should have no effect. Recall however, that we are adding Gaussian noise of the ambient dimension, which ‘thickens’ the data manifolds and makes their intrinsic dimension weakly dependent on $D$.. In particular, observe that an intermediate value of $p$, say $p=14$, is optimal when the ambient dimension is $50$ but that a smaller power weighting ($p=2$) is more appropriate when the ambient dimension is $10$. When the ambient dimension is $100$, no power weighting performs well, which is likely because for such a large value of $D$, and such a small amount of data per cluster, the noise drowns out any cluster structure.

7 Conclusions and Future Directions

In this paper we argued that $p$-wspm’s are well-suited to the problem of clustering high dimensional data when the data is sampled from a disjoint union of low dimensional manifolds. We showed that spectral clustering with a $p$-wspm outperforms spectral clustering with Euclidean distance, and using Algorithm 2 the increase in computational burden is negligible. From the results of §6 it is clear that the geometry of the data manifolds influences which power weighting is optimal, and it would be of interest to analyze this further. Also of interest would be to extend our work to more general data models, for example those that only require the sampling distributions $\mu_{a}$ to be supported “near” $\mathcal{M}_{a}$, or those that allow for intersections between the data manifolds.

8 Acknowledgements

The first author gratefully acknowledges the support of the Department of Mathematics, The University of Georgia, where the first author was a graduate student while this work was completed. The second author thanks the Department of Mathematics, The University of Michigan for their support. Both authors thank the anonymous reviewer for many useful suggestions.