{'fraction_non_alphanumeric': 0.03570747015542503, 'fraction_numerical': 0.005640253599372363, 'mean_word_length': 4.572728346921895, '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': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Path Planning with Divergence-Based Distance Functions

Distance functions between points in a domain are sometimes used to automatically plan a gradient-descent path towards a given target point in the domain, avoiding obstacles that may be present. A key requirement from such distance functions is the absence of spurious local minima, which may foil such an approach, and this has led to the common use of harmonic potential functions. Based on the planar Laplace operator, the potential function guarantees the absence of spurious minima, but is well known to be slow to numerically compute and prone to numerical precision issues. To alleviate the first of these problems, we propose a family of novel divergence distances. These are based on f-divergence of the Poisson kernel of the domain. We define the divergence distances and compare them to the harmonic potential function and other related distance functions. Our first result is theoretical: We show that the family of divergence distances are equivalent to the harmonic potential function on simply-connected domains, namely generate paths which are identical to those generated by the potential function. The proof is based on the concept of conformal invariance. Our other results are more practical and relate to two special cases of divergence distances, one based on the Kullback-Leibler divergence and one based on the total variation divergence. We show that using divergence distances instead of the potential function and other distances has a significant computational advantage, as, following a pre-processing stage, they may be computed up to an order of magnitude faster than the others when taking advantage of certain sparsity properties of the Poisson kernel. Furthermore, the computation is"embarrassingly parallel", so may be implemented on a GPU with up to three orders of magnitude speedup.


INTRODUCTION
ATH planning in a planar domain containing obstacles is an important problem in robotic navigation. The main challenge is for an autonomous agent to move from one point (the source) in the domain to another (the target) along a realistic path which avoids the obstacles, where the path is determined automatically and efficiently based only on knowledge of the domain and local information related to the current position of the agent. This important problem has attracted much attention in the robotics community and is the topic of ongoing research. A significant family of path planning algorithms is based on socalled potential functions, which are inspired by the physics of electrical force fields, first proposed in the late 1980's by Khatib [14] and developed by Kim and Khosla [15], Rimon and Koditschek [26], and Connolly and Grupen [6] soon after. The idea is, given the target point, to construct a scalar function on the domain, such that a path to the target point from any other source point may be obtained by following the negative gradient of the function. While elegant, Koren and Borenstein [17] have identified a number of significant pitfalls that these methods may encounter, the most important being the presence of so-called "trap" situations -the presence of local minima in the potential function. To avoid this, the scalar function must have a global minimum (typically zero-valued) at the target, and be void of local minima elsewhere in the domain. The presence of "spurious" local minima could be fatal, since the gradient vanishes and the agent becomes "stuck" there. Other critical points, such as saddles, are undesirable but not fatal, since a negative gradient can still be detected by "probing" around the point.
Designing and computing potential functions for planar domains containing obstacles has been a topic of intense activity for decades. Perhaps the most elegant type of potential function is the harmonic function [15,6], which has very useful mathematical properties, most notably the guaranteed absence of spurious local minima. Alas, the main problems preventing widespread use of these types of potential functions are the high complexity of computing the function, essentially the solution of a very large system of linear equations, and the fact that very high precision numerical methods are required, as the functions are almost constant, especially in regions distant from the target. This paper addresses the first of these issues. We describe a family of new functions, which, while quite distinct from the harmonic potential function, generate exactly the same gradient-descent paths. However, they do this at a tiny fraction of the computational cost.
In practical path planning scenarios, the planar domain is described by a set of polygons representing the domain boundary and the obstacles within, which can be quite complicated. A good potential function should be "shape-aware", in the sense that it should produce paths which naturally circumvent the obstacles. The agent is armed with an automatic algorithm relying on auxiliary data structures which, given its current position in the domain, can efficiently compute the direction in which it should proceed towards the target. As we shall see later, there is a tradeoff between space and time complexity in achieving this goal.
Although the classical term is "potential function", in this paper we use the more generic term "distance function" for the guiding scalar function. We believe this is more appropriate, as in a sense, the function measures a scalar distance value between the source and the target, which the path-planner tries to decrease as it advances towards the target. Although not identical to the classical shortest-path distance (also known as "geodesic" distance), this distance also takes into account the geometry of the domain and the obstacles.
The rest of this paper is organized as follows. We start with a mathematical analysis of a number of distance functions on continuous planar domains: In Section II we consider distance functions used for path-planning which are common in the literature, based on various forms of the Laplace operator, including the classical harmonic potential function. In Section III we introduce our new family of divergence distance functions, and show (in the Appendix) that they are all equivalent to the Green's function and the Poincaré metric on the disk. In Section IV we direct our attention to the more practical case of a discretized planar domain and provide explicit algebraic expressions and computation methods for the distance functions. There we show how divergence distance functions may be computed much faster than any of the traditional distances. In Section V we provide more experimental results and insights. We conclude in Section VI with a summary and open questions.

II. LAPLACIAN-BASED DISTANCE FUNCTIONS
The Green's function Classical potential functions are based on harmonic functions, which satisfy the second-order linear differential Laplace equation is the Laplace operator, also called the Laplacian. These are particularly attractive since a harmonic function satisfies a "minimum/maximum principle" -it obtains its minimum and maximum on the domain boundary, implying that the domain interior contains no local extrema. Beyond this, harmonic functions have many other "nice" properties and have been studied extensively for decades. Rather than providing a detailed exposition here, we refer the interested reader to the book by Axler et al. [1] and the related text by Garnett and Marshall [11], which contain a wealth of information, including all the classical results we use here. Denote by Ω the open domain, by Ω its boundary, by the source point, and by the target point. Note that the mathematical translation of "obstacles" in the domain, is to "holes" in Ω. If there are no obstacles, Ω is simply connected and has a single exterior boundary loop. If there are obstacles, Ω is multiply connected having a single exterior boundary loop and multiple interior boundary loops. For convenience, we identify the plane 2 with the complex field ℂ, and much of our notation and formula will use complex number algebra. For example, a point ( , ) ∈ 2 is identified with the point = + ∈ ℂ, its conjugate is ̅ = − , its absolute value is | | = √ ̅ = � 2 + 2 , its Wirtinger derivatives The harmonic potential function is derived from the classical Green's function, which plays an important role as a kernel function in solving the Poisson equation. A Green's function with Dirichlet boundary conditions is a symmetric scalar bivariate function on the domain: ∀ ∈ Ω or ∈ Ω where the Laplacian may be applied to either of the two variables.
In the plane with no boundary, it is well known that the Green's function (or fundamental solution) is the logarithmic potential: When Ω has a boundary, it is also well known that for any fixed , may be written as: where ℎ is a harmonic function of on Ω. We say that the resulting has a pole at . In practice, ℎ is obtained by solving the following Laplace equation with Dirichlet boundary conditions: .
Given a target , the potential function, whose negative gradient is followed, is actually the negative of ( , ), which is harmonic on all of Ω except , and goes from 0 on the boundary to −∞ at Note the very important fact that must be computed anew if is changed. See an example of a Dirichlet Green's function and the gradient-descent path it generates in Fig. 1 (top left). These paths tend to stay away from boundaries. Note that a Dirichlet Green's function will typically be very flat (i.e. almost constant) away from the target point, which, as we will see later, may cause numerical issues when computed. An interesting variant of the potential function is the lesserknown Green's function with Neumann boundary conditions -(also called the Neumann Green's function or just Neumann's function of the domain). Here the situation is a little more complicated [10]: the Laplace equation becomes the more general Poisson equation and the vanishing boundary conditions are replaced by vanishing normal derivatives on the boundary. By the normal derivative of a function, we mean the component of its gradient vector in the outward direction normal to the boundary: where is the area of Ω. It is obvious from (5) that for all ≠ ∈ Ω, ( , ) is superharmonic (i.e. has a positive Laplacian) in , and is subharmonic for = . Since superharmonic functions also have a minimum principle [1], ( , ) has no local minima in Ω. It is less obvious, but rather easy to prove, that where is the area element by . As opposed to the Dirichlet Green's function , the Neumann Green's function will not be symmetric and will have a healthy spread of positive and negative values over the domain. Another important distinction between the paths generated by the Neumann Green's and those generated by the Dirichlet Green's function is that the Neumann Green's path tends to be attracted to domain boundaries. An example of a Neumann Green's function and the gradient-descent path it generates may be seen in Fig. 1

(bottom left).
We conclude by noting that it is possible to build hybrid Green's functions by mixing Dirichlet and Neumann boundary conditions to form so-called Robin's boundary conditions. The interested reader is referred to Garrido et al. [12] for a qualitative comparison of the gradient-descent paths generated by these types of Green functions and the pure Dirichlet and Neumann cases.

The heat kernel
Just as the Green's function, which is a kernel function for the Poisson equation, may be used for path planning, other kernels for more complicated linear differential equations may also be used. A useful kernel is the so-called heat kernel [3], which is the generic solution to the heat equation, having two spatial parameters and and a temporal parameter : lim →0 ( , , ) = ( − ) In the plane with no boundary, the heat kernel (or fundamental solution) is the Gaussian: When Ω has a boundary, the heat kernel of Ω , as with the Green's functions, may be defined with either Dirichlet or Neumann boundary conditions similarly to (2) or (5).
The heat kernel turns out to have properties which make it useful as a distance function. As shown by Varadhan [28], for a given target point and very small , the gradients of ( , , ) parallel the gradients of the domain geodesics to . Crane et al. [7] take advantage of this to efficiently compute geodesic distances within a domain. In practice, ( , , ) may be computed for a given and target by solving the following differential equation, in which the derivative by has been approximated by a single backward Euler step: Conveniently, the solution to this particular differential equation also has a maximum principle [25], hence the distance function, which is actually the negative of , is void of local minima within the domain. We denote by the distance function based on Dirichlet boundary conditions, and by that based on Neumann boundary conditions. See an example of and and the paths they generate for = 10 −3 in the second column of Fig. 1. Note that both decay very rapidly away from the target, which, as with the Dirichlet Green's function, may lead to numerical issues. The main qualitative difference between the two, as with the Greens' functions, is that tends to be repelled from the boundaries and tends to be attracted to boundaries. However, as noted by Crane et al. [7], in regions distant from any boundary their behavior is very similar. The observant reader will notice that the Green's functions may also be obtained as a special case of the solution to (7) as → ∞.

Spectral distances
A number of other shape-aware distance functions have been proposed in the literature based on the eigenstructure of the Laplacian on the domain. The generic form of such a distance function is: where ( , ) is the -th eigenfunction and eigenvalue of ∇ 2 . The most well-known such distance function is the so-called diffusion distance [5], which, similarly to the heat kernel-based distance, has a temporal parameter . Here ( ) = exp(−2 ). Other distances are the resistance distance ( ) = 1/ and the biharmonic distance ( ) = 1 2 defined by Lipman et al. [19]. The diffusion distance is equivalent to the 2 distance between two heat kernels with parameter centered at and , respectively. The resistance distance may not be welldefined using (8), since the infinite sum may not converge, so the correct way to define it is as an extremal length, which is known to be a conformal invariant [9]. As we shall see later, the resistance distance has a very simple equivalent definition in the discrete case [16], which has been used extensively in network analysis. For an extensive discussion and comparison between these spectral distances, the interested reader is referred to Patanè and Spagnuolo [22].
The resistance and biharmonic distances may also be defined in terms of the kernel function of their respective linear differential operators, the Laplacian ∇ 2 or the bi-Laplacian ∇ 4 , namely, if ♦ is the operator and is the solution to: then the distance is: Note that for the resistance distance is no other than the Green's function. We will denote the resistance distance by and the biharmonic distance by . Due to the infinite values ( , ) and ( , ) of the Green's function, defining through (10) is quite unstable, resulting in noisy distances, as can be seen in the third column of Fig. 1. In contrast, defining through (10) works well in most cases, although this distance has two distinct disadvantages: (1) it does not come with a guarantee of no spurious local minima, and (2) it is quite expensive to compute for a given target = , since ( , ) is required for any . As demonstrated by Lipman et al [19], the spectral definition (8) implies an efficient computation of a spectral distance by discarding the eigenvectors with small eigenvalues, but this is still quite expensive. The third column of Fig. 1 shows the gradient-descent path generated by the biharmonic distance.

III. DIVERGENCE DISTANCES
In this section we introduce a new family of distance functions which, although different from all the distance functions mentioned in the previous section, generate gradient-descent paths identical to the Dirichlet Green's function . When this happens, we say that the distance functions are equivalent to . It happens because their gradient vector field is identical in direction, but not in magnitude, to the gradient vector field of . Moreover, as we will see, the new distance functions are much easier to compute. Our point of departure is Green's theorem [1], which implies that the unique harmonic function ℎ( ) on Ω having given boundary conditions ( Ω) may be computed as the following boundary integral: where, as in the previous section, is the Dirichlet Green's function and is its normal derivate at the boundary point .
The kernel function ( , ) = ( , ) relating boundary points to interior points is called the Poisson kernel [1] of Ω. Strictly speaking, the term Poisson kernel is usually used when the domain is the unit disk, but here we will be more liberal with the terminology. ( , ) may also be obtained directly as the solution to: meaning that for all ∈ Ω, ( , ) may be thought of as a probability density on Ω. For a subset of the boundary, the quantity is called the harmonic measure of [9], and can be interpreted as the probability that a random walk in Ω starting at will exit Ω for the first time through . Conversely, if we know the Poisson kernel of Ω, the Dirichlet Green's function may be obtained from (3) and (11) as the boundary integral: To illustrate, we note that in the special case where Ω is the unit disk, we have closed formulae for the Green's functions and Poisson kernel: , (bottom) . Note that both paths are identical to the path generated by .
Unfortunately, for any other non-trivial domain, there are no closed formula for these functions, and they are typically computed numerically, as we will elaborate on in the next section. We now define a new family of distance functions using boundary integrals. These are based on the concept of -divergence introduced by Csiszár [8], which is used in statistics and information theory as a way to measure the difference between two probability measures. Assume that : ℝ + → ℝ is a continuous convex function such that (1) = 0. The divergence distance between and is the -divergence of the two harmonic measures ( , ) and ( , ): The most important property of motivates its use as a distance function for our purposes: ( , ) ≥ 0 and ( , ) = 0 iff = In the Appendix we provide an extensive mathematical analysis of divergence distances, and the most important result there is the Equivalence Theorem (Theorem A6), which proves that all distance functions based on f-divergences are equivalent to the Dirichlet Green's function, i.e. generate exactly the same gradient-descent paths. The advantage, however, of working with divergence distances is that they may be computed much more efficiently. Computation of involves a pre-processing stage where the Poisson kernels ( ,⋅) are computed for all . Then, given arbitrary and , computing ( , ) becomes a simple integral.
To be concrete, in the sequel we focus on two specific members of the divergence distance family which we believe are sufficiently different and representative: The first, called the Total Variation divergence, corresponds to ( ) = |1 − |, and is simply the 1 distance between the values of the Poisson kernel on and : The range of is [0,2]. The second distance function is the Kullback-Leibler divergence [18] between the Poisson kernels: This divergence (sometimes called the discrimination information or relative entropy) corresponds to ( ) = − log , and its range is [0, ∞).
In contrast to , which is a harmonic function everywhere except at the target point, all divergence distances are subharmonic, i.e. have non-negative Laplacians. This is a direct consequence of the convexity of : The last equality is because ( , ) is a real harmonic function of . In general, a subharmonic distance function could be bad news, as a subharmonic function possesses only a maximum principle, which prevents the occurrence of spurious local maxima, but not of spurious local minima. However, as the Equivalence Theorem shows, is equivalent to , namely, has no spurious extrema and has gradient vector fields identical in direction to the gradient vector field of . Consequently, they generate gradient-descent paths which are identical to those generated by , as demonstrated in Fig. 1. The proof of the Equivalence Theorem is based on the fact that all three distance functions are conformal invariants of the domain.
Furthermore, Section A2 of the Appendix shows that in the special case of the unit disk, the paths generated by and are hyperbolic circles, which are symmetric in the source and the target, namely, swapping the roles of these two points results in the same path, which is a desirable feature, and carries over to general domains. The same cannot be said for the other distances not equivalent to . Some of these paths for the unit disk are shown in Fig. A1 in the Appendix.

IV. DISCRETE DISTANCE FUNCTIONS
Sections II and III described a number of scalar functions, all based on the Laplacian, which could be used as domain-aware distance functions to generate gradient-descent paths to a given target point in the domain. While , , , , and have been shown to have the necessary properties, the proof relies heavily on the fact that the domain is continuous and it is not obvious that these properties carry over to practical implementations where the domain is discretized. In this section, we explore the discrete setting. As in the continuous setting, the key player is the Laplacian and the related notion of harmonicity, which fortunately is well defined also in the discrete setting.
Assume a planar domain Ω discretized into finite elements by a triangulation having vertices in its vertex set , of which are in the set of boundary vertices and are in the set of interior vertices. The discrete Laplacian of Ω is a linear operator on the triangulation vertices -a sparse symmetric × matrix whose non-zero entries corresponding to triangulation edges are given by the cotangent formula [23], and each diagonal entry is the negative of the sum of the corresponding row entries: where and are the angles opposite the edge between vertices and .
This formula is the classical conformal Laplacian. It is well known that if is a Delaunay triangulation [2] of the domain, the entries of corresponding to interior edges are guaranteed to be positive, which is a critical property of a Laplacian.
For the matrix to best approximate the continuous Laplacian in the plane, a more accurate version is the non-symmetric matrix: where is the area of the dual Voronoi cell [2] of the -th vertex. In practice, is approximated as one third of the areas of the triangles incident on the -th vertex. In a more compact form, this can be written as: is a vector of the areas. The literature also contains a "normalized" non-symmetric version of the Laplacian: = −1 , = ( ) which is common in probabilistic scenarios, e.g. random walk theory. Obviously has a unit diagonal. If, without loss of generality, we assume the following block structure for : then the notion of harmonicity of a real function defined on the vertices of the triangulation (represented as a vector) is: + = 0 where is the subvector of corresponding to the interior vertices and is the subvector of corresponding to the boundary vertices.

The Dirichlet Green's function
The Green's function of with Dirichlet boundary conditions is the function (matrix) : × → , such that: The main difference between this definition and that of the continuous case (1) or (2) is that the matrix is not neccesarily symmetric and its values are all finite, also on the diagonal (since the continuous singular Dirac delta function is replaced by the discrete finite Kronecker delta). It is easy to see that the unique solution to (19) is the × matrix: = −1 (20) so, given an interior vertex ∈ , it is possible to compute − the Green's function of with pole at − or alternatively, the -th column (or row) of − by solving the following linear system: where is the set of interior vertices excluding , and is a column vector of − 1 zeros, whose solution is: It is interesting to note that is equivalent to the harmonic potential function obtained by solving the Laplace equation on with boundary conditions = 1 and = : in (22) up to the constant factor . Both are discrete harmonic on . In practice, we obtain by directly solving (21). The advantage of this is that Cholesky pre-factorization of enables to efficiently solve (21) for different 's without inverting by back-substitution [24]. This is because changing the identity of changes only the right hand side of the equation.
It is also interesting to note that using either of , or when computing or merely changes the individual columns of these matrices by constant factors, giving essentially the same result for all target vertices. It is computationally advantageous to use the symmetric matrix, as this allows using efficient sparse Cholesky factorization.

The Neumann Green's function
With Neumann boundary conditions, the requirement of vanishing boundary values is replaced by the requirement of vanishing "normal derivative" values at the boundary vertices. In practice, in the discrete setting, this is approximated by the extension of the harmonic requirement also to the boundary vertices, effectively eliminating any boundary conditions. Thus, in principle, we should have = −1 but since is singular, the closest alternative is: where + is the pseudo-inverse of , is a column vector of the vertex areas and is a constant column vector of 1 √ , which is equivalent to: Unfortunately, this system is singular because of , but recalling that the rows of sum to 0, it is equivalent to the nonsingular system where is the set of all vertices excluding an arbitrary vertex . For a given target vertex ≠ , we immediately see that − the -th column of -is superharmonic at every vertex but , where it is subharmonic.
is obtained from (25) by solving a slightly simpler but symmetric system (for ) where 2 is the vector of squared entries of . The advantage of using (26) is that only the right hand side depends on , therefore similarly to (21), the system may be solved efficiently for different 's without inverting by Cholesky pre-factorization of , and back-substitution for different right hand sides. The only caveat is that can never be the target vertex.
It is worthy to note that despite the implicit Neumann boundary conditions, the normal derivative of , as computed in (24) or (25), does not vanish on the boundary. For example, when the domain is the unit disk and = 0, (15) implies that and are identical and their normal derivative on the boundary is 1 . When moves away from the origin, the normal derivative of on the boundary diminishes. The left column of Fig. 3 shows gradient-descent paths generated by and for a lowresolution triangulation of the domain of Fig. 1.

Gradient-descent path generation
In the discrete setting, the source point and target point are vertices of the triangulation and the distance function is a scalar per triangulation vertex. The simplest way to generate a gradient-descent path from to is to move from a vertex to the neighbor of that vertex which decreases the distance value the most (assuming this neighbor exists because it is not a local minimum). The resulting path is then along edges from to . However, a smoother path may be generated by cutting through triangles, at the price of a more complex computation. This is achieved by realizing that the three values of the distance function on a triangle define a unique gradient vector for that triangle. Starting at , the steepest gradient at an adjacent triangle may be followed into that triangle, exiting at one of its edges into a new triangle. (For some vertices, especially in a low-resolution mesh, such an adjacent triangle may not exist. In that case, we simply pick the adjacent edge having the largest gradient of the distance, follow it to an adjacent vertex and start the process there.) The gradient of the new triangle may then be followed from that point, and the process repeats. If the path intersects an edge incident on the target, the path terminates along that edge at the target. As with the simple edge-based path tracing procedure, should any of the vertices in the vicinity of the path be a local minimum, this process may get stuck. Assuming that this does not happen, the fact that is a global minimum guarantees that this path will eventually reach . Fig. 2 compares between the two types of paths. All the paths in the figures in the sequel relating to the discrete case were generated using the second method.

The heat kernel distance
Computing the heat kernel distance in the discrete setting involves solving a linear system, the equivalent to (7) = ( − ) −1 (29) As pointed out by Crane et al. [7] (and is the case with Green's functions), the resulting linear system is equivalent (up to a constant factor) to a linear system that can be pre-factored and back-substituted on demand for different 's. The second column of Fig. 3 shows the resulting gradient-descent paths on a triangulation as compared to those generated by the discrete Green's functions. As mentioned above, the Green's function are obtained when → ∞ (and in this case is no longer full rank, thus the pseudo-inverse should be used).

The spectral distances
Resistance distance is defined between two vertices in a graph as the effective resistance between these two vertices if each edge is modeled as a unit resistor [16]. This is sometimes called the commute time, as it is equivalent to the expected time to perform a random walk from one vertex to the other and back (averaged over all possible paths). The analog in the continuous case is called extremal length and is known to be a conformal invariant [9].
In the discrete setting, if we view the triangulation as a planar graph, the resistance distance between vertices and has the following simple expression [19]: = + − 2 , = ( ) + (30) In contrast to both Green's functions, the resistance distance cannot be solved efficiently for a given , since the diagonal value is required for the different 's, which in turns implies that the entire matrix is required, as it is difficult to compute the diagonal without computing the entire matrix.
Similarly, the discrete biharmonic distance between and has the following analogous expression using the pseudo-inverse of the squared Laplacian [19]: T E T E The divergence distances As described above, the and distances are defined via the Poisson kernel. In the discrete setting, the Poisson kernel is a × ( − ) non-negative matrix with unit sum rows (one "Poisson kernel" vector for each triangulation vertex). The -th column of represents the -th boundary vertex, and is a harmonic function on the triangulation with boundary value vector which is the binary indicator vector of that boundary vertex. With the notation of (18), this means where is the ( − ) × ( − ) identity matrix. Thus: Note that identical results are obtained if , or are used to define . Given the Poisson kernel vectors, the and distances are × matrices with the following entries: Figure 3 shows the paths generated using the and distances in the discrete setting. Note that they are identical to the paths generated by in the same figure.

V. EXPERIMENTAL RESULTS
We have implemented all the distance functions described in this paper on a variety of domains: convex and non-convex, with and without holes. The domains were discretized at various resolutions using a Delaunay triangulation generated by the Triangle software package [27].

Computational Complexity
One of the practical problems with the classical Green's functions and , as with the heat kernel distances and , given a target point , is that they are difficult to compute. In the discrete setting, as elaborated in the previous section, they involve solving a sparse linear system, which may be solved using either iterative methods, e.g. Gauss-Seidel, which are slow but require modest space, or direct methods [24], which are faster but have more storage requirements. Indeed, in the early days, direct methods on large linear systems were prohibitive. With the advent of efficient numerical solvers, sufficient resources and parallel processing hardware, these have become more accessible and a popular method to solve sparse linear systems is by prefactoring the positive-definite system matrix using e.g. Cholesky decomposition, and then solving for an arbitrary right hand side using back-substitution, which makes for a relatively efficient ( ) method. However, it also means that for a given target vertex of the triangulation, the distance function for the entire vertex set is solved for, which is obviously unnecessary, as the distance at only the vertices relevant to the path is really required. Using the or distances avoids this unnecessary computation, as once the Poisson kernel vectors have been computed in a pre-processing stage (independent of the identity of ), the boundary sum definition

Green
Heat Spectral Divergence = biharmonic distance (with Neumann boundary conditions). Note that the path based on resistance distance encounters a local minimum almost immediately, thus gets "stuck". (right) Divergence distances: (top) and (bottom) . Domain is color-coded with function value. Curves are function contours (level sets). Note that the paths generated by and are visually identical to the path generated by .
(34), (35) allows to compute the distance between any two vertices relatively efficiently (in �√ � or less, depending on sparsity) on demand. We have implemented prefactorization and backsubstitution required for computing on the CPU using the Intel Math Kernel Library (MKL) [13] and on the GPU using the NVIDIA cuSOLVER package [21]. Computation of the and distances on the CPU is straightforward and on the GPU may be done easily in CUDA thanks to its embarrassingly parallel nature.
When computing on the CPU, the and distances are computed only for relevant vertices along the generated path, whereas the distances are inevitably computed for the entire mesh. When computing on the GPU, the entire distance field is computed for a given even for and since, as in many parallel computing scenarios, it costs almost the same as computing the values only along the path. Furthermore, the sparse nature of the Poisson kernel vectors may be used to accelerate the computations. Tables 1 and 2 show CPU and GPU performance timings in computing paths on different domains using the , and distances. The domains are a variety of shapes and sizes, where the boundary size (in vertices) ranges between 2.5% (the disk) and 33% (the maze) of the total triangle mesh size. Computing and / efficiently both require preprocessing which is independent of the target vertex , and online processing which is dependent on . For the preprocessing consists of Cholesky decomposition, and the online processing is back-substitution. For and , the preprocessing consists of computation of the Poisson kernel vectors, and the online processing is computing the sums (34) or (35). The tables list also the preprocessing times, but since this is per-formed only once per input domain, they are of minor importance. The more important figure is the online processing time. When using the CPU, on all domains but the maze (see Fig. 7), we observe a speedup in the online processing time of anywhere between 2x and 100x. On the maze, in which a long path may effectively cover the entire domain, the online computation of and was 3x slower than the online computation of . However, when employing the approximation to possible because of the sparsity of the Poisson kernels, a speedup of 22x was obtained for the maze, and anywhere between 15x to 180x on the other domains. The sparsity in this maze was the best (99%) and in the disk was the worst (0%). On the other domains it was in the vicinity of 90%.
When using the GPU, the speedup of the online computation of and relative to is quite significant. Interestingly, probably due to the inherent serial nature of both the decomposition and the back-substitution operations, their GPU implementations seem to be much slower than their CPU implementation, further boosting the speedups of the online computation of and relative to . For it is anywhere between 42x and 780x, and for anywhere between 35x and 1100x. The sparse approximation for did not improve this significantly.

Precision
One of the other practical problems with the Dirichlet Green's function is that it is very flat, almost constant, when far away from the target point, as evident in Figs. 1 and 3. This makes its computation highly sensitive to numerical precision issues and in some cases can even make the result unusable. Wray et. al [25], who use an iterative Gauss-Seidel method to compute ,  show how to modify the iteration so that essentially log is computed with sufficient precision (of course, this new method is far superior to computing and simply applying the logarithm to the result). It is based on the well-known principle that it is best to compute the logarithm of a product of very small numbers as the sum of the individual logarithms. Our distance function, which we have shown to be equivalent to the Dirichlet Green's function, suffers from the same precision problem since the Poisson kernel vectors are computed similarly to . As a remedy, the method of Wray et. al [29] may be applied here as well to obtain the vectors log . These are easily combined in the sum (36), after the matrix is sparsified by eliminating very small values below a threshold. We found that a threshold of 1/√ is sufficient for an acceptable approximation, resulting in a typical 90% sparsity pattern.

Low resolution discretizations
Obviously the higher the resolution of the triangulation, the closer the discrete setting approaches the continuous case, for which we can prove the equivalence of our distance functions and to the Dirichlet Green's function . For low resolution triangulations, it is not obvious that a similar equivalence result holds, even if the discrete functions in principle have properties analogous to their continuous counterparts. The main issue is spurious local minima. The discrete Green's functions and , as the heat kernels and , have maximum principles, so will never contain spurious local minima. It is not so obvious what happens with the resistance distance , the biharmonic distance , and the and divergence distances. It seems that , and all suffer from this affliction when the resolution of the discretization is too low. See Fig. 4. We did not encounter spurious minima at all in any of our experiments with .   Figures 5 and 6 compare the different paths generated on a high resolution discretization of some interesting 2D shapes. The results are quite consistent with what we observed on other domains: The equivalence of , and , the very flat character of , , and , and the noisiness of . Maze-type domains, as explored in the robotics community, are especially challenging inputs. This is because they are typically "one-dimensional" or "long and skinny" simply-connected domains, and the paths generated are quite long. Fig. 7 shows our results on the maze featured in Fig. 2 of Wray et al. [29]. On these types of domains, the fact that any interior point is quite close to some boundary point has two interesting and quite significant consequences: 1) ( , ) ≈ log ( , ), 2) Most of the Poisson kernel vector ( ) is negligible, thus may be approximated well as a sparse vector, significantly saving time in the computation of . As mentioned above, for these types of examples, it is easy to run into precision problems when computing some of the distance functions. The heat functions decay so rapidly that they are effectively zero at any precision, yielding no gradient. The Dirichlet Green's function also decays very rapidly but the floating-point mechanism is still able

