extfieldsofstudy:
Computer Science
Mathematics

text:
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).


BP-MAC: Fast Authentication for Short Messages

Resource-constrained devices increasingly rely on wireless communication for the reliable and low-latency transmission of short messages. However, especially the implementation of adequate integrity protection of time-critical messages places a significant burden on these devices. We address this issue by proposing BP-MAC, a fast and memory-efficient approach for computing message authentication codes based on the well-established Carter-Wegman construction. Our key idea is to offload resource-intensive computations to idle phases and thus save valuable time in latency-critical phases, i.e., when new data awaits processing. Therefore, BP-MAC leverages a universal hash function designed for the bitwise preprocessing of integrity protection to later only require a few XOR operations during the latency-critical phase. Our evaluation on embedded hardware shows that BP-MAC outperforms the state-of-the-art in terms of latency and memory overhead, notably for small messages, as required to adequately protect resource-constrained devices with stringent security and latency requirements.

• Since caching comes with a significant memory trade-off, we introduce memory optimizations realizing a smaller memory footprint than BP-MAC's closest contender (UMAC [15]) for messages smaller than 12 to 43 bytes, depending on tag lengths.
• We evaluate BP-MAC's performance on different processor architectures using the Zolertia RE-Mote and Z1 boards to show BP-MAC's capability to reduce the overhead of integrity protection by more than an order of magnitude for small messages.
Availability Statement. The source code underlying this paper is available at: https://github.com/fkie-cad/bpmac

MESSAGE AUTHENTICATION CODES
In the following, we formally introduce Message Authentication Codes (MACs) (Sec. 2.1) and universal hashing (Sec. 2.2) as a basis for efficient, secure authentication schemes. Then, we discuss related work on optimizing MACs for low-latency scenarios (Sec. 2.3).

Formal Definition
A MAC scheme is composed of two algorithms, Sig and Vrfy , that generate and verify authentication tags with the help of a pre-shared secret key [8]. A message , extended by a tag , computed as = Sig ( ), allows the receiver of a message to verify the integrity of the received message. Upon reception of a message-tag pair ( , ), the receiver uses the algorithm Vrfy ( , ) to verify whether the message or the tag has been manipulated in transit by an attacker not possessing the secret key . In most popular MAC schemes, Vrfy ( , ) computes the tag * = Sig ( ) based on the received message and returns whether * is identical to the received tag . Consequently, MAC schemes are considered secure if, even under a chosen-message attack, an attacker cannot create an existential MAC forgery, i.e., a tag for a previously unseen message , that Vrfy ( , ) would accept. Achieving this property typically requires the execution of (for resource-constrained devices) processing-intensive cryptographic algorithms.

The Carter-Wegman MAC Construction
The Carter-Wegman construction allows building efficient and secure MAC schemes, such as UMAC [15] and Poly1305 [6], which gained popularity due to their superior performance compared to traditional schemes. Carter-Wegman MAC constructions add an additional nonce to the parameters of Sig and Vrfy , which natively protects against replay attacks [8]. These constructions first compute the digest of the message with a hash function unknown to the attacker. Then, this digest is masked by the encrypted nonce. More formally, a Carter-Wegman MAC construction computes the tag for a message and a nonce as = 1 ( ) ⊕ 2 ( ), where is a pseudorandom function covering the same output space as the universal hash function [8]. A universal hash function is a (not necessarily cryptographic) secret hash function for which it is hard to find a collision, i.e., finding , ′ such that ( ) = ( ′ ), as long as the attacker does not know a single output of [8].

