Mis understanding of the Heat Method

Question

I have been trying to understand a paper in CG for a while, called the Heat Method by Ken

Many things have clicked but I don't fully understand it yet. In particular.

In the following $u$ is a vector with dimensions equal to the number of vertices in a mesh.

The paper states that solving $(id - t\Delta)u_t = u_0$ is a time discretized approximation of the heat flow on a surface. Rewriting gives: $u_t = u_0 + t\Delta u_t $ or in English, the heat flow at time $t$ is equal to the heat flow at time 0 plus the laplacian of $u$ at time $t$. Which is just a backwards Euler method pretty straightforward.

Then we get that: $(M - tL_C)u = δγ,$ which, to summarize, is the same as above except we have discretized the medium as a triangular mesh.

$M$ is a square diagonal matrix where the entry $m_{i,i}$ is 2 times one third of the sum of the areas of the triangles that contain vertex $i$. In math $m_{i,i} = 2\sum_j A_j / 3$ where $A_j$ are the areas of the triangles containing vertex $i$. And $L_C$ is also a square diagonal matrix where $(L_C u)_i = \sum_j (\cot \alpha_{ij} + \cot \beta_{ij})(u_j-u_i)$ which in short, is an approximation of the laplace operator at vertex $i$.

Thus manipulating the formula: $(id - tM^{-1}L_C)u = M^{-1}δγ,$ Which is the spatially discretized form of: $(id - t\Delta)u_t = u_0$.

So, this makes sense algebraically, but now comes the part I don;t understand.

Why does this work? Why would solving that equation give the correct heat diffusion? Let me explain what I mean. In my head, a geodesic distance (or heat flow, which for the purposes of this method are the same thing) is highly dependent on the shape of a mesh. So in order to know the geodesic distance at vertex $i$ I must first know the geodesic distance at the vertices that come before it relative to the source.

This mehtod however seems to imply that, given an arbitrary mesh and an arbitrary point $p_0$ on that mesh. I can grab any arbitrary point $p_1$ on the mesh and tell you what the geodesic distance from $p_1$ to $p_2$ is without having to look at the full connectivity of the mesh to determine a connectivity graph.

I am not entirely sure if what i am asking is clear. I understand the algebra of the problem, but I cannot link how solving this equation gives you the correct heat flow everywhere on a mesh. Why can you do it in parallel without any regards for the specific connectivity information?

Answer 1

When solving the equation $(M - tL_C)u = \delta_\gamma$, you effectively have to invert the operator: $$ u = (M - tL_C)^{-1} \delta_\gamma $$ Note that while the individual operators $M$ and $L_C$ are only local, containing information about the individual vertices and edges of the mesh, the inverse operator is decidedly not local. Inversion is a global operation that incorporates information across the whole matrix being inverted, which means information across the whole mesh. (Note that $M$ has a trivial inverse since it's diagonal, but $L_C$ definitely doesn't.)

In practice, you might not literally compute the inverse matrix but rather use Gaussian elimination or some such to solve the $u$ for a given $\delta_\gamma$. In that case, it's the solving process that integrates information across the whole mesh. Effectively $L_C$ is encoding a bunch of constraints between neighboring vertices, and the solver has to satisfy all the constraints across the whole mesh at once.

(Note that it's not a trivially parallelizable job. There are methods for using parallelism to accelerate solving large linear systems, but they're going to be multi-pass methods that transfer information up and down between different "scales" of the problem. Think like parallel prefix scan, or FFT—that sort of thing.)

BTW, another perspective on this equation is that it is solving for the fixed point of $u$ in $$ u = M^{-1}(\delta_\gamma + tL_Cu) $$ This is yet another way of algebraically rewriting that same equation. In this form it looks a lot like the rendering equation $L = L_e + \int L \, f_{\text{brdf}}$! It has a similar structure, where $\delta_\gamma$ is like the emitted radiance and $L_C$ is like the scattering at surfaces. Just as in the rendering equation, you're looking for a global equilibrium solution. (And if rendering was discretized, we could similarly solve it by turning the scattering operator into a matrix and throwing the whole thing into a linear solver.)