Green
Heat Spectral Divergence to yield a meaningful gradient.
is not as lucky in this respect.

VI. CONCLUSION AND OPEN QUESTIONS
We have introduced a new family of planar distance functions based on the concept of f-divergences [8] which may be used to generate gradient-descent paths useful in robotic path planning. These distance functions make use of the concept of Poisson kernels (or harmonic measure) ( , ) on the domain boundary: where is a convex function such that (1) = 0. We show that all distance functions based on f-divergences are equivalent to the Dirichlet Green's function, i.e. generate exactly the same gradient-descent paths, but may be computed much more efficiently. We have focused on two specific members of the family which are sufficiently different and representative: (Total Variance) is the simple 1 distance between two Poisson kernel vectors and is the more sophisticated Kullback-Leibler divergence between two probability vectors. A natural generalization of the distance is to use ( ) = |1 − | for integer ≥ 1. The special case = 1 coincides with , and the special case = 2 with the well-known Chi-Squared distance. A natural generalization of the distance ( ( ) = − log or ( ) = log ) is the -divergence for ≠ ±1: There are many more examples of interesting and populardivergences in the literature, e.g. the Helliger divergence: It remains to be seen which -divergence is optimal for pathplanning, where a major concern is its sensitivity to discretization.
Our Equivalence Theorem (Theorem A6 in the Appendix) has been proven for the case of a simply connected domain using the conformal invariance property. The analogous proof for the more general case of multiply connected domains seems to be more complicated and alludes us presently. However, we are confident, based on our experimental results, that the claim is true for this case too.
We have proved the Equivalence Theorem for the continuous case, and it has been experimentally observed to be true also when the domain is discretized by a dense triangulation and a discrete Laplacian used instead of a continuous one. We have detected a number of cases where spurious minima emerge in low resolutions. We speculate that there exists a condition on the minimal triangulation resolution such that the claim holds also for the discrete case above that resolution. This is an interesting topic for future research.
All of the theory and experimentation in this paper has been on planar domains. However, there exists a conformal theory (the so-called uniformization theorems [20]) also between manifold surfaces in three dimensions. Thus we speculate that results similar to ours could be obtained for path-planning on these types of surfaces, in particular terrains.
Finally, we note that another application of the conformal invariant distance functions could be in shape matching, as explored by Bronstein and Bronstein [4].