Efficient MACs for Low Latency Scenarios
Optimized computations of established primitives such as AES or SHA256 have been proposed to address the need for fast message authentication on resource-constrained devices. Such approaches include hardware acceleration [26] or the preprocessing of authentication for predictable parts of future messages [13]. However, these approaches still fail to enable secure sub-millisecond communication as many CPS scenarios demand [11,18]. Consequently, a range of novel MAC algorithms based on new cryptographic primitives has been proposed to reduce processing and latency overheads for resource-constrained devices, including SipHash [4], TuLP [12], Chaskey [21], and LightMAC [19]. MergeMAC [2] saves valuable time during latency-critical phases by extracting predictable parts of future messages to authenticate them in advance. However, missioncritical CPS must rely on scrutinized and time-proven cryptographic primitives to ensure the maximum possible security guarantees.
Therefore, competing approaches target efficient Carter-Wegman MAC schemes, whose security can be reduced to underlying cryptographic primitives (e.g., AES). To do so, UMAC [15] and VMAC [14] use custom hash functions optimized for 32-bit and 64-bit architectures, respectively. The universal hash function employed by Poly1305 evaluates a polynomial with the message as coefficients over the finite field Z 2 130 −5 [6]. Although highly efficient, these schemes optimize for messages of at least a few hundred bytes [4]. PMAC [7] fragments messages to parallelize and thus speed up tag computations. Nevertheless, typical resource-constrained devices with single-core processors do not benefit from such optimizations.
Another branch of research tackles the reduction of bandwidth overhead introduced by MACs schemes via progressive MACs (Pro-MACs [3]). ProMACs reduce the tag sizes by aggregating authentication over several messages [3,16,17,23,24]. Thereby, ProMACs protect each message with immediately reduced security, allowing their optimistic processing without waiting for the reception of subsequent messages. Still, a message's authenticity is reinforced with subsequent messages to ensure eventual strong security [3]. Like traditional authentication, ProMACs also benefit from faster underlying MAC schemes to reduce processing latency.

BP-MAC: A BITWISE PRECOMPUTED MAC
To remove the bottleneck of cryptographic processing for message authentication and thus enable low-latency transmissions of small messages on resource-constrained devices using conventional security primitives, we propose Bitwise Precomputable MAC (BP-MAC). BP-MAC combines the idea of preprocessing message-independent computations with a Carter-Wegman MAC construction specifically designed for small messages. We begin by discussing preprocessing as the underlying principle of BP-MAC (Sec. 3.1) before detailing BP-MAC's construction (Sec. 3.2). Then, we present memory optimizations (Sec. 3.3) and discuss the security of BP-MAC (Sec. 3.4).

Precomputations for Fast Authentication
BP-MAC takes advantage of precomputations since the processing loads in CPS scenarios are not evenly distributed over time [2,13]. Instead, devices periodically have idle times until new data has to be transmitted, received, or processed. Hence, the idea behind BP-MAC is to offload computationally expensive operations into those idling phases to reduce computations during latency-critical phases to a minimum. Therefore, a naïve approach would be to precompute and store authentication tags for each possible message. However, precomputing authentication information, e.g., for all possible, rather nonce n masking tag message m Figure 1: BP-MAC achieves high-speed tag computations during the latency-critical phase by only XORing the bit tag depending on each bit value and the masking tag to generate a secure authentication tag. The masking tag is independent of the message and can thus be precomputed during idle times. Bit tags are precomputed once for each unique key.
tiny two-byte messages would then already consume 1 MB of memory (assuming the recommended tag length of 16 bytes), a number that is exponentially growing for longer messages. Hence, this naïve approach is infeasible for resource-constrained CPS devices. In contrast, the core idea of BP-MAC is to precompute tags for individual bits, which are then efficiently combined when a message needs authentication. As each bit can only represent either 0 or 1, BP-MAC's memory consumption grows linearly with message lengths and amounts to only 512 bytes for all possible two-byte long messages. These tags can be precomputed once a new key is exchanged or even provided by the (more powerful) communication partner. In addition, BP-MAC masks the combined tags with a secret masking tag to prevent replay attacks and enable their secure aggregation. Such tags can be conveniently precomputed between transmissions since they are independent of the actual message.

