{"id":1951,"date":"2025-10-20T22:47:37","date_gmt":"2025-10-20T22:47:37","guid":{"rendered":"https:\/\/summergeometry.org\/sgi2025\/?p=1951"},"modified":"2025-10-27T14:51:03","modified_gmt":"2025-10-27T14:51:03","slug":"the-log-diffusion-solver","status":"publish","type":"post","link":"https:\/\/summergeometry.org\/sgi2025\/the-log-diffusion-solver\/","title":{"rendered":"The Log Diffusion Solver"},"content":{"rendered":"\n
\n \"Teaser\n
Figure 1:<\/strong> Comparison of distance estimation by heat methods. Standard diffuse-then-log (middle) suffers from numerical issues when the heat is small (large geodesic distances), leading to noisy and inaccurate isolines. Our proposed log-space diffusion (right) performs diffusion directly in the log-heat domain, producing smooth and accurate estimates, closely matching the ground-truth distances (left). It is worth noting that our solver outperforms others primarily for smaller time step values. This example was done at a time step of 0.01 with 50 time steps.<\/figcaption>\n <\/figure>\n\n

SGI Mentor: Justin Solomon\n

SGI Fellows: Shree Singhi, Nhu (Jade) Tran, Haojun Qiu, Eva Kato\n\n

Introduction<\/h2>\n

The heat equation is one of the canonical diffusion PDEs: it models how “heat” (or any diffusive quantity) spreads out over space as time progresses. Mathematically, let \\(u(x,t) \\in \\mathbb{R}_{\\geq 0}\\) denote the temperature at spatial location \\(x \\in \\mathbb{R}^{d}\\) and time \\(t \\in \\mathbb{R}\\). The continuous heat equation is given by<\/p>\n \n

\n \\[\\frac{\\partial u}{\\partial t} = \\Delta u,\\]\n <\/div>\n \n

where \\(\\Delta u (x, t) = \\sum_{i=1}^{d} \\frac{\\partial^2 u}{\\partial x_i^2}\\) is the Laplace operator representing how local temperature gradients drive heat from hotter to colder regions. On a surface mesh, the spatial domain is typically discretized so that the PDE reduces to a system of ODEs in time, which can be integrated using stable and efficient schemes such as backward Euler. Interestingly, in many applications\u2014including geodesic distance computation 1<\/a><\/sup> and optimal transport 2<\/a><\/sup>\u2014it is often more useful to work with the logarithm<\/em> of the heat rather than the heat itself. For example, Varadhan’s formula states that the geodesic distance \\(\\phi\\) between two points \\(x\\) and \\(x^*\\) on a curved domain can be recovered from the short-time asymptotics of the heat kernel:<\/p>\n \n

\n \\[\\phi(x, x^{*}) = \\lim_{t \\to 0} \\sqrt{-4t \\,\\log k_{t,x^{*}}(x)},\\]\n <\/div>\n \n

where \\(k_{t,x^{*}}(x) = u(x,t)\\) is the heat distribution at time \\(t\\) resulting from an initial Dirac delta at \\(x^*\\), i.e., \\(u(x,0) = \\delta(x-x^{*})\\). In practice, this leads to a two-step workflow: (1) simulate the heat equation to obtain \\(u_t\\), and (2) take the logarithm to recover the desired quantity. However, this approach can be numerically unstable in regions where \\(u\\) is very small: tiny relative errors in \\(u\\) can be greatly magnified in the log domain. For example, an absolute error of \\(10^{-11}\\) in \\(u\\) becomes an error of roughly \\(25.33\\) in \\(\\log u\\).<\/p>\n \n

A natural idea is therefore to evolve the logarithm of the heat directly. Let \\(u>0\\) satisfy the heat equation and define \\(w = \\log u\\). One can show that \\(w\\) obeys the nonlinear PDE 3<\/a><\/sup><\/p>\n \n

\n \\[\\frac{\\partial w}{\\partial t} = \\Delta w + \\|\\nabla w\\|_2^2.\\]\n <\/div>\n \n

This formulation has the appealing property that it avoids taking the logarithm after<\/em> simulation, thereby reducing the amplification of small relative errors in \\(u\\) when \\(u\\) is tiny. However, the price of this reformulation is the additional nonlinear term \\(\\|\\nabla w\\|_2^2\\), which complicates time integration as each step requires iterative optimization, and introduces a discrepancy with the discretized heat equation. In practice, evaluating \\(\\|\\nabla w\\|_2^2\\) requires a discrete gradient operator, and the resulting ODE is generally not<\/em> equivalent to the one obtained by first discretizing the heat equation for \\(u\\) (which involves only a discrete Laplacian) up to a log transform. This non-equivalence stems from (1) there is a new spatially discretized gradient operator, and (2) the fact that the logarithm is nonlinear and does not commute with discrete differential operators.<\/p>\n \n

Motivated by this discrepancy, we take a different approach that preserves the exact relationship to the standard discretization. We first discretize the heat equation in space, yielding the familiar linear ODE system for \\(u\\). We then apply the transformation \\(w = \\log u\\) to this discrete<\/em> system, producing a new ODE for \\(w\\). Now, by construction, the new ODE is equivalent to the heat ODE up to a log transform of variable. This still allows us to integrate directly in log-space, which allows more accurate and stable computation. We show that this formulation offers improved numerical robustness in regions where \\(u\\) is very small, and demonstrate its effectiveness in applications such as geodesic distance computation.<\/p>\n\n \n \n \n \n \n \n
Table 1:<\/strong> Side-by-side procedures for computing \\(\\log u\\) in heat diffusion, note the highlighted steps are the same for standard heat diffusion and ours.<\/caption>\n
Standard (post-log)<\/th>\n Continuous log-PDE<\/th>\n Proposed discrete-log<\/th>\n <\/tr>\n
1. Start with the heat PDE for \\(u\\).<\/td>\n 1. Start with the log-heat PDE for \\(w=\\log u\\).<\/td>\n 1. Start with the heat PDE for \\(u\\).<\/td>\n <\/tr>\n
2. Discretize in space \\(\\Rightarrow\\) ODE for \\(u\\).<\/td>\n 2. Discretize in space for \\(\\nabla\\) and \\(\\Delta\\) \\(\\Rightarrow\\) ODE for \\(w\\).<\/td>\n 2. Discretize in space \\(\\Rightarrow\\) ODE for \\(u\\).<\/td>\n <\/tr>\n
3. Time-integrate in \\(u\\)-space.<\/td>\n 3. Time-integrate in \\(w\\)-space and output \\(w\\).<\/td>\n 3. Apply change of variables \\(w=\\log u\\) \\(\\Rightarrow\\) ODE for \\(w\\).<\/td>\n <\/tr>\n
4. Compute \\(w=\\log u\\) as a post-process.<\/td>\n <\/td>\n 4. Time-integrate in \\(w\\)-space and output \\(w\\).<\/td>\n <\/tr>\n <\/table>\n\n

Background<\/h2>\n

We begin by reviewing how previous approaches are carried out, focusing on the original heat method 1<\/a><\/sup> and the log-heat PDE studied in 3<\/a><\/sup>. In particular, we discuss two key stages of the discretization process: (1) spatial discretization<\/em>, where continuous differential operators such as the Laplacian and gradient are replaced by their discrete counterparts, turning the PDE into an ODE; and (2) time discretization<\/em>, where a numerical time integrator is applied to solve the resulting ODE. Along the way, we introduce notation that will be used throughout.<\/p>\n \n

Heat Method<\/h3>\n

Spatial discretization.<\/strong> In the heat method 1<\/a><\/sup>, the heat equation is discretized on triangle mesh surfaces by representing the heat values at mesh vertices as a vector \\(u \\in \\mathbb{R}^{|V|}\\), where \\(|V|\\) is the number of vertices. A corresponding discrete Laplace operator is then required to evolve \\(u\\) in time. Many Laplacian discretizations are possible, each with its own trade-offs 4<\/a><\/sup>. The widely used cotangent Laplacian<\/em> is written in matrix form as \\(M^{-1}L\\), where \\(M, L \\in \\mathbb{R}^{|V| \\times |V|}\\). The matrix \\(L\\) is defined via<\/p>\n \n

\n \\[L_{ij} = \\begin{cases} \\frac{1}{2}(\\cot \\alpha_{ij} + \\cot \\beta_{ij}) \\quad & i \\neq j \\\\ -\\sum_{i \\neq j}L_{ij} & i = j \\end{cases},\\]\n <\/div>\n \n

where \\(i,j\\) are vertex indices and \\(\\alpha_{ij}, \\beta_{ij}\\) are the two angles opposite the edge \\((i,j)\\) in the mesh. The mass matrix \\(M = \\mathrm{diag}(m_1, \\dots, m_{|V|})\\) is a lumped (diagonal) area matrix, with \\(m_i\\) denoting the area associated with vertex \\(i\\). Note that \\(L\\) has several useful properties: it is sparse, symmetric, and each row sums to zero. This spatial discretization converts the continuous PDE into an ODE for the vector-valued heat function \\(u(t)\\):<\/p>\n \n

\n \\[\\frac{\\mathrm{d}u}{\\mathrm{d}t} = M^{-1}L\\,u.\\]\n <\/div>\n \n

Time discretization and integration.<\/strong> The ODE is typically advanced in time using a fully implicit backward Euler scheme, which is unconditionally stable. For a time step \\(\\Delta t = t_{n+1} – t_n\\), backward Euler yields:<\/p>\n \n

\n \\[(M – \\Delta t\\,L)\\,u_{n+1} = M\\,u_n,\\]\n <\/div>\n \n

where \\(u_n\\) is given and we solve for \\(u_{n+1}\\). Thus, each time stepping requires solving a sparse linear system with system matrix \\(M – \\Delta t\\,L\\). Since this matrix remains constant throughout the simulation, it can be factored once (e.g., by Cholesky decomposition) and reused at every iteration, greatly improving efficiency.<\/p>\n \n

Log Heat PDE<\/h3>\n

Spatial discretization.<\/strong> In prior work 3<\/a><\/sup> on the log-heat PDE, solving this requires discretizing not only the Laplacian term but also the additional nonlinear term \\(\\|\\nabla w\\|_2^{2}\\). After spatial discretization, the squared gradient norm of a vertex-based vector \\(w \\in \\mathbb{R}^{|V|}\\) is computed as<\/p>\n \n

\n \\[g(w) = W \\left( \\sum_{k=1}^{3} (G_k w) \\odot (G_k w) \\right) \\in \\mathbb{R}^{|V|}\\]\n <\/div>\n \n

where \\(G_{k} \\in \\mathbb{R}^{|F| \\times |V|}\\) is the discrete gradient operator on mesh faces for the \\(k\\)-th coordinate direction, and \\(W \\in \\mathbb{R}^{|V| \\times |F|}\\) is an averaging operator that maps face values to vertex values. With this discretization, the ODE for \\(w(t) \\in \\mathbb{R}^{|V|}\\) becomes<\/p>\n \n

\n \\[\\frac{\\mathrm{d}w}{\\mathrm{d}t} = M^{-1}L\\,w + g(w)\\]\n <\/div>\n \n

Time discretization and integration.<\/strong> If we apply fully implicit Euler, we would find \\(w^{n+1}\\) by solving a nonlinear equation. To avoid this large nonlinear solve, prior work 3<\/a><\/sup> splits the problem into two substeps:<\/p>\n \n