Conformal Mappings and Invariance
A distance function on a domain Ω is a non-negative real function : Ω × Ω → ℝ such that ∀ ∈ Ω, ( , ) = 0. In this Appendix we focus on conformally invariant distance functions in a simply connected domain Ω, namely a distance function Ω such that Ω ( , ) = (Ω) � ( ), ( )� where , ∈ Ω and : Ω → Ω ′ is a conformal mapping. In a nutshell, we say that a complex-valued function is conformal iff it is holomorphic (sometimes called analytic), namely a function of alone (and not of ̅ ), and, in addition, its derivative never vanishes in Ω. Conformal mappings are a cornerstone of complex analysis and have many important properties. We refer the interested reader to the classic textbook by Nehari [20] on this subject. In general the superscript of makes explicit which domain is referred to. When this is obvious, the superscript may be omitted.

The Green's function
A simple and well-known example of a conformally invariant function is − the Green's function of Ω with Dirichlet boundary conditions mentioned in Section 2: Theorem A1 [20,11]: is a conformal invariant, namely, if is a conformal mapping on Ω and , ∈ Ω, then Ω ( , ) = (Ω) � ( ), ( )� QED Although not strictly a distance function, is used as such for path planning based on gradient-descent.
In the special case where the domain is the unit disk, there is another well-known distance function which is a conformal invariant -the hyperbolic distance. This distance is actually a metric: The divergence distances Let Ω be a general simply connected domain and let Ω ( , ) be the Poisson kernel of Ω, where is a point on Ω, namely, for each ∈ Ω, Ω ( , ) = Ω ( ) is the solution to the following Laplace equation with Dirichlet boundary conditions: ∇ Ω = 0 Ω | Ω = ( ) Ω ( , ) is also known as the density of the harmonic measure [11] of Ω, namely, for a subset ⊂ Ω of the boundary, the harmonic measure of relative to is: Ω ( , ) = � Ω ( , ) Ω ( , ) is well known to be equal to the probability (density) that a random walk starting at will first exit Ω through . As expected, by definition: Harmonic measure is well known to be, like the Green's function, a conformal invariant. Namely if is a conformal map on Ω, then Ω ( , ) = (Ω) ( ( ), ( )) We focus now on the family of divergence distance functions, namely, those based on applying f-divergences [8] to the harmonic measure: where : ℝ + → ℝ is a strictly convex real function such that (1) = 0.
The distance function is obtained as the special case ( ) = − log( ) and the distance function as the special case ( ) = |1 − |.

QED
We are mostly interested in the behavior of the gradient of the divergence distance. We start by stating a simple relationship between the gradients of a conformally invariant distance function on two conformally related domains: Theorem A3: If is a conformal mapping on Ω and is a conformally invariant distance function (from a fixed point ∈ Ω) on Ω, then ∇ Ω ( ) = ∇ (Ω) ( ( )) ���� ( )

Proof:
Recall that the gradient of as a complex number is ∇ = 2 � , and apply the chain rule for a holomorphic func-

tion. QED
A gradient-descent path between two points is generated by following the (direction of the) gradient of the distance function, which is possible iff that gradient does not vanish. We now show that, if it exists, the gradient-descent path of a conformally invariant distance function is also a conformal invariant: Theorem A4: Denote by Ω ( , ) the gradient-descent path between and in a simply-connected domain Ω, generated by a conformally invariant distance function . If is a conformal mapping on Ω and , ∈ Ω, then:

The unit disk and =
We may analyze the behavior of the divergence distance functions and obtain exact expressions for these and other distance functions and their gradients when Ω = D (the unit disk). We start with the canonical case that = 0, and call this domain D 0 .