BP-MAC's Design in Detail
BP-MAC achieves fast authentication by only executing the efficient aggregation of previously precomputed tags during the latencycritical phase. We illustrate the different steps of BP-MAC's design in Figure 1 and discuss them in detail in the following. Figure 1, bit tags are intermediate tags computed for each message bit. To ensure BP-MAC's security and prevent collisions, all bit tags have to be a unique pseudorandom number that must be kept secret. BP-MAC computes bit tags as the digest of a pseudorandom function by concatenating bit index and bit value. Concretely, BP-MAC uses AES-128 1 as a pseudorandom function. Here, 1 is a secret key shared between both communicating parties due to its increased efficiency over other functions, e.g., HMAC-SHA256. Thus, if the value of the third bit in the second byte is zero, the corresponding bit tag is AES-128 1 (10||0). Overall, there exist 2 · -bit tags that have to be precomputed for a message with bits. However, this only occurs once a new key is established and could even be offloaded to a more powerful device if necessary.

Masking
Tags. The purpose of masking tags is twofold. First, they provide replay protection by incorporating a nonce into the final tag. Second, these pseudorandom tags prevent the leakage of information concerning individual bit tags and are thus crucial for the security of BP-MAC. As nonce, BP-MAC uses a zero-initialized counter that is incremented for each newly computed tag.
Computing the masking tag from a pseudorandom function and the nonce is analogous to UMAC [15] and relies on AES-128 2 . For tags with a length between 9 and 16 bytes, we use the leading bytes of the AES-encrypted nonce, where 2 (≠ 1) is also a secret key shared between both parties. For smaller tags, we use the cached 16-byte long AES output for multiple masking tags, as long as no output bytes are reused. Thus, for tags of 4 bytes, only every fourth nonce has to be encrypted to reduce processing overhead, while we can use the unused bytes of the cached output to blind the remaining tags. As nonces can be predicted, masking tags can be precomputed and thus do not add to the latency-critical delay.

Tag Computation and Verification.
After discussing the two precomputation phases generating the bit and masking tags, we now explain the computation of the -th authentication tag for the message . As shown in Figure 1, the final step of tag computation consists of XORing the intermediate tags. This observation also manifests in the formal definition of computing tag : The second part is the masking tag that is XORed with the bit tags that correspond to the message. Here | | denotes the length of the message and [ ] denotes the value of the -th bit. Thus, all computations besides these final XOR operations are independent of the message that should be authenticated. Therefore, BP-MAC's construction enables the quick computation of new tags, especially for short messages which require a low number of XOR operations.
As usual for deterministic MAC algorithms, the tag verification happens by computing the tag * for the received message and then comparing it to the received tag . If both tags are identical, the receiver can conclude, with high confidence, that the received message has not been altered in transit.

Memory Optimizations for BP-MAC
The naïve realization of BP-MAC requires the storage of 2 · | | bit tags. For a fixed-length messages, we can, however, half this memory overhead by precomputing the default tag of a message composed of only zeros. To support the same optimization for variable-length messages, we must pad the message before authenticating it. Then, we only need to store how to alter the tag for a bit set to one with so-called bitflip tags, of which we need | | in total (one for each bit). We sketch the tag computation with this optimization in Algorithm 1, highlighting in blue the computations executed during the latency-critical phase. Already before a new message is ready, we can compute the masking tag and XOR it with the default tag t default to generate the tag of an all-zero message as: ⊲ Padding for the message return t Then, to compute the tag for the actual message, we need to change the tag for each bit that is not zero as assumed for the default tag t default . To efficiently realize these changes, we can precompute bitflips tags for all bits in a message as: By XORing these bitflip tags to the default tag t default , we compute the tag as if the bit at the corresponding position is set. Finally, we pad the message to a fixed length to be able to differentiate between a shorter message and one with trailing zeros. BP-MAC uses the Padding method 2 as described by ISO/IEC 9797-1 [1], appending a 1-bit at the end of the message followed by as many 0-bits as necessary to reach the desired message length. As the default tag t default already assumes that all bits are zero, we only need to incorporate a single further bitflip tag t for the first index after the end of the message. Thus, we only need to store | | + 3 tags (one for each bit, including padding, t nonce , and t default ). This procedure reduces the number of operations during the latency-critical phase for each 0-bit, which accounts for about half of encrypted traffic.

