{'fraction_non_alphanumeric': 0.036328489534728704, 'fraction_numerical': 0.010798560191974404, 'mean_word_length': 4.525782688766114, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SpectroMeter: Amortized Sublinear Spectral Approximation of Distance on Graphs

We present a method to approximate pairwise distance on a graph, having an amortized sub-linear complexity in its size. The proposed method follows the so called heat method due to Crane et al. The only additional input are the values of the eigenfunctions of the graph Laplacian at a subset of the vertices. Using these values we estimate a random walk from the source points, and normalize the result into a unit gradient function. The eigenfunctions are then used to synthesize distance values abiding by these constraints at desired locations. We show that this method works in practice on different types of inputs ranging from triangular meshes to general graphs. We also demonstrate that the resulting approximate distance is accurate enough to be used as the input to a recent method for intrinsic shape correspondence computation.


Introduction
Distances on graphs are used on a daily basis, in many hardware and software applications, sometimes without us being fully aware of it. The best known example is map navigation, where one searches the best route to traverse a city, minimizing either distance or commute time. Another common application is packet routing through the internet, where two of the most commonly used algorithms (OSPF [26] and IS-IS [27]) allow multiple machines to collaboratively find the best route for a packet.
Two other common domains of application of distance on graphs are computer graphics and geometry processing, where images and shapes are often described as graphs (perhaps, with some additional structure as in the case of meshes). For example, being able to calculate distances between points on a shape allows to represent the intrinsic information in a way invariant to bending or other extrinsic deformations.
While existing methods for calculating distances on graphs are very efficient, they are still at least linear in the graph size. This may become prohibitively expensive when the graph is very big or when the number of distance calculations grows.
Our contribution. We introduce a spectral method to approximate distances on graphs, given only a few of the eigen functions of the graph Laplacian, possibly sampled at a subset of the vertex set. We see two main contributions of this paper: First, we show how to perform pairwise distance approximation in sub-linear time in the size of the graph (excluding some preprocessing). Second, unlike the previous method [11], we do not rely on a discrete divergence operator, which allows extension to types of graphs where such an operator is not trivial to define.

Distance maps on graphs
There exists a wealth of literature dedicated to distance computation on graphs. A useful dichotomy is to categorize such algorithms into heap based methods that gradually propagate distances over the graph, usually accessing its vertices in an unpredictable order, and stack based methods operating on the entire vertex set or on pre-fixed subsets thereof.

Heap based methods
Heap-based approaches for the computation of distances on graphs have been dominant until, probably, the last decade. One of the first such algorithms was suggested by Dijksrta [14], solving the shortest path problem for any directed graph with n vertices in O(n log n). Mitchell et al. [24] extended this algorithm to polyhedral meshes, providing an exact solution for the discrete geodesic distance computation in O(n 2 log n). It was later shown by Surazhsky et al. [33] that an approximation with rather low worstcase error bounds can be achieved in O(n log n).
Since graphs often constitute a discretization of a continuous domain, the question of consistency with the continuous geodesic distances naturally arises. While discrete geodesics are often inconsistent, "sampling theorems" have been derived guaranteeing the convergence of the discrete solution to the continuous one under certain conditions [5]. When a manifold is discretized as a graph satisfying these conditions, the discrete geodesic distances computed on it are guaranteed to converge to the continuous ones. As an alternative to solving a discrete geodesic problem, another line of approaches suggests to discretize directly the partial differential equation governing the distance map on the underlying continuous domain, ∇d = 1 s.t. d| source = 0, known as the eikonal equation. Sethian [30] and Tsitsiklis [36] have independently developed an efficient eikonal solver having O(n log n) runtime, which is essentially a continuous variant of Dijkstra's algorithm. This method, known as the fast marching method (FMM), was later extended from regular grids to other types of graphs like triangular meshes [19] and point clouds [22], just to mention a few. A faster version was introduced by Yatziv et al. [42], solving at O(n) at the cost of quantizing the distances.
Even though the mentioned methods differ in the types of data handled and the specific calculations performed, they are all dynamic programming algorithms and share the property of gradually propagating the distance calculation from the source(s) to the rest of the vertices. One potential advantage of this is that one can stop the calculation before covering the entire graph once the destination is achieved, saving time for 'close' distances. However, as a result of their inherent sequential structure and datadependent access to memory, these methods are hard to parallelize and implement efficiently on modern computer architectures. Furthermore, they usually cannot reuse the calculations from previous runs.

Stack based methods
The other category of algorithms, in contrast to the former, operate on the entire vertex set or fixed subsets thereof.
One line of methods that provide a fast approximation of FMM are called fast sweeping, and can be dated back to Danielsson [13] who proposed to calculate a complete Euclidean distance map using only four 'sweeps' of a 2D plane in O(n). The method was later expanded to other types of manifolds, as formally suggested by Zhao [43], after it was implicitly used in other works [44,34,35]. Even though these methods have O(n) complexity, their main drawback is that for some manifolds it might require several repetitions to achieve sufficient accuracy, up to O(2 n ) for an exact solution, as shown by Hysing et al. [17]. A highly parallel variant called parallel marching was suggested by Weber et al. [40], who proposed to decompose the surface into regular grids, solving each in a manner similar to [13]. Even though parallel marching requires a low-distortion parameterization into quadrilateral patches, it is still one of the fastest implementations of surface distance computation to date.
The newest line of works in this category follows the heat method by Crane et al. [11], which will be covered in greater detail in the sequel.
Finally, we would like to mention a recent work by Aflalo et al. [3], where spectral-based distance approximations are also performed as the input of a subspace parameterized multidimensional scaling (MDS) problem. While this approach performs some approximation in a manner seemingly similar to ours, we will show in the sequel that it does not scale to general pairwise distance approximations.

Background
Graph Laplacian. Let G = (V, E) be an undirected graph with |V | = n vertices. Let the graphs edges be further equipped with nonnegative weights {w ij }, defining the n × n weighted adjacency matrix of the graph as Additionally, the vertex degree matrix is an n × n diagonal matrix, whose elements are usually defined as The unnormalized Laplacian of a graph is defined as the n×n matrix L u = A−W. As with all Laplacians, it is easy to observe that the sum of all rows and columns is zero. For an undirected graph, the Laplacian is a symmetric positive semidefinite matrix. For other properties of the Laplacian, the reader is referred to [25]. One of the ways to normalize the Laplacian adopted here is the so-called random walk Laplacian [21,10], defined as The name stems from its relation to the row-stochastic transition probability matrix P = I − L rw .
Mesh Laplacian. A triangular mesh is a specific type of graph, realizing a homogeneous simplicial complex of order 2. Since meshes typically represent continuous geometric objects (manifolds), specific constructions of mesh Laplacians exist that consistently discretize some properties of the Laplace-Beltrami operator (LBO) of the underlying continuous manifold. Mesh Laplacians have been extensively used in geometry processing; for a review, the reader is referred to [31].
α β i j We adopt the popular cotangent weight scheme first formalized in its current form by Meyer et al. [23]. Under this scheme, the n × n stiffness matrix L c is constructed as an unnormalized Laplacian whose edge weights are with α and β being the two angles facing the edge (i, j), and the diagonal of degree matrix is the same as in (1). This scheme also defines the mass matrix M whose elements are Graph geodesics. A path on the graph is a sequence of p vertices P = {v 1 , . . . , v p } ∈ V , where each consecutive pair is one of the graph edges, i.e. (v i , v i+1 ) ∈ E. The length of a path P is the sum of all its edge weights, A shortest path (also known as minimum geodesic on a continuous manifold) between two vertices v and v is a path starting at v 1 = v and ending at v p = v whose length is minimum over all such paths. Note that a minimizer might not be unique as there might be several such paths.
As mentioned previously, sometimes graphs are used as a discretization of a continuous manifold. In these cases the above definition of a geodesic has to be relaxed in order to allow convergence of the discrete solution to the continues one [5]. Paths on a mesh, for example, will not be restricted to the edges and are allowed to traverse the faces.
Discretized gradient. We follow the definition of the intrinsic gradient for triangular meshes from [11], for the ease of comparison. Other realizations exist, e.g., by Pokrass et al. [28].
Let f be a piece-wise linear function on the mesh defined by the values f i on the vertices. Let δ be a triangle formed w.l.o.g. by the vertices 1, 2 and 3 with the corresponding function values f 1 , f 2 and f 3 . The gradient of f on δ is constant and is defined by where n s is the normal to δ, and e i is the edge opposing vertex i oriented counterclockwise. In a general (non-mesh) graph, we can no longer rely on faces, and define the gradient on edges. Specifically, an edge e = (v i , v j ) whose vertices have scalar values

The heat method
A recent seminal work by Crane et al. [11] proposed a new approach to distance calculation dubbed as the heat method. This method, which is extended in this paper, is summarized in Algorithm 1. While we focus on the specific realization of the heat method on triangular meshes, it is more general and can be applied to other types of data, as mentioned in [11].
The key idea behind the heat method stems from Varadhan's equation [37]: it can be shown that the solution u t of the heat equation input : Laplacian L, Gradient ∇, Divergence ∆, time value t 0 and origin point(s) output: Distance map d 1. Solve the heat equation (4) from the source point(s) for time value t, obtaining u t .
2. Calculate the gradient of the result u t , and normalize to unit length according to (5), obtaining g. (6) for the divergence of the normalized gradient g, obtaining d.

Solve Poisson equation
Algorithm 1: The Heat method [11] describing the heat distribution at time t is a function whose gradient direction coincides with that of a distance map from the same source (yet, unlike the distance map, the gradient magnitude ∇u t is not necessarily unit). The initial heat placement u 0 at t = 0 is a discretized Dirac distribution, having 1 in the source location(s) and 0 elsewhere. Following this observation, we can normalize the gradient of u t to unit length and solve Poisson's equation for the distance map d from the same source such that its Laplacian coincides with the divergence of the normalized gradient g We did not specify here the exact implementation of the divergence operator for triangular mesh, and refer the reader to Section 3.2 in [11] for details. Note however, that the discretization or even a consistent definition of this operator on more general graphs is a limitation of the heat method.
In what follows, we describe an alternative technique for distance map calculation similar in spirit to the heat method, yet relying only on the gradient but not on the divergence.

Proposed method -SpectroMeter
The proposed method which we dub SpectroMeter can be thought of as a variant of the heat method [11] in the spectral (frequency) domain. The details of the method are given in the following section and are summarized in Algorithm 2. The proposed method relies on several spectral properties of the Laplacian that we briefly highlight next.
Laplacian Spectrum. Spectral graph theory is a vast field with numerous and ramified applications [25,12]. For example, a recent tutorial by von Luxburg [39]  2. Calculate ∇h t the gradient of h t at these locations and normalize to unit length g = ∇h/ ∇h .
3. Find coefficients a that best approximate g, i.e. g ≈ k i=1 a i ∇φ i , according to (8).
4. Set the first coefficient(s) a i corresponding to constant φ i 's so that the result is non-negative (or has zero value at the source).

Synthesize approximate distance at desired loca
Algorithm 2: SpectroMeter -the proposed sublinear distance map approximation algorithm. The full (nonsublinear) flavor of this algorithm simply uses all the graph points rather than the given O(k) samples.
and graph-cuts. Being essentially an averaging operator, graph Laplacian excels at promoting smooth functions on the graph [31], and specifically, it has been shown that the Laplacian eigenbasis is optimal in the 2 sense for the approximation of smooth functions [2]. We denote the eigendecomposition of a Laplacian L as Lφ i = λ i φ i , with λ i and φ i being the eigenvalues and the corresponding eigenfunctions. We will assume ordering by increasing value of λ i . It is easy to show that the smallest eigenvalue λ 0 has the value 0 with a constant corresponding eigenvector and has multiplicity equal to the number of connected components in the graph.
Preprocessing. As with the other distance methods [30,11], some of the calculations have to be done for any input, and can be pre-computed and stored to reduce the amortized cost of distance computation.
The proposed technique, being a spectral method, relies on first performing eigendecomposition of the graph Laplacian. Specifically, we require the first smallest k n eigenvalues and the corresponding eigenfunctions. These low-frequency harmonics can be reused for many other tasks relying on the Laplacian spectrum. We henceforth denote the truncated eigenbasis collectively by the matrix Φ = (φ 1 , . . . , φ k ). Next, we also calculate the (intrisic) gradients of all the eigenvectors (the constant one(s) have vanishing gradients), denoting the results collectively as ∇Φ = (∇φ 1 , . . . , ∇φ k ).
Finally, we sample the graph at O(k) locations, preferably using farthest point sampling [15] (FPS) based on an 'exact' method like FMM. This sampling is used for the sublinear version of the algorithm, as detailed in the sequel.
Heat in the spectrum. The major difference between [11] and the proposed method is that the solution is done in the spectral domain. This has two implications. First, the heat equation (4) can be solved using the eigenbasis. One way to do this is to apply the Woodbury matrix inversion identity [41] to (4), which will still cost at least O(n 2 ) and might be inaccurate. We chose to approximate the solution directly through the spectrum (following, e.g., [32]) as the truncated serieŝ The functionĥ t (x s , y) approximates as the heat kernel, which measures the amount of heat that flows from a source point x s to point y after time t. For the case of multiple (or non-point) sources, the solution for a set of source points is the average of the heat kernel of all source points,ĥ t (χ s , y) = 1 s iĥ t (x i , y). We will henceforth denote the values ofĥ t (x s , y) at all vertices y, and averaged over all sources χ s in the graph as h t .
In essence (7) is the spectral approximation of the solution of the 'exact' heat diffusion (4), but is more efficient as it can be done in O(nk) instead of O(n 2 ) with k n. The approximation is accurate as long as the time constant t is sufficiently big and the residual of the truncated series (7) is small. This formulation allows further acceleration down to O(k 2 ) as described in the sequel.
The second implication is that we no longer need to solve the full n × n Poisson equation (6), but rather find a set of k Fourier coefficients a defining a low-rank approximation of the functionĝ = Φa, such that its gradient ∇d = ∇Φa is the closest in the 2 sense to the normalized gradientĝ = ∇h t / ∇h t of the approximate heat kernel, resulting in the following closed-form expression where † marks pseudo inverse. Note that this formulation does not rely on the divergence operator. Since the coefficient(s) of the constant eigenfunction(s) cannot be determined from (8) (as the corresponding gradients vanish), we set them to make the entire map non-negative or zero at the source.
Random walk. The only difference for general (nonmesh) graphs is that the heat equation is replaced with random walk process, where the notion of 'heat' is replaced with probability of arrival at a location y after t steps. We follow the definition from [39] for this probability where here {λ i , φ i } are the eigenvalues and the corresponding eigenfunctions of the random-walk laplacian L rw .
Sublinear approximation. Both the approximation of the heat kernel (7) and the solution of (8) are O(nk). However, we observe that both (7) and (8) are highly overdetermined, in the case of k n, and still can be solved unambiguously if only a subset of O(k) graph vertices is considered. This reduces the overall complexity of the method to O(k 2 ), which is sublinear in n. This excludes, of course, the pre-computation mentioned above which costs no less than o(nk) for eigen decomposition, then O(nk log n) for FPS, and finally O(k 3 ) for the matrix inversion in (8).
It is important to emphasize that this sublinear complexity is for a single pairwise distance, and increasing the number of sources or destinations will increase the overall complexity. For example, approximating the distance from one point to the entire vertex set will still cost O(nk). That said, empirical results show that our method is still considerably faster than all mentioned methods, even in the latter case.
In practice we sample the graph using FPS based on one of the 'exact' method mentioned previously. The minimum amount of samples is the number of non-constant eigenfunctions. We use an amount of distinct samples slightly higher than what makes ∇Φ in (8) numerically full rank.
A similar sampling is performed by Aflalo et al. [3], but is not limited by the basis size O(k), nor can it be generalized to full distance computation, as shown below.
Finding optimal value of t. In [11], the time value t in the heat equation was set to be the squared average edge length, based on discretization arguments. This selection seem to work in practice for most shapes even though there is no worst-case guarantee, and still failing in some isolated cases (e.g. in Figure 1 we had to use ten times the latter value).
As shown in [8], the spectral heat kernel in (7) does not depend on the discretization, but does not scale trivially with the global scaling of the shape. A scale agnostic version of (7) would then bê  Table I in [11] or from the TOSCA dataset [8]. Notice that in most examples a uniform setting m = 8 × 10 −3 is close to the shape-specific optimal parameter value (black dots) and yields mean error below 5% (dashed line). This figure was inspired by its counterpart in [11].
where tr(M) is the shape area. In a manner similar to [11], we also choose define the time value t by empirically setting a multiplier m t = m tr(M).
Testing on several shapes from various sources, we set m = 8 × 10 −3 , as can be drived from Figure 2. Note, that while in [11] the time value is discretization dependent, our time value depends only on the total area of the shape.

SMDS [3] FMM [19]
Comparison to spectral MDS. As mentioned previously, one part of the Spectral MDS (SMDS) method [3] includes distance calculations seemingly similar the Spec-troMeter. This part of the SMDS algorithm includes interpolation of distance maps sampled on a sparse set of points, based on a projection on the Laplacian eigenbasis. While [3] report very low reconstruction errors, these are true only for maps whose origin is close to one of the former samples. If one would use the same method with origin far from these samples, the result will look similar to the result of (7), i.e. non-unit magnitude gradients with the correct direction (as can be seen in the inset).

Experimental evaluation
In what follows, we assess the actual performance of the discussed approximation method compared to the 'exact' methods, and highlight some of its limitations. Two fla-  Table 1: Average error of projecting FMM distance on the Laplacian basis of the cat shape (26 × 10 3 vertices) from TOSCA [8] shown in Figure 1. Shown for different number of harmonics are the 2 and ∞ errors in % of intrinsic diameter as well as the relative error in %. The ∞ error does not decay fast due to the Gibbs phenomenon [16]. The errors are averaged over several runs with difference source points.
vors of SpectroMeter will be tested: full, which uses the entire vertex set, and sublinear, which uses only the predefined O(k) samples. For comparison, we used a Matlab mex implementation of FMM from the authors of [19], and a Matlab implementation of the Heat method accelerated by the CHOLMOD algorithm [9], used originally in [11].
We chose not to compare to the exact geodesic distance on polyhedra [24] in the large scale experiments, as its results are very similar to FMM result (less than 0.5% error), while the computational time is prohibitively slow. Errors are measured using the 2 and ∞ norms, normalized by the graph diameter. Another criterion we measure is relative error, defined as for a result d versus a reference value of d ref . We emphasize that latter error is somewhat biased against our method: errors of spectral representation usually concentrate in points where the function is non-differentiable, e.g., at the sources, where the relative error is most sensitive (small denominator). Nevertheless, we include these statistics here as it was used in previous studies [3,11,33] to which we compare.

Spectral representation accuracy
First, we test the accuracy of projecting (essentially, compressing) distance maps on the Laplacian eigenbasis. Table 1 shows reconstruction error of FMM distance maps, averaged over the shapes from Figure 2, using different numbers of harmonics. As suggested by [2], the 2 error is decaying very fast, and is quite negligible even for k = 100 harmonics. We attribute this to the fast decay of coefficients of distance maps, which can be shown to be quadratic for Euclidean domains [1]. The ∞ error does not decay as fast, which can be related to the Gibbs phenomenon [16].
Next, for comparison, we repeat the experiment only this time using SpectroMeter sublienar (full has very similar stats). The result, summarized in   tested, in all the following experiments we set k = 250.

Complexity analysis
In order to measure how all methods scale with the graph size, we ran them on an implicit shape discretized as a mesh, with gradually increasing vertex count. For each mesh, we measured the distance from 5% of the vertices to all the rest. Note that due to the size of the result, this task has quadratic complexity in the size of graph, even for our sublinear method. As can be seen in Figure 3, our approach outperforms the other ones, and the incline of the sub-linear flavor of our method is the lowest. If we include the eigendecomposition into the timing, both SpectroMeter implementations are still faster for large graphs. Table 3 shows error of the latter experiment compared to FMM. While the errors of SpectroMeter are roughly twice as high as the heat method, it scales more gracefully for larger graphs.

Distance on a nearest neighbor graph
One of the contributions of this work is that we do not need do define a divergence operator. This allows us to  apply the method on non-manifold graphs, for example a graph generated by running k nearest neighbors (k-NN) on data from a Euclidean domain. We generated a graph from 7-NN of 3 × 10 3 points in R 5 , which resulted in a graph Laplacian with similar sparsity to the meshes in other experiments. Figure 4 depicts distance maps from three source points, visualized in two dimensions using t-SNE [20]. The resulting maps are very similar to the exact ones generated using Dijkstra's method [14], but require a fraction of the computational cost.
In order to give a qualitative measure for the results, we sorted the vertices according to the resulting distance in each map separately, and measured Kendel's τ [18] as the permutation distance compared to Dijkstra's result. On average, projection onto the Laplacian eigenbasis gave 12% error, while both types of SpectroMeter gave 14% error. While this may seem high, note that random guess results in 50% error, and the probability of randomly achieving error below 30% vanishes quickly even for permutations with n = 50.

Shape correspondence
We now show our methods' result on a recent method for shape correspondence by Vestner et al. [38]. Due to lack of space, this section is not self-contained, as we only cover the parts relevant for this work. We refer the reader to the aforementioned manuscript for details. In a nutshell (skipping all theoretic justification), the latter method finds bijective correspondence by solving a linear assignment problem (LAP) over all the vertices Spectrally-Exact projected SpectroMeter SpectroMeter Dijkstra [14] Dijkstra [14] Full Sublinear Figure 4: Distance maps produced by four methods on a graph created using 7-NN in R 5 . The graph vertices are visualized in 2D using t-SNE [20]. Shown in each row is a results for a different source, color mapped from blue (small distance) to red (large distance).
where D 1 and D 2 are pairwise distance matrices, exp is element-wise exponent, Π 0 is an initial correspondence guess (possibly non bijective) and P is the set of all possible permutations of size n. LAP has an O(n 3 ) worse-case complexity, but can be solved in O(n 2 log n) operation on average using the auction method [6] (taking, e.g., a few seconds for a problem of size n = 10 3 ). In practice, the smooth geometry of the problem allows further multi-scale acceleration, which makes the LAP cost matrix sparse (zero entries indicate forbidden matches). In the particular case of triangular meshes, truncating the pairwise distance matrices at 5% of the diameter we can achieve about 99% sparsity, allowing solution of a large LAP problem of size n = 10 4 in about one minute. The main bottleneck that remains is the calculation of the full pairwise distance matrices. We tested the correspondence method on the FAUST dataset [7] containing 3D scans of human shapes with 7 × 10 3 vertices. We solve LAP using the sparse auction solver publicly provided by the authors of [4]. The initial guess input Π 0 was a obtained by using the current stateof-the-art correspondence method by Rodolà et al. [29] on a coarse version of the shape with 2×10 3 vertices, and interpolating the result to full resolution using nearest neighbors. Note that the initial guess method [29] also uses the Laplacian eigendecomposition, making Spectrometer's complexity even more attractive in this case: Average run times in seconds (excluding LAP) is 260 for FMM [19], 120 for the heat method [11], compared to 42 and 24 respectively for the full and sublinear flavors of SpectroMeter. Figure 5 shows results of the method [38] using several  Figure 5: Performance of the shape correspondence method from [38] on the FAUST dataset [7], based on different distance methods as input. Dashed lines indicate that the distance was projected onto the LBO basis (i.e. compressed). SperctroMeter yields slightly lower results due to the fact it is an approximation, but can further scale up to very big shapes.
distance methods as input. We also tested the effect of projecting the distance matrix onto the Laplacian eigenbasis. While the performance is lower for Spectrometer compared to the exact distance, the overall method of [38] becomes more scalable as it can be run for shapes of very high vertex count.

Conclusion
We presented an algorithm allowing to approximate pairwise distances on sparsely-connected graphs (such as, for example, triangular meshes) in sub-linear time given the low-frequency eigenfunctions of the Laplacian. Empirical evidence exhibits very good 2 approximation errors and reasonable ∞ errors.
Similarly to both the heat method [11] and FMM [30], the obtained distance map might violate the triangle inequality or distance symmetry. For similar reasons, it is possible that these distance maps cannot be used to compute the geodesics themselves via back-tracking, due to possible existence of local extrema along the path.
The influence of boundary conditions that has been ignored here deserves attention. Additionally, we would like to devise a criterion more robust than 2 for the solution of (8). Lastly, but most importantly, rigorous theoretical bounds on the average and worst-case approximation errors shall be developed.