Discrete Laplace operator

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.

For the discrete equivalent of the Laplace transform, see Z-transform.

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.

Definitions

Graph Laplacians

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 ${\displaystyle G=(V,E)}$ be a graph with vertices ${\displaystyle V}$ and edges ${\displaystyle E}$. Let ${\displaystyle \phi \colon V\to R}$ be a function of the vertices taking values in a ring. Then, the discrete Laplacian ${\displaystyle \Delta }$ acting on ${\displaystyle \phi }$ is defined by

${\displaystyle (\Delta \phi )(v)=\sum _{w:\,d(w,v)=1}\left[\phi (v)-\phi (w)\right]}$

where ${\displaystyle d(w,v)}$ is the graph distance between vertices w and v. Thus, this sum is over the nearest neighbors of the vertex v. For a graph with a finite number of edges and vertices, this definition is identical to that of the Laplacian matrix. That is, ${\displaystyle \phi }$ can be written as a column vector; and so ${\displaystyle \Delta \phi }$ is the product of the column vector and the Laplacian matrix, while ${\displaystyle (\Delta \phi )(v)}$ is just the v'th entry of the product vector.

If the graph has weighted edges, that is, a weighting function ${\displaystyle \gamma \colon E\to R}$ is given, then the definition can be generalized to

${\displaystyle (\Delta _{\gamma }\phi )(v)=\sum _{w:\,d(w,v)=1}\gamma _{wv}\left[\phi (v)-\phi (w)\right]}$

where ${\displaystyle \gamma _{wv}}$ is the weight value on the edge ${\displaystyle wv\in E}$.

Closely related to the discrete Laplacian is the averaging operator:

${\displaystyle (M\phi )(v)={\frac {1}{\deg v}}\sum _{w:\,d(w,v)=1}\phi (w).}$

Mesh Laplacians

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 manifold triangle mesh, the Laplace-Beltrami operator of a scalar function ${\displaystyle u}$ at a vertex ${\displaystyle i}$ can be approximated as

${\displaystyle (\Delta u)_{i}\approx {\frac {1}{2A_{i}}}\Sigma _{j}(\cot \alpha _{ij}+\cot \beta _{ij})(u_{i}-u_{j}),}$

where the sum is over all adjacent vertices ${\displaystyle j}$ of ${\displaystyle i}$, ${\displaystyle \alpha _{ij}}$ and ${\displaystyle \beta _{ij}}$ are the two angles opposite of the edge ${\displaystyle ij}$, and ${\displaystyle A_{i}}$ is the vertex area of ${\displaystyle i}$; that is, one third of the summed areas of triangles incident to ${\displaystyle i}$. The above cotangent formula can be derived using many different methods among which are piecewise linear finite elements, finite volumes (see [2] for a derivation), and discrete exterior calculus (see ).

To facilitate computation, the Laplacian is encoded in a matrix ${\displaystyle L\in \mathbb {R} ^{|V|\times |V|}}$ such that ${\displaystyle Lu=(\Delta u)_{i}}$. Let ${\displaystyle C}$ be the (sparse) cotangent matrix with entries

${\displaystyle C_{ij}={\begin{cases}-{\frac {1}{2}}(\cot \alpha _{ij}+\cot \beta _{ij})&ij{\text{ is an edge}},\\-\sum \limits _{k\in N(i)}C_{ik}&i=j,\\0&{\text{otherwise,}}\end{cases}}}$ Where ${\displaystyle N(i)}$ denotes the neighborhood of ${\displaystyle i}$.

And let ${\displaystyle M}$ be the diagonal mass matrix ${\displaystyle M}$ whose ${\displaystyle i}$-th entry along the diagonal is the vertex area ${\displaystyle A_{i}}$. Then ${\displaystyle L=M^{-1}C}$ is the sought discretization of the Laplacian.

A more general overview of mesh operators is given in.[3]

Finite differences

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

${\displaystyle \Delta f(x,y)\approx {\frac {f(x-h,y)+f(x+h,y)+f(x,y-h)+f(x,y+h)-4f(x,y)}{h^{2}}},}$

where the grid size is h in both dimensions, so that the five-point stencil of a point (x, y) in the grid is

