Computing geodesic distances with diffusion

Question

I am trying to solve an APSP (All-Pair Shortest Path) problem on a weighted graph. This graph is actually a 1, 2 or 3 dimensional grid, and the weights on each edge represent the distance between its two vertices. What I want to have is the geodesic graph distance (shortest path through the graph), for every pair of vertices.

I want a diffusion-based method, because it is faster than a Dijkstra, or a Floyd-Warshall algorithm. I'm trying to achieve this using the heat equation: $$\frac{du}{dt} = \Delta u.$$ In the end, my application needs a kernel of the form $\exp(-d^2/\gamma)$ with $d$ the graph geodesic distance.

My hope is that since the solution is supposed to be the Green function for diffusion:

$$ u(t,x,y) = \left(\frac{1}{4\pi t}\right)^{-\frac{dim}{2}}\exp\left(\frac{-d^2(x,y)}{4t}\right),$$

then I can directly use that solution (with a few adjustements to get rid of the factor in front) as my kernel, and the parameter $\gamma$ will be adjusted by adjusting $t$.

I haven't been able to do something that works yet, and I would love some help. I have tried many things so far, and there are multiple problems that arise. It is difficult and long to explain all of them in one question, so I will first explain what I think is the beginning of a good approach, and then ask a few general questions.

---

In the same way it is done in the first step of the Geodesic in Heat algorithm by Crane et al., with a backward Euler step, I can solve the linear system: $$(Id - tL)u = u_0 \tag{1}$$ with $t$ the diffusion step, $L$ the laplacian matrix, and $u_0$ a Dirac at one of the vertices.

Solving equation (1) actually gives a kernel of the form $\exp(-d/\gamma)$, which is not desired. Therefore I have to do K subiterations in time, and solve K times: $$(Id - \frac{t}{k}L)u = u_0 \iff Mu = u_0$$ which gives $u = M^{-1}...M^{-1}u_0 \iff u = M^{-K}u_0$.

As K increases, the kernel is supposed to converge to a square one $\exp(-d^2/\gamma)$.

---

Now here come the questions :

• Should I use a graph Laplacian, or a finite differences Laplacian ? AFAIU, a graph laplacian is normalized to have 1 in the diagonal, whereas a FE Laplacian has the degree in the diagonal, and is multiplied by $\frac{1}{h^2}$

• How do I embed the graph weights in the Laplacian, so that the distance I get in the solution is the graph geodesic distance ? I want to be able to predict what will be the range of values of $d(x,y)$ in the solution, with regard to the range of the weights, and the parameters $t$, $K$, and $n$ the number of points in one direction (total domain size: $N = n^{dim}$).

• Which boundary conditions should I use in the Laplacian ? I feel like I shouldn't impose the function values (Dirichlet) at the boundary, because that would mean imposing the highest distance. Or should I ? I tried homogeneous Neumann conditions, and homogeneous Dircihlet conditions, but in either cases I get some distortion near the boundaries of the parabola $d(x,y)^2$ (which I check by computing the log of the solution $u(t)$, and substracting the minimum).

• Should I use a diffusion equation instead ? : $\frac{\partial u}{\partial t}=\nabla\cdot\left(\kappa\nabla u\right)$, since the diffusion is spatially dependent ?

References :

• Crane et al. 2013 : Geodesics in heat: A new approach to computing distance based on heat flow

Answer 1

Combinatorial Laplacian ? It depends on what you expect from your solution. Something reasonable to expect is that your solution should be independent of the way you mesh your surface (or as independent as possible). If you want that, then clearly a combinatorial graph Laplacian is not what you need, since the result would be completely different if you change the mesh.

Discrete or Continous ? A possibility is to consider that your surface is not discrete, it is rather a continuous object. The surface is parameterized by a set of discrete variables (points and function values at the vertice). For instance, a triangle $T$ is defined by:

$$ T = \{ \lambda_1 p_1 + \lambda_2 p_2 + \lambda_3 p_3\ |\ \lambda_1 + \lambda_2 + \lambda_3 = 1; \lambda_1 \ge 0; \lambda_2 \ge 0; \lambda_3 \ge 0\} $$

where $p_1, p_2, p_3$ denote the three vertices of the triangle. Note that the triangle contains an infinity of points (it is continuous !) but it is fully described by its three vertices. Similarly, when considering piecewise linear functions defined by their values at the vertices, you consider continuous functions (they are $C^0$) described by a discrete vector of their values at the vertices. It is the function space that is discrete, not the geometry.

Discretizing with FEM With this viewpoint (classical FEM), by projecting the equation onto the basis on piecewise linear functions (or other ones if you want), you can construct a FEM discretization of your equation and solve it. Refer to a classical textbook on how to do that, for instance [5] (also available on its author's webpage as pdf files in English and French), that fully explains the derivation for the Laplacian projected onto piecewise linear functions. Some care needs to be taken: too ofen computer graphics papers mix the discretization of the operator (stiffness matrix) with the "dot product" of the used function basis (mass matrix), by combining them into a single matrix (that they call "discrete Laplacian"). To fully understand what is going on, you will need to consider them separate, and use the mass matrix each time you want to project a function onto the function basis.

Varadhan's theorem and topology: Something that needs to be mentioned about the "heat method" [3] that you cite: it is based on a theorem proved by Varadhan [1]. The original article by Varadhan proves the result by using a parameterization of the surface, and by "pulling-back" all computations in parameter space (using the chain rule). Since it uses a parameterization, the proof is only valid for objects that have a non-degenerate parameterization. In particular this imposes that the object under consideration is homeomorphic to a disk. This restriction on Varadhan's theorem is also mentioned in [2], page 400: "Varadhan's formula works within the injectivity radius, what happens when y moves to the cut locus of x ... dramatic changes take place ... strange events occur at antipodal points". The 'cut locus' is the zone where you would glue a toplogical disk to form an object of arbitrary topology (or put differently, if you grow a 'Voronoi cell' from a point, this is where the boundary of the Voronoi cell will self-collide). If you want to use something else than a topological disk, then you are on your own (no mathematical guarantee). However, [3] reports some empirical results tending to show that what you obtain is reasonable, so for practical applications it may do the job. Now if you want to prove theorems, you need to know that they will only work for topological disks (unless you find a way of extending them to more general cases).

[1] Varadhan 1967, Diffusion processes in a small time interval, Comm. Pure and Applied Math, 20, 659-685

[2] A panoramic view of Riemannian geometry, Marcel Berger, Springer, 2003

[3] Geodesics in heat, K. Crane et.al, 2013

[5] https://global.oup.com/academic/product/numerical-analysis-and-optimization-9780199205219?cc=fr&lang=en&