Security Discussion
BP-MAC utilizes the secure Carter-Wegman MAC construction [9,25] in which a tag for a message and a nonce is computed as where is a pseudorandom function covering the same output space as the universal hash function (cf. Sec. 2.2). For , BP-MAC uses the same AES-based procedure as UMAC [15], with the sole difference that it computes the result in advance. Meanwhile, the fragmentation of messages into individually authenticatable bits is a special case of the established universal hash function ⊕ [8]. Thus, BP-MAC is secure as long as the underlying cryptographic primitive (i.e., AES) is secure. Moreover, BP-MAC can easily exchange this primitive, if ever considered not secure enough. In the following, we define the considered threat model and subsequently discuss possible attack scenarios.

Threat
Model. The attacker's goal is to alter traffic such that the recipient of a message accepts the modified message is genuine. To achieve this goal, the attacker can observe and alter the (plaintext) messages and exchanged tags to either learn the secret key used by BP-MAC or directly alter a message and the corresponding tag. The information to realize such an attack can be extracted directly from the observed traffic or through sidechannel information by observing the timing of individual packets. However, we explicitly do not consider an attacker with (partial) control over one or both communicating entities that could try to access keying information through a e.g., cache side-channel attacks. However, BP-MAC does not prevent the use of common mitigation techniques that protect against such attacks [20].

Resilience to Key Recovery
Attacks. BP-MAC is not susceptible to key recovery attacks, as otherwise a key recovery attack against the underlying cryptographic primitive, i.e., AES, would exist. By overhearing transmissions, an attacker learns = ⊕ 1 ( ) ⊕ 2 ( ) for a known message and nonce . An attacker cannot learn 2 from such tags, as otherwise there would exist a key recovery attack against 2 , i.e., AES-128: After learning a digest = 2 ( ), the attacker could generate a key ′ , consider all unique messages as nonces, and compute a BP-MAC tag ′ for arbitrary new messages ′ as ′ = ⊕ ′ ( ′ ) ⊕ 2 ( ). Thus, if an attack that recovers 2 from exists, this attack could be used to recover 2 from simple 2 digests. Similarly, 1 cannot be recovered by an attacker. Even if an attacker would learn the output of ⊕ 1 ( ), a recovery of 1 would imply the existence of a key recovery attack against AES-128: The attacker can compute ⊕ 1 ( 1 || · · · || ) for observed ciphertexts = ( ) by considering the first part of an observed plaintext to be the index and the second part to be a truncated message, i.e., = || . Then, if a key recovery attack would exist against ⊕ , this same attack would also enable recovering for AES plaintext-ciphertext pairs.

Unforgeability of Authentication
Tags. However, a successful attack does not need to recover the authentication keys to attack a MAC scheme. In many cases, it suffices if an attacker can forge authentication tags and thus inject harmful messages into a protected communication. Nevertheless, BP-MAC is also secure against such attacks. First, a message is hashed by ⊕ , with an unpredictable result for the attacker if they did not observe previous digests of ⊕ . Therefore, each output of ⊕ is masked by a unique pseudorandom number, ensuring that an attacker cannot learn a single output of ⊕ . As the nonce changes for each message, no information about the internal state of BP-MAC is revealed to an outsider. Thus, BP-MAC is also secure in the presence of an attacker that merely wants to generate valid tags for altered or new messages.

