scipy.sparse.linalg.

LaplacianNd#

class scipy.sparse.linalg.LaplacianNd(*args, **kwargs)[source]#

The grid Laplacian in N dimensions and its eigenvalues/eigenvectors.

Construct Laplacian on a uniform rectangular grid in N dimensions and output its eigenvalues and eigenvectors. The Laplacian L is square, negative definite, real symmetric array with signed integer entries and zeros otherwise.

Parameters:

grid_shapetuple: A tuple of integers of length N (corresponding to the dimension of the Lapacian), where each entry gives the size of that dimension. The Laplacian matrix is square of the size np.prod(grid_shape).
boundary_conditions{‘neumann’, ‘dirichlet’, ‘periodic’}, optional: The type of the boundary conditions on the boundaries of the grid. Valid values are 'dirichlet' or 'neumann'``(default) or ``'periodic'.
dtypedtype: Numerical type of the array. Default is np.int8.

Notes

Compared to the MATLAB/Octave implementation [1] of 1-, 2-, and 3-D Laplacian, this code allows the arbitrary N-D case and the matrix-free callable option, but is currently limited to pure Dirichlet, Neumann or Periodic boundary conditions only.

The Laplacian matrix of a graph (scipy.sparse.csgraph.laplacian) of a rectangular grid corresponds to the negative Laplacian with the Neumann conditions, i.e., boundary_conditions = 'neumann'.

All eigenvalues and eigenvectors of the discrete Laplacian operator for an N-dimensional regular grid of shape grid_shape with the grid step size h=1 are analytically known [2].

References

[1]

lobpcg/blopex

[2]

“Eigenvalues and eigenvectors of the second derivative”, Wikipedia https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative

Examples

>>> import numpy as np
>>> from scipy.sparse.linalg import LaplacianNd
>>> from scipy.sparse import diags, csgraph
>>> from scipy.linalg import eigvalsh

The one-dimensional Laplacian demonstrated below for pure Neumann boundary conditions on a regular grid with n=6 grid points is exactly the negative graph Laplacian for the undirected linear graph with n vertices using the sparse adjacency matrix G represented by the famous tri-diagonal matrix:

>>> n = 6
>>> G = diags(np.ones(n - 1), 1, format='csr')
>>> Lf = csgraph.laplacian(G, symmetrized=True, form='function')
>>> grid_shape = (n, )
>>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
>>> np.array_equal(lap.matmat(np.eye(n)), -Lf(np.eye(n)))
True

Since all matrix entries of the Laplacian are integers, 'int8' is the default dtype for storing matrix representations.

>>> lap.tosparse()
<DIAgonal sparse array of dtype 'int8'
    with 16 stored elements (3 diagonals) and shape (6, 6)>
>>> lap.toarray()
array([[-1,  1,  0,  0,  0,  0],
       [ 1, -2,  1,  0,  0,  0],
       [ 0,  1, -2,  1,  0,  0],
       [ 0,  0,  1, -2,  1,  0],
       [ 0,  0,  0,  1, -2,  1],
       [ 0,  0,  0,  0,  1, -1]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True

Any number of extreme eigenvalues and/or eigenvectors can be computed.

>>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
>>> lap.eigenvalues()
array([-4., -3., -3., -1., -1.,  0.])
>>> lap.eigenvalues()[-2:]
array([-1.,  0.])
>>> lap.eigenvalues(2)
array([-1.,  0.])
>>> lap.eigenvectors(1)
array([[0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829]])
>>> lap.eigenvectors(2)
array([[ 0.5       ,  0.40824829],
       [ 0.        ,  0.40824829],
       [-0.5       ,  0.40824829],
       [-0.5       ,  0.40824829],
       [ 0.        ,  0.40824829],
       [ 0.5       ,  0.40824829]])
>>> lap.eigenvectors()
array([[ 0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,
         0.40824829],
       [-0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,
         0.40824829],
       [ 0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,
         0.40824829],
       [-0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,
         0.40824829],
       [ 0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,
         0.40824829],
       [-0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,
         0.40824829]])

The two-dimensional Laplacian is illustrated on a regular grid with grid_shape = (2, 3) points in each dimension.

>>> grid_shape = (2, 3)
>>> n = np.prod(grid_shape)

Numeration of grid points is as follows:

>>> np.arange(n).reshape(grid_shape + (-1,))
array([[[0],
        [1],
        [2]],

       [[3],
        [4],
        [5]]])

Each of the boundary conditions 'dirichlet', 'periodic', and 'neumann' is illustrated separately; with 'dirichlet'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 20 stored elements and shape (6, 6)>
>>> lap.toarray()
array([[-4,  1,  0,  1,  0,  0],
       [ 1, -4,  1,  0,  1,  0],
       [ 0,  1, -4,  0,  0,  1],
       [ 1,  0,  0, -4,  1,  0],
       [ 0,  1,  0,  1, -4,  1],
       [ 0,  0,  1,  0,  1, -4]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-6.41421356, -5.        , -4.41421356, -3.58578644, -3.        ,
       -1.58578644])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

with 'periodic'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 24 stored elements and shape (6, 6)>
>>> lap.toarray()
    array([[-4,  1,  1,  2,  0,  0],
           [ 1, -4,  1,  0,  2,  0],
           [ 1,  1, -4,  0,  0,  2],
           [ 2,  0,  0, -4,  1,  1],
           [ 0,  2,  0,  1, -4,  1],
           [ 0,  0,  2,  1,  1, -4]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-7., -7., -4., -3., -3.,  0.])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

and with 'neumann'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 20 stored elements and shape (6, 6)>
>>> lap.toarray()
array([[-2,  1,  0,  1,  0,  0],
       [ 1, -3,  1,  0,  1,  0],
       [ 0,  1, -2,  0,  0,  1],
       [ 1,  0,  0, -2,  1,  0],
       [ 0,  1,  0,  1, -3,  1],
       [ 0,  0,  1,  0,  1, -2]])
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-5., -3., -3., -2., -1.,  0.])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

Attributes:

H: Hermitian adjoint.
T: Transpose this linear operator.

Methods

toarray()	Construct a dense array from Laplacian data
tosparse()	Construct a sparse array from Laplacian data
eigenvalues(m=None)	Construct a 1D array of m largest (smallest in absolute value) eigenvalues of the Laplacian matrix in ascending order.
eigenvectors(m=None):	Construct the array with columns made of m eigenvectors (`float`) of the `Nd` Laplacian corresponding to the m ordered eigenvalues.
.. versionadded:: 1.12.0