Skip to main content
`Ctrl`+`K`
scipy.sparse.linalg.

# LaplacianNd

class scipy.sparse.linalg.LaplacianNd(*args, **kwargs)

  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_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`.

  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]], 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([-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 |