Timing
Side-Channel Attack against BP-MAC. The presented memory optimizations for BP-MAC execute some computations only for bits set to one, potentially enabling timing side-channel attacks. In particular, the time a message authentication computation takes might hint at the number of bits set to one in the authenticated message. If the message is sent in plaintext, this information is already accessible to third parties and does not leak confidential information or help in recovering keys. However, when authenticating encrypted traffic, it is crucial operating in the encrypt-then-MAC mode [8] not to reveal information about the plaintext.

PERFORMANCE EVALUATION
BP-MAC reduces the delay of message authentication for short, mission-critical wireless communication. We validate this claim

Latency Measurements
We implement BP-MAC for Contiki-NG and evaluate its performance on two different architectures to ensure the observed benefits generalize across them: Zolertia Z1 (MSP430 @ 16 MHz, 16-bit CPU, 8 kB RAM) and Zolertia RE-Mote (ARM Cortex-M3 @ 32 MHz, 32bit CPU, 16 kB RAM). To assess the performance of BP-MAC, we have to compare it to other MAC schemes with similar security guarantees, such as UMAC, VMAC, or Poly1305. We choose to compare BP-MAC against an optimized implementation of UMAC 1 as it introduces the least processing overhead for short messages (<32 bytes) even on an unfavorable architecture [14].

BP-MAC on the Zolertia Z1.
First, we compare BP-MAC's and UMAC's performance on the 16-bit MSP430 processor of the Zolertia Z1, a typical architecture for resource-constrained CPS devices. Here, we report on the time for computing 16-byte tags for messages of varying lengths. Figure 2 depicts the means and 99% confidence intervals for the respective computations of one tag. We measured 100 tag computations and derive from this the time of one tag computation due to a too low clock resolution. We repeated each measurement 30 times. Note that Sig and Vrfy require one tag computation each, such that the actual delay introduced into one secure transmission is twice the reported time.
As expected, the time required to compute UMAC tags is independent of the message length on the analyzed scale. In contrast, BP-MAC's bitwise processing introduces a linear dependency between message length and processing time. We observe that BP-MAC significantly outperforms UMAC, such that even for 32 bytes long messages, the overall processing time of BP-MAC is still 40 % lower. For 1 byte long messages, the difference is even more extreme, as BP-MAC induces a delay of only 86 for one tag computation, compared to 3.2 ms for UMAC. BP-MAC is thus significantly faster than UMAC for short messages on a 16-bit architecture.
We repeated the same measurements for shorter tags (12,8, and 4 bytes). There, we observed similar behavior, whereas the absolute processing time of BP-MAC became shorter by 24.8 − 35.4%, 38.6 − 54.0%, and 51.1 − 66.5% for increasingly shorter tags as shorter tags require less XOR operations. In turn, UMAC also becomes faster, primarily through its ability to reuse AES computations to mask multiple tags for 8 and 4-byte tags. However, the irregular need to compute AES blocks introduces jitter to the processing of UMAC,  Figure 3: Even on the faster Zolertia RE-Mote, BP-MAC realizes faster tag computations than UMAC in the latencycritical phase for messages shorter than 21 bytes and less overall processing for messages shorter than 10 bytes.
which is undesirable in low-latency communication [10]. For BP-MAC, these bursty computations only occur in the preprocessing phase and thus do not influence actual transmission delays.

