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