In mathematics, the discrete Laplace operator is an analog of the continuous Laplace operator, defined so that it has meaning on a graph or a discrete grid. For the case of a finite-dimensional graph (having a finite number of edges and vertices), the discrete Laplace operator is more commonly called the Laplacian matrix.
The discrete Laplace operator occurs in physics problems such as the Ising model and loop quantum gravity, as well as in the study of discrete dynamical systems. It is also used in numerical analysis as a stand-in for the continuous Laplace operator. Common applications include image processing,[1] where it is known as the Laplace filter, and in machine learning for clustering and semi-supervised learning on neighborhood graphs.
There are various definitions of the discrete Laplacian for graphs, differing by sign and scale factor (sometimes one averages over the neighboring vertices, other times one just sums; this makes no difference for a regular graph). The traditional definition of the graph Laplacian, given below, corresponds to the negative continuous Laplacian on a domain with a free boundary.
Let
G=(V,E)
V
E
\phi\colonV\toR
\Delta
\phi
(\Delta\phi)(v)=\sumw:d(w,v)=1\left[\phi(v)-\phi(w)\right]
where
d(w,v)
\phi
\Delta\phi
(\Delta\phi)(v)
If the graph has weighted edges, that is, a weighting function
\gamma\colonE\toR
(\Delta\gamma\phi)(v)=\sumw:d(w,v)=1\gammawv\left[\phi(v)-\phi(w)\right]
where
\gammawv
wv\inE
Closely related to the discrete Laplacian is the averaging operator:
(M\phi)(v)= | 1 |
\degv |
\sumw:d(w,v)=1\phi(w).
In addition to considering the connectivity of nodes and edges in a graph, mesh Laplace operators take into account the geometry of a surface (e.g. the angles at the nodes). For a two-dimensional manifold triangle mesh, the Laplace-Beltrami operator of a scalar function
u
i
(\Deltau)i\equiv
1 | |
2Ai |
\sumj(\cot\alphaij+\cot\betaij)(uj-ui),
where the sum is over all adjacent vertices
j
i
\alphaij
\betaij
ij
Ai
i
i
To facilitate computation, the Laplacian is encoded in a matrix
L\inR|V| x |V|
Lu=(\Deltau)i
C
Cij=\begin{cases}
1 | |
2 |
(\cot\alphaij+\cot\betaij)&ijisanedge,thatisj\inN(i),\\ -\sum\limitskCik&i=j,\\ 0&otherwise \end{cases}
Where
N(i)
i
And let
M
M
i
Ai
L=M-1C
A more general overview of mesh operators is given in.[3]
Approximations of the Laplacian, obtained by the finite-difference method or by the finite-element method, can also be called discrete Laplacians. For example, the Laplacian in two dimensions can be approximated using the five-point stencil finite-difference method, resulting in
\Deltaf(x,y) ≈
f(x-h,y)+f(x+h,y)+f(x,y-h)+f(x,y+h)-4f(x,y) | |
h2 |
,
where the grid size is h in both dimensions, so that the five-point stencil of a point (x, y) in the grid is
\{(x-h,y),(x,y),(x+h,y),(x,y-h),(x,y+h)\}.
If the grid size h = 1, the result is the negative discrete Laplacian on the graph, which is the square lattice grid. There are no constraints here on the values of the function f(x, y) on the boundary of the lattice grid, thus this is the case of no source at the boundary, that is, a no-flux boundary condition (aka, insulation, or homogeneous Neumann boundary condition). The control of the state variable at the boundary, as f(x, y) given on the boundary of the grid (aka, Dirichlet boundary condition), is rarely used for graph Laplacians, but is common in other applications.
Multidimensional discrete Laplacians on rectangular cuboid regular grids have very special properties, e.g., they are Kronecker sums of one-dimensional discrete Laplacians, see Kronecker sum of discrete Laplacians, in which case all its eigenvalues and eigenvectors can be explicitly calculated.
In this approach, the domain is discretized into smaller elements, often triangles or tetrahedra, but other elements such as rectangles or cuboids are possible. The solution space is then approximated using so called form-functions of a pre-defined degree. The differential equation containing the Laplace operator is then transformed into a variational formulation, and a system of equations is constructed (linear or eigenvalue problems). The resulting matrices are usually very sparse and can be solved with iterative methods.
Discrete Laplace operator is often used in image processing e.g. in edge detection and motion estimation applications.[4] The discrete Laplacian is defined as the sum of the second derivatives Laplace operator#Coordinate expressions and calculated as sum of differences over the nearest neighbours of the central pixel. Since derivative filters are often sensitive to noise in an image, the Laplace operator is often preceded by a smoothing filter (such as a Gaussian filter) in order to remove the noise before calculating the derivative. The smoothing filter and Laplace filter are often combined into a single filter.[5]
For one-, two- and three-dimensional signals, the discrete Laplacian can be given as convolution with the following kernels:
1D filter:
2 | |
\vec{D} | |
x=\begin{bmatrix}1 |
&-2&1\end{bmatrix}
2D filter:
2 | |
D | |
xy |
=\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}
2 | |
D | |
xy |
2D filter:
2 | |
D | |
xy |
=\begin{bmatrix}0.25&0.5&0.25\\0.5&-3&0.5\\0.25&0.5&0.25\end{bmatrix}
3D filter:
2 | |
D | |
xyz |
first plane =
\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\end{bmatrix}
\begin{bmatrix}0&1&0\\1&-6&1\\0&1&0\end{bmatrix}
\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\end{bmatrix}
and using 27-point stencil by:[7]
first plane =
1 | |
26 |
\begin{bmatrix}2&3&2\\3&6&3\\2&3&2\end{bmatrix}
1 | |
26 |
\begin{bmatrix}3&6&3\\6&-88&6\\3&6&3\end{bmatrix}
1 | |
26 |
\begin{bmatrix}2&3&2\\3&6&3\\2&3&2\end{bmatrix}
D filter: For the element
a | |
x1,x2,...,xn |
2 | |
D | |
x1,x2,...,xn |
,
a | |
x1,x2,...,xn |
=\left\{\begin{array}{ll} -2n&ifs=n,\\ 1&ifs=n-1,\\ 0&otherwise, \end{array}\right.
where is the position (either, or) of the element in the kernel in the -th direction, and is the number of directions for which .
Note that the nD version, which is based on the graph generalization of the Laplacian, assumes all neighbors to be at an equal distance, and hence leads to the following 2D filter with diagonals included, rather than the version above:
2D filter:
2 | |
D | |
xy |
=\begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix}.
These kernels are deduced by using discrete differential quotients.
It can be shown[8] [9] that the following discrete approximation of the two-dimensional Laplacian operator as a convex combination of difference operators
2 | |
\nabla | |
\gamma |
=(1-\gamma)
2 | |
\nabla | |
5 |
+\gamma\nabla
2 | |
x |
=(1-\gamma)\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix} +\gamma\begin{bmatrix}1/2&0&1/2\\0&-2&0\\1/2&0&1/2\end{bmatrix}
for γ ∈ [0, 1] is compatible with discrete scale-space properties, where specifically the value γ = 1/3 gives the best approximation of rotational symmetry.[8] [9] [10] Regarding three-dimensional signals, it is shown[9] that the Laplacian operator can be approximated by the two-parameter family of difference operators
2 | |
\nabla | |
\gamma1,\gamma2 |
=(1-\gamma1-\gamma2)
2 | |
\nabla | |
7 |
+\gamma1
2 | |
\nabla | |
+3 |
+\gamma2
2 | |
\nabla | |
x 3 |
),
where
2 | |
(\nabla | |
7 |
f)0,= f-1,+f+1,+f0,+f0,+f0,+f0,-6f0,,
2 | |
(\nabla | |
+3 |
f)0,=
1 | |
4 |
(f-1,+f-1,+f+1,+f+1,+f-1,+f-1,+f+1,+f+1,+f0,+f0,+f0,+f0,-12f0,),
2 | |
(\nabla | |
x 3 |
f)0,=
1 | |
4 |
(f-1,+f-1,+f-1,+f-1,+f+1,+f+1,+f+1,+f+1,-8f0,).
A discrete signal, comprising images, can be viewed as a discrete representation of a continuous function
f(\barr)
\barr\inRn
f\inR
f
f(\barr)=\sumk\infk\muk(\barr)
fk\inR
f
K
\muk
K
\muk(\barr)=\mu(\barr-\barrk)
\mu
n
\barr=(x1,x2...x
T | |
n) |
\mu
n
f(\barr)
\nabla2f(\barrk)=\sumk'\infk'(\nabla2\mu(\barr-\barrk'
))| | |
\barr=\barrk |
K
f
fk
n
\barr
n
The spectrum of the discrete Laplacian on an infinite grid is of key interest; since it is a self-adjoint operator, it has a real spectrum. For the convention
\Delta=I-M
Z
[0,2]
[-1,1]
If the graph is an infinite square lattice grid, then this definition of the Laplacian can be shown to correspond to the continuous Laplacian in the limit of an infinitely fine grid. Thus, for example, on a one-dimensional grid we have
\partial2F | |
\partialx2 |
=\lim\epsilon
[F(x+\epsilon)-F(x)]-[F(x)-F(x-\epsilon)] | |
\epsilon2 |
.
This definition of the Laplacian is commonly used in numerical analysis and in image processing. In image processing, it is considered to be a type of digital filter, more specifically an edge filter, called the Laplace filter.
Suppose describes a temperature distribution across a graph, where is the temperature at vertex . According to Newton's law of cooling, the heat transferred from node to node is proportional to if nodes and are connected (if they are not connected, no heat is transferred). Then, for thermal conductivity ,
\begin{align}
d\phii | |
dt |
&=-k\sumjAij\left(\phii-\phij\right)\\ &=-k\left(\phii\sumjAij-\sumjAij\phij\right)\\ &=-k\left(\phii \deg(vi)-\sumjAij\phij\right)\\ &=-k\sumj\left(\deltaij \deg(vi)-Aij\right)\phij\\ &=-k\sumj\left(Lij\right)\phij. \end{align}
In matrix-vector notation,
\begin{align}
d\phi | |
dt |
&=-k(D-A)\phi\\ &=-kL\phi, \end{align}
which gives
d\phi | |
dt |
+kL\phi=0.
Notice that this equation takes the same form as the heat equation, where the matrix −L is replacing the Laplacian operator ; 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 as a linear combination of eigenvectors of L (so that ) with time-dependent coefficients,
Plugging into the original expression (because L is a symmetric matrix, its unit-norm eigenvectors are orthogonal):
\begin{align} 0={}&
d\left(\sumici(t)vi\right) | |
dt |
+kL\left(\sumici(t)vi\right)\\ {}={}&\sumi\left[
dci(t) | |
dt |
vi+kci(t)Lvi\right]\\ {}={}&\sumi\left[
dci(t) | |
dt |
vi+kci(t)λivi\right]\\ ⇒ 0={}&
dci(t) | |
dt |
+kλici(t),\\ \end{align}
whose solution is
ci(t)=ci(0)
-kλit | |
e |
.
As shown before, the eigenvalues 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 and the initial condition , the solution at any time t can be found.[12]
To find for each in terms of the overall initial condition , simply project onto the unit-norm eigenvectors ;
ci(0)=\left\langle\phi(0),vi\right\rangle
In the case of undirected graphs, this works because is symmetric, and by the spectral theorem, its eigenvectors are all orthogonal. So the projection onto the eigenvectors of is simply an orthogonal coordinate transformation of the initial condition to a set of coordinates which decay exponentially and independently of each other.
To understand , the only terms that remain are those where , since
\limt\toinfty
-kλit | |
e |
=\begin{cases} 0,&if&λi>0\\ 1,&if&λi=0 \end{cases}
In other words, the equilibrium state of the system is determined completely by the kernel of .
Since by definition, , the vector of all ones is in the kernel. If there are disjoint connected components in the graph, then this vector of all ones can be split into the sum of independent 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 for a graph with vertices
\limt\toinfty\phi(t)=\left\langlec(0),
v1 |
\right\rangle
v1 |
where
v1 |
=
1 | |
\sqrt{N |
For each element of , i.e. for each vertex in the graph, it can be rewritten as
\limt\toinfty\phij(t)=
1 | |
N |
N | |
\sum | |
i=1 |
ci(0)
In other words, at steady state, the value of 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.
This section shows an example of a function 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.
% Use 8 neighbors, and fill in the adjacency matrixdx = [- 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 endend
% BELOW IS THE KEY CODE THAT COMPUTES THE SOLUTION TO THE DIFFERENTIAL EQUATIONDeg = diag(sum(Adj, 2)); % Compute the degree matrixL = 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 matrixD = 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 eigenvectorsfor 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
Let
P\colonV → R
\phi
(P\phi)(v)=P(v)\phi(v).
Then
H=\Delta+P
If the number of edges meeting at a vertex is uniformly bounded, and the potential is bounded, then H is bounded and self-adjoint.
The spectral properties of this Hamiltonian can be studied with Stone's theorem; this is a consequence of the duality between posets and Boolean algebras.
On regular lattices, the operator typically has both traveling-wave as well as Anderson localization solutions, depending on whether the potential is periodic or random.
The Green's function of the discrete Schrödinger operator is given in the resolvent formalism by
G(v,w;λ)=\left\langle\deltav\left|
1 | |
H-λ |
\right|\deltaw\right\rangle
\deltaw
\deltaw(v)=\deltawv
For fixed
w\inV
λ
(H-λ)G(v,w;λ)=\deltaw(v).
Certain equations involving the discrete Laplacian only have solutions on the simply-laced Dynkin diagrams (all edges multiplicity 1), and are an example of the ADE classification. Specifically, the only positive solutions to the homogeneous equation:
\Delta\phi=\phi,
"Twice any label is the sum of the labels on adjacent vertices,"are on the extended (affine) ADE Dynkin diagrams, of which there are 2 infinite families (A and D) and 3 exceptions (E). The resulting numbering is unique up to scale, and if the smallest value is set at 1, the other numbers are integers, ranging up to 6.
The ordinary ADE graphs are the only graphs that admit a positive labeling with the following property:
Twice any label minus two is the sum of the labels on adjacent vertices.In terms of the Laplacian, the positive solutions to the inhomogeneous equation:
\Delta\phi=\phi-2.