${\displaystyle \{(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.

Finite-element method

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.

Image processing

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]

Implementation via operator discretization

For one-, two- and three-dimensional signals, the discrete Laplacian can be given as convolution with the following kernels:

1D filter: ${\displaystyle {\vec {D}}_{x}^{2}={\begin{bmatrix}1&-2&1\end{bmatrix}}}$,

2D filter: ${\displaystyle \mathbf {D} _{xy}^{2}={\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}}}$.

${\displaystyle \mathbf {D} _{xy}^{2}}$ corresponds to the (Five-point stencil) finite-difference formula seen previously. It is stable for very smoothly varying fields, but for equations with rapidly varying solutions more stable and isotropic form of the laplacian operator is required,[6] such as the nine-point stencil, which includes the diagonals:

2D filter: ${\displaystyle \mathbf {D} _{xy}^{2}={\begin{bmatrix}0.25&0.5&0.25\\0.5&-3&0.5\\0.25&0.5&0.25\end{bmatrix}}}$,

3D filter: ${\displaystyle \mathbf {D} _{xyz}^{2}}$ using seven-point stencil is given by:

first plane = ${\displaystyle {\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\end{bmatrix}}}$; second plane = ${\displaystyle {\begin{bmatrix}0&1&0\\1&-6&1\\0&1&0\end{bmatrix}}}$; third plane = ${\displaystyle {\begin{bmatrix}0&0&0\\0&1&0\\0&0&0\end{bmatrix}}}$.

and using 27-point stencil by:[7]

first plane = ${\displaystyle {\frac {1}{26}}{\begin{bmatrix}2&3&2\\3&6&3\\2&3&2\end{bmatrix}}}$; second plane = ${\displaystyle {\frac {1}{26}}{\begin{bmatrix}3&6&3\\6&-88&6\\3&6&3\end{bmatrix}}}$; third plane = ${\displaystyle {\frac {1}{26}}{\begin{bmatrix}2&3&2\\3&6&3\\2&3&2\end{bmatrix}}}$.

nD filter: For the element ${\displaystyle a_{x_{1},x_{2},\dots ,x_{n}}}$ of the kernel ${\displaystyle \mathbf {D} _{x_{1},x_{2},\dots ,x_{n}}^{2},}$

${\displaystyle a_{x_{1},x_{2},\dots ,x_{n}}=\left\{{\begin{array}{ll}-2n&{\text{if }}s=n,\\1&{\text{if }}s=n-1,\\0&{\text{otherwise,}}\end{array}}\right.}$

where x_i is the position (either −1, 0 or 1) of the element in the kernel in the i-th direction, and s is the number of directions i for which x_i = 0.

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: ${\displaystyle \mathbf {D} _{xy}^{2}={\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

${\displaystyle \nabla _{\gamma }^{2}=(1-\gamma )\nabla _{5}^{2}+\gamma \nabla _{\times }^{2}=(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.[10] Regarding three-dimensional signals, it is shown[9] that the Laplacian operator can be approximated by the two-parameter family of difference operators

${\displaystyle \nabla _{\gamma _{1},\gamma _{2}}^{2}=(1-\gamma _{1}-\gamma _{2})\,\nabla _{7}^{2}+\gamma _{1}\,\nabla _{+^{3}}^{2}+\gamma _{2}\,\nabla _{\times ^{3}}^{2}),}$

where

${\displaystyle (\nabla _{7}^{2}f)_{0,0,0}=f_{-1,0,0}+f_{+1,0,0}+f_{0,-1,0}+f_{0,+1,0}+f_{0,0,-1}+f_{0,0,+1}-6f_{0,0,0},}$

${\displaystyle (\nabla _{+^{3}}^{2}f)_{0,0,0}={\frac {1}{4}}(f_{-1,-1,0}+f_{-1,+1,0}+f_{+1,-1,0}+f_{+1,+1,0}+f_{-1,0,-1}+f_{-1,0,+1}+f_{+1,0,-1}+f_{+1,0,+1}+f_{0,-1,-1}+f_{0,-1,+1}+f_{0,+1,-1}+f_{0,+1,+1}-12f_{0,0,0}),}$

${\displaystyle (\nabla _{\times ^{3}}^{2}f)_{0,0,0}={\frac {1}{4}}(f_{-1,-1,-1}+f_{-1,-1,+1}+f_{-1,+1,-1}+f_{-1,+1,+1}+f_{+1,-1,-1}+f_{+1,-1,+1}+f_{+1,+1,-1}+f_{+1,+1,+1}-8f_{0,0,0}).}$

Implementation via continuous reconstruction

A discrete signal, comprising images, can be viewed as a discrete representation of a continuous function ${\displaystyle f({\bar {r}})}$, where the coordinate vector ${\displaystyle {\bar {r}}\in R^{n}}$ and the value domain is real ${\displaystyle f\in R}$. Derivation operation is therefore directly applicable to the continuous function, ${\displaystyle f}$. In particular any discrete image, with reasonable presumptions on the discretization process, e.g. assuming band limited functions, or wavelets expandable functions, etc. can be reconstructed by means of well-behaving interpolation functions underlying the reconstruction formulation,[11]

${\displaystyle f({\bar {r}})=\sum _{k\in K}f_{k}\mu _{k}({\bar {r}})}$

where ${\displaystyle f_{k}\in R}$ are discrete representations of ${\displaystyle f}$ on grid ${\displaystyle K}$ and ${\displaystyle \mu _{k}}$ are interpolation functions specific to the grid ${\displaystyle K}$. On a uniform grid, such as images, and for bandlimited functions, interpolation functions are shift invariant amounting to ${\displaystyle \mu _{k}({\bar {r}})=\mu ({\bar {r}}-{\bar {r}}_{k})}$ with ${\displaystyle \mu }$ being an appropriately dilated sinc function defined in ${\displaystyle n}$-dimensions i.e. ${\displaystyle {\bar {r}}=(x_{1},x_{2}...x_{n})^{T}}$. Other approximations of ${\displaystyle \mu }$ on uniform grids, are appropriately dilated Gaussian functions in ${\displaystyle n}$-dimensions. Accordingly the discrete Laplacian becomes a discrete version of the Laplacian of the continuous ${\displaystyle f({\bar {r}})}$

${\displaystyle \nabla ^{2}f({\bar {r}}_{k})=\sum _{k'\in K}f_{k'}(\nabla ^{2}\mu ({\bar {r}}-{\bar {r}}_{k'}))|_{{\bar {r}}={\bar {r}}_{k}}}$

which in turn is a convolution with the Laplacian of the interpolation function on the uniform (image) grid ${\displaystyle K}$. An advantage of using Gaussians as interpolation functions is that they yield linear operators, including Laplacians, that are free from rotational artifacts of the coordinate frame in which ${\displaystyle f}$ is represented via ${\displaystyle f_{k}}$, in ${\displaystyle n}$-dimensions, and are frequency aware by definition. A linear operator has not only a limited range in the ${\displaystyle {\bar {r}}}$ domain but also an effective range in the frequency domain (alternatively Gaussian scale space) which can be controlled explicitly via the variance of the Gaussian in a principled manner. The resulting filtering can be implemented by separable filters and decimation (signal processing)/pyramid (image processing) representations for further computational efficiency in ${\displaystyle n}$-dimensions. In other words, the discrete Laplacian filter of any size can be generated conveniently as the sampled Laplacian of Gaussian with spatial size befitting the needs of a particular application as controlled by its variance. Monomials which are non-linear operators can also be implemented using a similar reconstruction and approximation approach provided that the signal is sufficiently over-sampled. Thereby, such non-linear operators e.g. Structure Tensor, and Generalized Structure Tensor which are used in pattern recognition for their total least-square optimality in orientation estimation, can be realized.

Spectrum

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 ${\displaystyle \Delta =I-M}$ on ${\displaystyle Z}$, the spectrum lies within ${\displaystyle [0,2]}$ (as the averaging operator has spectral values in ${\displaystyle [-1,1]}$). This may also be seen by applying the Fourier transform. Note that the discrete Laplacian on an infinite grid has purely absolutely continuous spectrum, and therefore, no eigenvalues or eigenfunctions.

Theorems

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

${\displaystyle {\frac {\partial ^{2}F}{\partial x^{2}}}=\lim _{\epsilon \rightarrow 0}{\frac {[F(x+\epsilon )-F(x)]+[F(x-\epsilon )-F(x)]}{\epsilon ^{2}}}.}$

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.

Discrete Schrödinger operator

Let ${\displaystyle P\colon V\rightarrow R}$ be a potential function defined on the graph. Note that P can be considered to be a multiplicative operator acting diagonally on ${\displaystyle \phi }$

${\displaystyle (P\phi )(v)=P(v)\phi (v).}$

Then ${\displaystyle H=\Delta +P}$ is the discrete Schrödinger operator, an analog of the continuous Schrödinger operator.

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.

Discrete Green's function

The Green's function of the discrete Schrödinger operator is given in the resolvent formalism by

${\displaystyle G(v,w;\lambda )=\left\langle \delta _{v}\left|{\frac {1}{H-\lambda }}\right|\delta _{w}\right\rangle }$

where ${\displaystyle \delta _{w}}$ is understood to be the Kronecker delta function on the graph: ${\displaystyle \delta _{w}(v)=\delta _{wv}}$; that is, it equals 1 if v=w and 0 otherwise.

For fixed ${\displaystyle w\in V}$ and ${\displaystyle \lambda }$ a complex number, the Green's function considered to be a function of v is the unique solution to

${\displaystyle (H-\lambda )G(v,w;\lambda )=\delta _{w}(v).}$

ADE classification

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:

${\displaystyle \Delta \phi =\phi ,}$

in words,

"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:

${\displaystyle \Delta \phi =\phi -2.}$

The resulting numbering is unique (scale is specified by the "2"), and consists of integers; for E₈ they range from 58 to 270, and have been observed as early as 1968.[12]

See also

• Spectral shape analysis
• Electrical network
• Kronecker sum of discrete Laplacians
• Discrete calculus

References

1. Leventhal, Daniel (Autumn 2011). "Image processing" (PDF). University of Washington. Retrieved 2019-12-01.
2. Crane, K., de Goes, F., Desbrun, M., and Schröder, P. (2013). "Digital geometry processing with discrete exterior calculus". ACM SIGGRAPH 2013 Courses. SIGGRAPH '13. 7. ACM, New York, NY, USA. pp. 7:1–7:126. doi:10.1145/2504435.2504442.CS1 maint: multiple names: authors list (link)
3. Reuter, M. and Biasotti, S. and Giorgi, D. and Patane, G. and Spagnuolo, M." (2009). "Discrete Laplace-Beltrami operators for shape analysis and segmentation". Computers & Graphics. 33 (3): 381–390df. CiteSeerX 10.1.1.157.757. doi:10.1016/j.cag.2009.03.005.CS1 maint: multiple names: authors list (link)
4. Forsyth, D. A.; Ponce, J. (2003). "Computer Vision". Computers & Graphics. 33 (3): 381–390. CiteSeerX 10.1.1.157.757. doi:10.1016/j.cag.2009.03.005.
5. Matthys, Don (Feb 14, 2001). "LoG Filter". Marquette University. Retrieved 2019-12-01.
6. Provatas, Nikolas; Elder, Ken (2010-10-13). Phase-Field Methods in Materials Science and Engineering (PDF). Weinheim, Germany: Wiley-VCH Verlag GmbH & Co. KGaA. p. 219. doi:10.1002/9783527631520. ISBN 978-3-527-63152-0.
7. O’Reilly, H.; Beck, Jeffrey M. (2006). "A Family of Large-Stencil Discrete Laplacian Approximations in Three Dimensions". International Journal For Numerical Methods in Engineering: 1–16.
8. Lindeberg, T., "Scale-space for discrete signals", PAMI(12), No. 3, March 1990, pp. 234–254.
9. Lindeberg, T., Scale-Space Theory in Computer Vision, Kluwer Academic Publishers, 1994, ISBN 0-7923-9418-6.
10. Patra, Michael; Karttunen, Mikko (2006). "Stencils with isotropic discretization error for differential operators". Numerical Methods for Partial Differential Equations. 22 (4): 936–953. doi:10.1002/num.20129. ISSN 0749-159X.
11. Bigun, J. (2006). Vision with Direction. Springer. doi:10.1007/b138918. ISBN 978-3-540-27322-6.
12. Bourbaki, Nicolas (1968), "Chapters 4–6", Groupes et algebres de Lie, Paris: Hermann

• T.Sunada, Discrete geometric analysis, Proceedings of Symposia in Pure Mathematics (ed. by P. Exner, J. P. Keating, P. Kuchment, T. Sunada, A. Teplyaev), 77(2008), 51-86

External links

• Definition and application to spectral gap