<p>When solving the equation <span class="math-container">$(M - tL_C)u = \delta_\gamma$</span>, you effectively have to invert the operator:
<span class="math-container">$$
u = (M - tL_C)^{-1} \delta_\gamma
$$</span>
Note that while the individual operators <span class="math-container">$M$</span> and <span class="math-container">$L_C$</span> 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 <span class="math-container">$M$</span> has a trivial inverse since it's diagonal, but <span class="math-container">$L_C$</span> definitely doesn't.)</p>
<p>In practice, you might not literally compute the inverse matrix but rather use Gaussian elimination or some such to solve the <span class="math-container">$u$</span> for a given <span class="math-container">$\delta_\gamma$</span>. In that case, it's the solving process that integrates information across the whole mesh. Effectively <span class="math-container">$L_C$</span> is encoding a bunch of constraints between neighboring vertices, and the solver has to satisfy all the constraints across the whole mesh at once.</p>
<p>(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 &quot;scales&quot; of the problem. Think like parallel prefix scan, or FFT—that sort of thing.)</p>
<p>BTW, another perspective on this equation is that it is solving for the fixed point of <span class="math-container">$u$</span> in
<span class="math-container">$$
u = M^{-1}(\delta_\gamma + tL_Cu)
$$</span>
This is yet another way of algebraically rewriting that same equation. In this form it looks a lot like the rendering equation <span class="math-container">$L = L_e + \int L \, f_{\text{brdf}}$</span>! It has a similar structure, where <span class="math-container">$\delta_\gamma$</span> is like the emitted radiance and <span class="math-container">$L_C$</span> 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.)</p>

CC BY-SA 4.0
2020-09-14T00:51:31.607Z
2020-09-14T00:51:31.607Z
Nathan Reed
<p>I have been trying to understand a <a href="https://www.cs.cmu.edu/%7Ekmcrane/Projects/HeatMethod/paperCACM.pdf" rel="nofollow noreferrer">paper in CG</a> for a while, called the Heat Method by Ken</p>
<p>Many things have clicked but I don't fully understand it yet. In particular.</p>
<p>In the following <span class="math-container">$u$</span> is a vector with dimensions equal to the number of vertices in a mesh.</p>
<p>The paper states that solving <span class="math-container">$(id - t\Delta)u_t = u_0$</span> is a time discretized approximation of the heat flow on a surface. Rewriting gives: <span class="math-container">$u_t = u_0 + t\Delta u_t $</span> or in English, the heat flow at time <span class="math-container">$t$</span> is equal to the heat flow at time 0 plus the laplacian of <span class="math-container">$u$</span> at time <span class="math-container">$t$</span>. Which is just a backwards Euler method pretty straightforward.</p>
<p>Then we get that: <span class="math-container">$(M - tL_C)u = δγ,$</span> which, to summarize, is the same as above except we have discretized the medium as a triangular mesh.</p>
<p><span class="math-container">$M$</span> is a square diagonal matrix where the entry <span class="math-container">$m_{i,i}$</span> is 2 times one third of the sum of the areas of the triangles that contain vertex <span class="math-container">$i$</span>. In math <span class="math-container">$m_{i,i} = 2\sum_j A_j / 3$</span> where <span class="math-container">$A_j$</span> are the areas of the triangles containing vertex <span class="math-container">$i$</span>. And <span class="math-container">$L_C$</span> is also a square diagonal matrix where <span class="math-container">$(L_C u)_i = \sum_j (\cot \alpha_{ij} + \cot \beta_{ij})(u_j-u_i)$</span> which in short, is an approximation of the laplace operator at vertex <span class="math-container">$i$</span>.</p>
<p>Thus manipulating the formula: <span class="math-container">$(id - tM^{-1}L_C)u = M^{-1}δγ,$</span> Which is the spatially discretized form of: <span class="math-container">$(id - t\Delta)u_t = u_0$</span>.</p>
<p>So, this makes sense algebraically, but now comes the part I don;t understand.</p>
<p>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 <span class="math-container">$i$</span> I must first know the geodesic distance at the vertices that come before it relative to the source.</p>
<p>This mehtod however seems to imply that, given an arbitrary mesh and an arbitrary point <span class="math-container">$p_0$</span> on that mesh. I can grab any arbitrary point <span class="math-container">$p_1$</span> on the mesh and tell you what the geodesic distance from <span class="math-container">$p_1$</span> to <span class="math-container">$p_2$</span> is without having to look at the full connectivity of the mesh to determine a connectivity graph.</p>
<p>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?</p>

CC BY-SA 4.0
2020-09-13T16:10:57.370Z
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
CC BY-SA 4.0
2020-09-13T16:18:11.377Z
$L_C$ here is not a diagonal matrix btw, it contains a nonzero entry for each neighboring pair of vertices $i, j$.
Nathan Reed
CC BY-SA 4.0
2020-09-13T18:55:08.910Z
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
CC BY-SA 4.0
2020-09-13T19:00:47.003Z
@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
<mathematics><computational-geometry><mesh><geometry><gradient>
Mis understanding of the Heat Method
Makogan
# Mis understanding of the Heat Method

(Asked by: username_0 on 2020-09-12)

I have been trying to understand a [paper in CG](https://www.cs.cmu.edu/%7Ekmcrane/Projects/HeatMethod/paperCACM.pdf) 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?

**Tags:** mathematics, computational-geometry, mesh, geometry, gradient

> **Comments:**
> **username_1**: 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.  
> **username_1**: $L_C$ here is not a diagonal matrix btw, it contains a nonzero entry for each neighboring pair of vertices $i, j$.  
> **username_0**: 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.  
> **username_0**: @username_1, 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.  

---

## 1 answers

### Answer by username_1 on 2020-09-14. Score: 1 | Accepted 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.)


---

