bundles / scipy 1.18.0 / scipy / sparse / linalg / _special_sparse_arrays / LaplacianNd

class

`scipy.sparse.linalg._special_sparse_arrays:LaplacianNd`

source: /scipy/sparse/linalg/_special_sparse_arrays.py:10

Signature

class LaplacianNd(grid_shape, *, boundary_conditions = neumann, dtype = <class 'numpy.int8'>)

Members

Summary

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

Extended Summary

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_shape : tuple: 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'.
dtype : dtype: Numerical type of the array. Default is np.int8.

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

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].

Examples

import numpy as np
from scipy.sparse.linalg import LaplacianNd
from scipy.sparse import diags_array, 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_array(np.ones(n - 1), offsets=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()
np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
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()
lap.eigenvalues()[-2:]
lap.eigenvalues(2)

✓

array([-1.,  0.])

lap.eigenvectors(1)
lap.eigenvectors(2)
lap.eigenvectors()

✗

array([[ 0.40825,  0.5    ,  0.28868,  0.28868,  0.5    ,  0.40825],
       [-0.40825,  0.     , -0.57735,  0.57735,  0.     ,  0.40825],
       [ 0.40825, -0.5    ,  0.28868,  0.28868, -0.5    ,  0.40825],
       [-0.40825,  0.5    ,  0.28868, -0.28868, -0.5    ,  0.40825],
       [ 0.40825,  0.     , -0.57735, -0.57735,  0.     ,  0.40825],
       [-0.40825, -0.5    ,  0.28868, -0.28868,  0.5    ,  0.40825]])

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()
np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
np.array_equal(lap.tosparse().toarray(), lap.toarray())

✓

True

lap.eigenvalues()

✗

array([-6.41421, -5.     , -4.41421, -3.58579, -3.     , -1.58579])

eigvals = eigvalsh(lap.toarray().astype(np.float64))
np.allclose(lap.eigenvalues(), eigvals)
np.allclose(lap.toarray() @ lap.eigenvectors(),
            lap.eigenvectors() @ np.diag(lap.eigenvalues()))

✓

True

with ``'periodic'``

lap = LaplacianNd(grid_shape, boundary_conditions='periodic')

✓

lap.tosparse()
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())
np.array_equal(lap.tosparse().toarray(), lap.toarray())
lap.eigenvalues()
eigvals = eigvalsh(lap.toarray().astype(np.float64))
np.allclose(lap.eigenvalues(), eigvals)
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()
np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
np.array_equal(lap.tosparse().toarray(), lap.toarray())
lap.eigenvalues()
eigvals = eigvalsh(lap.toarray().astype(np.float64))
np.allclose(lap.eigenvalues(), eigvals)
np.allclose(lap.toarray() @ lap.eigenvectors(),
            lap.eigenvectors() @ np.diag(lap.eigenvalues()))

✓

True

Aliases

scipy.sparse.linalg.LaplacianNd