Python SciPy linalg.LaplacianNd用法及代碼示例

本文簡要介紹 python 語言中 scipy.sparse.linalg.LaplacianNd 的用法。

用法:

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

N 維度中的網格拉普拉斯及其特征值/特征向量。

在均勻矩形網格上構造拉普拉斯算子N維數並輸出其特征值和特征向量。拉普拉斯算子L是平方、負定、實對稱數組，具有帶符號整數條目，否則為零。

參數 ：：

grid_shape： 元組

長度為N(對應於拉帕算子的維度)的整數元組，其中每個條目給出該維度的大小。拉普拉斯矩陣是大小為 np.prod(grid_shape) 的平方。

boundary_conditions： {‘neumann’, ‘dirichlet’, ‘periodic’}，可選

網格邊界上的邊界條件類型。有效值為 'dirichlet' 或 'neumann'``(default) or ``'periodic' 。

dtype： 類型

數組的數值類型。默認為 np.int8 。

注意：

與 1-、2-和 3-D 拉普拉斯算子的 MATLAB/Octave 實現 [1] 相比，此代碼允許任意 N-D 情況和 matrix-free 可調用選項，但目前僅限於純 Dirichlet、Neumann 或僅周期性邊界條件。

矩形網格圖 ( scipy.sparse.csgraph.laplacian ) 的拉普拉斯矩陣對應於具有諾伊曼條件的負拉普拉斯矩陣，即 boundary_conditions = 'neumann' 。

離散拉普拉斯算子的所有特征值和特征向量N形狀的維規則網格grid_shape與網格步長h=1是分析已知的[2]。

參考：

[1]

https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/laplacian/laplacian.m

[2]

“二階導數的特征值和特征向量”，維基百科https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative

例子：

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

下麵演示的具有 n=6 網格點的規則網格上的純諾依曼邊界條件的一維拉普拉斯正好是具有 n 頂點的無向線性圖的負圖拉普拉斯，使用著名的稀疏鄰接矩陣 G 表示tri-diagonal矩陣：

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

由於拉普拉斯算子的所有矩陣項都是整數，因此 'int8' 是用於存儲矩陣表示的默認數據類型。

>>> lap.tosparse()
<6x6 sparse matrix of type '<class 'numpy.int8'>'
    with 16 stored elements (3 diagonals) in DIAgonal format>
>>> 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

可以計算任意數量的極端特征值和/或特征向量。

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

二維拉普拉斯算子在規則網格上進行說明，每個維度都有 grid_shape = (2, 3) 點。

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

網格點計數如下：

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

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

分別示出了邊界條件'dirichlet'、'periodic'和'neumann'；與'dirichlet'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet')
>>> lap.tosparse()
<6x6 sparse array of type '<class 'numpy.int8'>'
    with 20 stored elements in Compressed Sparse Row format>
>>> 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

與'periodic'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
>>> lap.tosparse()
<6x6 sparse array of type '<class 'numpy.int8'>'
    with 24 stored elements in Compressed Sparse Row format>
>>> 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

並使用'neumann'

>>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
>>> lap.tosparse()
<6x6 sparse array of type '<class 'numpy.int8'>'
    with 20 stored elements in Compressed Sparse Row format>
>>> 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

屬性 ：：

H

埃爾米特伴隨。

T

轉置該線性運算符。