Editing Discrete Laplace operator (section)

==Discrete heat equation==
Suppose <math display="inline">\phi</math> describes a temperature distribution across a [[graph (discrete mathematics)|graph]], where <math display="inline">\phi_i</math> is the temperature at vertex <math display="inline">i</math>. According to [[Newton's law of cooling]], the heat transferred from node <math display="inline">i</math> to node <math display="inline">j</math> is proportional to <math display="inline">\phi_i - \phi_j</math> if nodes <math display="inline">i</math> and <math display="inline">j</math> are connected (if they are not connected, no heat is transferred). Then, for thermal conductivity <math display="inline">k</math>,

:<math>\begin{align}
  \frac{d \phi_i}{d t}
    &= -k \sum_j A_{ij} \left(\phi_i - \phi_j \right) \\
    &= -k \left(\phi_i \sum_j A_{ij} - \sum_j A_{ij} \phi_j \right) \\
    &= -k \left(\phi_i \ \deg(v_i) - \sum_j A_{ij} \phi_j \right) \\
    &= -k \sum_j \left(\delta_{ij} \ \deg(v_i) - A_{ij} \right) \phi_j \\
    &= -k \sum_j \left(L_{ij} \right) \phi_j.
\end{align}</math>

In matrix-vector notation,
:<math>\begin{align}
  \frac{d\phi}{dt} &= -k(D - A)\phi \\
                   &= -kL \phi,
\end{align}</math>

which gives
:<math>\frac{d \phi}{d t} + kL\phi = 0.</math>

Notice that this equation takes the same form as the [[heat equation]], where the matrix −''L'' is replacing the Laplacian operator <math display="inline">\nabla^2</math>; hence, the "graph Laplacian".

To find a solution to this differential equation, apply standard techniques for solving a first-order [[matrix differential equation]].  That is, write <math display="inline">\phi</math> as a linear combination of eigenvectors <math display="inline">\mathbf{v}_i</math> of ''L'' (so that <math display="inline">L\mathbf{v}_i = \lambda_i \mathbf{v}_i</math>) with time-dependent coefficients, <math display="inline">\phi(t) = \sum_i c_i(t) \mathbf{v}_i.</math>

Plugging into the original expression (because ''L'' is a symmetric matrix, its unit-norm eigenvectors <math display="inline">\mathbf{v}_i</math> are orthogonal):

:<math>\begin{align}
              0 ={} &\frac{d\left(\sum_i c_i(t) \mathbf{v}_i\right)}{dt} + kL\left(\sum_i c_i(t) \mathbf{v}_i\right)  \\
              {}={} &\sum_i \left[\frac{dc_i(t)}{dt} \mathbf{v}_i + k c_i(t) L \mathbf{v}_i\right] \\
              {}={} &\sum_i \left[\frac{dc_i(t)}{dt} \mathbf{v}_i + k c_i(t) \lambda_i \mathbf{v}_i\right] \\
  \Rightarrow 0 ={} &\frac{dc_i(t)}{dt} + k \lambda_i c_i(t), \\
\end{align}</math>

whose solution is
:<math>c_i(t) = c_i(0) e^{-k \lambda_i t}.</math>

As shown before, the eigenvalues <math display="inline">\lambda_i</math> of ''L'' are non-negative, showing that the solution to the diffusion equation approaches an equilibrium, because it only exponentially decays or remains constant. This also shows that given <math display="inline">\lambda_i</math> and the initial condition <math display="inline">c_i(0)</math>, the solution at any time ''t'' can be found.<ref>{{cite book
| last                  = Newman
| first                 = Mark
| author-link            = Mark Newman
| title                 = Networks: An Introduction
| year                  = 2010
| publisher             = Oxford University Press
| isbn                  = 978-0199206650
}}</ref>

To find <math display="inline">c_i(0)</math> for each <math display="inline">i</math> in terms of the overall initial condition <math display="inline">\phi(0)</math>, simply project <math display="inline">\phi(0)</math> onto the unit-norm eigenvectors <math display="inline">\mathbf{v}_i</math>;
: <math>c_i(0) = \left\langle \phi(0), \mathbf{v}_i \right\rangle </math>.
This approach has been applied to quantitative heat transfer modelling on unstructured grids.<ref name="Yavari2020">{{cite journal
 | title = Computational heat transfer with spectral graph theory: Quantitative verification 
 | author1=Yavari, R. | author2=Cole, K. D. | author3=Rao, P. K. 
 | journal = International Journal of Thermal Sciences 
 | volume = 153
| pages = 106383
 | year = 2020
| doi=10.1016/j.ijthermalsci.2020.106383
| doi-access = free| bibcode=2020IJTS..15306383C }}</ref>
<ref name="Cole2022">{{cite journal
 | title = Discrete Green's functions and spectral graph theory for computationally efficient thermal modeling
 | author1=Cole, K. D. | author2=Riensche, A. | author3=Rao, P. K. 
 | journal = International Journal of Heat and Mass Transfer
 | volume = 183
| pages = 122112
 | year = 2022
| doi=10.1016/j.ijheatmasstransfer.2021.122112
| s2cid=244652819 | doi-access=free| bibcode=2022IJHMT.18322112C }}</ref>

In the case of undirected graphs, this works because <math display="inline">L</math> is symmetric, and by the [[spectral theorem]], its eigenvectors are all orthogonal.  So the projection onto the eigenvectors of <math display="inline">L</math> is simply an orthogonal coordinate transformation of the initial condition to a set of coordinates which decay exponentially and independently of each other.

=== Equilibrium behavior ===
To understand <math display="inline">\lim_{t \to \infty}\phi(t)</math>, the only terms <math display="inline"> c_i(t) = c_i(0) e^{-k \lambda_i t}</math> that remain are those where <math display="inline">\lambda_i = 0</math>, since
: <math>\lim_{t\to\infty} e^{-k \lambda_i t} = \begin{cases}
  0, & \text{if} & \lambda_i > 0 \\
  1, & \text{if} & \lambda_i = 0
\end{cases}</math>

In other words, the equilibrium state of the system is determined completely by the [[Kernel (linear algebra)|kernel]] of <math display="inline">L</math>.

Since by definition, <math display="inline">\sum_{j}L_{ij} = 0</math>, the vector <math display="inline">\mathbf{v}^1</math> of all ones is in the kernel. If there are <math display="inline">k</math> disjoint [[Connected component (graph theory)|connected components]] in the graph, then this vector of all ones can be split into the sum of <math display="inline">k</math> independent <math display="inline">\lambda = 0</math> eigenvectors of ones and zeros, where each connected component corresponds to an eigenvector with ones at the elements in the connected component and zeros elsewhere.

The consequence of this is that for a given initial condition <math display="inline">\phi(0)</math> for a graph with <math display="inline">N</math> vertices
: <math>\lim_{t\to\infty}\phi(t) = \left\langle \phi(0), \mathbf{v^1} \right\rangle \mathbf{v^1}</math>

where
: <math>\mathbf{v^1} = \frac{1}{\sqrt{N}} [1, 1, \ldots, 1] </math>

For each element <math display="inline">\phi_j</math> of <math display="inline">\phi</math>, i.e. for each vertex <math display="inline">j</math> in the graph, it can be rewritten as
: <math>\lim_{t\to\infty}\phi_j(t) = \frac{1}{N} \sum_{i = 1}^N \phi_i(0) </math>.

In other words, at steady state, the value of <math display="inline">\phi</math> converges to the same value at each of the vertices of the graph, which is the average of the initial values at all of the vertices.  Since this is the solution to the heat diffusion equation, this makes perfect sense intuitively.  We expect that neighboring elements in the graph will exchange energy until that energy is spread out evenly throughout all of the elements that are connected to each other.

===Example of the operator on a grid===
[[File:Graph Laplacian Diffusion Example.gif|thumb|This GIF shows the progression of diffusion, as solved by the graph laplacian technique.  A graph is constructed over a grid, where each pixel in the graph is connected to its 8 bordering pixels.  Values in the image then diffuse smoothly to their neighbors over time via these connections.  This particular image starts off with three strong point values which spill over to their neighbors slowly.  The whole system eventually settles out to the same value at equilibrium.]]

This section shows an example of a function <math display="inline">\phi</math> diffusing over time through a graph.  The graph in this example is constructed on a 2D discrete grid, with points on the grid connected to their eight neighbors.  Three initial points are specified to have a positive value, while the rest of the values in the grid are zero.  Over time, the exponential decay acts to distribute the values at these points evenly throughout the entire grid.

The complete Matlab source code that was used to generate this animation is provided below.  It shows the process of specifying initial conditions, projecting these initial conditions onto the eigenvalues of the Laplacian Matrix, and simulating the exponential decay of these projected initial conditions.

<syntaxhighlight lang="matlab">
N = 20; % The number of pixels along a dimension of the image
A = zeros(N, N); % The image
Adj = zeros(N * N, N * N); % The adjacency matrix

% Use 8 neighbors, and fill in the adjacency matrix
dx = [- 1, 0, 1, - 1, 1, - 1, 0, 1];
dy = [- 1, - 1, - 1, 0, 0, 1, 1, 1];
for x = 1:N
    for y = 1:N
        index = (x - 1) * N + y;
        for ne = 1:length(dx)
            newx = x + dx(ne);
            newy = y + dy(ne);
            if newx > 0 && newx <= N && newy > 0 && newy <= N
                index2 = (newx - 1) * N + newy;
                Adj(index, index2) = 1;
            end
        end
    end
end

% BELOW IS THE KEY CODE THAT COMPUTES THE SOLUTION TO THE DIFFERENTIAL EQUATION
Deg = diag(sum(Adj, 2)); % Compute the degree matrix
L = Deg - Adj; % Compute the laplacian matrix in terms of the degree and adjacency matrices
[V, D] = eig(L); % Compute the eigenvalues/vectors of the laplacian matrix
D = diag(D);

% Initial condition (place a few large positive values around and
% make everything else zero)
C0 = zeros(N, N);
C0(2:5, 2:5) = 5;
C0(10:15, 10:15) = 10;
C0(2:5, 8:13) = 7;
C0 = C0(:);

C0V = V'*C0; % Transform the initial condition into the coordinate system
% of the eigenvectors
for t = 0:0.05:5
    % Loop through times and decay each initial component
    Phi = C0V .* exp(- D * t); % Exponential decay for each component
    Phi = V * Phi; % Transform from eigenvector coordinate system to original coordinate system
    Phi = reshape(Phi, N, N);
    % Display the results and write to GIF file
    imagesc(Phi);
    caxis([0, 10]);
     title(sprintf('Diffusion t = %3f', t));
    frame = getframe(1);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);
    if t == 0
        imwrite(imind, cm, 'out.gif', 'gif', 'Loopcount', inf, 'DelayTime', 0.1);
    else
        imwrite(imind, cm, 'out.gif', 'gif', 'WriteMode', 'append', 'DelayTime', 0.1);
    end
end
</syntaxhighlight>