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 sizenp.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
.
- Attributes:
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 theNd
Laplacian corresponding to the m ordered eigenvalues... versionadded:: 1.12.0
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 sizeh=1
are analytically known [2].References
[1][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 withn
vertices using the sparse adjacency matrixG
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