# Mis understanding of the Heat Method

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?

• It's still a discrete timestep method, right, so it's not going to solve the whole thing in one step. It's not going between arbitrary points far apart on the mesh, but only between neighboring ones. The neighbor information is encoded in the Laplace operator. I presume you would still have to iterate this equation many times to get to an equilibrium solution for the whole mesh. – Nathan Reed Sep 13 '20 at 16:10
• $L_C$ here is not a diagonal matrix btw, it contains a nonzero entry for each neighboring pair of vertices $i, j$. – Nathan Reed Sep 13 '20 at 16:18
• Ah i missunderstood, i thought the entire sum was put inside the diagonal. But thanks to your comment I see that you can also put each term of the sum in an entry of a matrix and get the same result. – Makogan Sep 13 '20 at 18:55
• @NathanReed, The method actually only executes one step, because it's not interested in actually modelling the heat diffusion over time. It only uses it to model geodesics, as heat diffusion and geodesics are intimately related. Bot even if you were to simulate heat diffusion, my confusion remains. How can you parallelize the solution in this fashion? Why can you parallelize the solution? Why don;t you need to look at the full conectivity? Why can you just focus on the local connectivy and that;s enough? That's what doesn;t click for me. – Makogan Sep 13 '20 at 19:00

## 1 Answer

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.)