BP-MAC on the Zolertia RE-Mote.
To show the performance of BP-MAC on a more powerful device, we furthermore compare BP-MAC and UMAC on the Zolertia RE-Mote. While it still constitutes an embedded device, its ARM Cortex-M3 is significantly more powerful than the Zolertia Z1. Also, its 32-bit processor is precisely the architecture targeted by UMAC. Hence, we repeat our previous measurements and report on the results in Figure 3. Still, BP-MAC outperforms UMAC for small messages, partially by more than one order of magnitude. To be precise, the latencycritical processing of BP-MAC is faster for messages not longer than 21 bytes. Furthermore, even the overall processing overhead is smaller for messages shorter than 10 bytes. These results thus show that BP-MAC is particularly suited for authenticating short messages, even in a worst-case comparison to UMAC. For shorter tags, this trend continues, with BP-MAC outperforming UMAC by 17.2−24.4%, 22.8−32.0%, and 35.2−48.3% during the latency-critical phase for tags of sizes 12, 8, and 4, respectively. The tipoff points below which BP-MAC outperforms UMAC are 26, 18, and 15-byte messages for increasingly shorter tags.
Concluding, BP-MAC enables secure communication based on the established security of AES with low processing overhead for small messages across different processor architectures. However, the performance gap between BP-MAC and UMAC is much narrower on the Zolertia RE-Mote. This behavior is not due to a worse performance by BP-MAC, but instead due to UMAC being specifically optimized for the 32-bit processor in this scenario.

BP-MAC and
Hardware-accelerated Cryptography. Due to the rising importance of low-latency security in CPS scenarios, lowpowered devices increasingly provide hardware accelerators for cryptographic operations. However, using the Zolertia RE-Mote's sha-256 accelerator for HMAC computations shows that even hardware-accelerated cryptography can introduce significant processing overheads. Indeed, our measurements reveal that the computation of one HMAC-SHA256 tag takes approximately 220 , independent of tag and payload lengths for small payloads. Thus, in these cases, hardware-accelerated HMAC-SHA256 is even slower than a software implementation of UMAC. Still, accelerators for specialized MAC schemes can be constructed to achieve even faster cryptographic operations. Fittingly, BP-MAC's heavy reliance on XOR operations and its high parallelizability by individually processing each bit promise highly efficient hardware implementations.

Memory Overhead
Fast MAC schemes based on universal hashing, such as BP-MAC and UMAC, typically trade memory usage for high processing speeds, a limited resource on low-power devices. Hence, we compare BP-MAC's and UMAC's memory footprint by compiling Contiki-NG as a native Linux application. We then use Valgrind's massif tool to analyze the peak memory footprint of both schemes for varying message and tag lengths. We depict our results in Figure 4. We notice that the memory footprint of UMAC is constant for different message sizes and hardly changes across tag lengths (ranging from 1.4 kB to 1.6 kB). In contrast, BP-MAC's memory footprint increases with the sizes of tags and messages. Overall, BP-MAC's memory footprint increases by eight times the tag lengths when messages become one byte longer, which we expected since an additional bitflip tag has to be stored for each additional message bit. Consequently, BP-MAC's memory footprint is favorable for small messages up until a tipoff point where UMAC becomes, in turn, more resource-efficient. These tipoff points lie at a message size of 12, 15, 21, and 43 bytes for 16,12,8,4 byte long tags, respectively.
Overall, our evaluation thus shows that BP-MAC outperforms UMAC in terms of processing latency and memory footprint for small messages while providing the same security guarantees, i.e., reducible to the security of AES. Thus, BP-MAC enables secure communication in critical CPS scenarios with a significantly smaller impact on communication latency than state-of-the-art approaches.

CONCLUSION
CPSs rely on low-latency wireless communications with stringent security guarantees. In this context, a significant challenge is to ensure the authenticity and integrity of exchanged messages, which typically leads to computationally intensive calculations in latencycritical phases, e.g., once a new message is ready for transmission. While numerous MAC schemes aim to reduce latency, none address small messages of only a few bytes, an essential characteristic of many latency-critical CPS scenarios. To fill this gap, we propose a new Carter-Wegman MAC scheme with AES-based security (BP-MAC), optimizing the extensive use of preprocessing with a universal hash function for the aforementioned small messages. Thus, BP-MAC reduces latency-critical computations to a few XOR operations, which overall leads to more than a ten-fold reduction in processing overhead compared to the state-of-the-art. Consequently, BP-MAC enables integrity protection of small messages in latency-critical scenarios, even for CPS devices with limited computational and memory resources.