使用 SciPy 模块构建稀疏的 2D 拉普拉斯矩阵答案

【问题标题】：Buildin a sparse 2D laplacian matrix using SciPy modules使用 SciPy 模块构建稀疏的 2D 拉普拉斯矩阵
【发布时间】：2019-05-13 11:53:45
【问题描述】：

我需要构造如下所示的 2D 拉普拉斯算子：

，其中

，I 是单位矩阵。到目前为止，我已经使用diags method of scipy 完成了它，但我想知道是否有更聪明的方法可以使用block_diag method 来完成它。有人试过用这种方法构建二维拉普拉斯算子吗？

我目前的创建方法是通过这个函数：

from scipy.sparse import diags
# Defining the size of the matrix
nx = 3
ny = 3
N  = nx*ny
main_diag = [-4.0 for i in xrange(N)]
side_diag = []
for i in xrange(1,N):
    if i%4 == 0:
        side_diag.append(0)
    else:
        side_diag.append(1)
up_down_diag = [1.0 for i in xrange(N-4)]
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

【问题讨论】：

请显示您当前的方法，以免我们不小心再次提出相同的建议:)
``diads` 很好，但是您应该能够使用 numpy 方法来创建 diagonals 列表的元素。我也使用coo 进行这种布局。
运行代码时出现错误：ValueError: Diagonal length (index 3: 5 at offset 3) doesn't agree with matrix size (9, 9).
Numerical Recipes pages 1028-1029 有很好的拉普拉斯矩阵的图表和解释。（注意：一些函数有参数nx, ny，一些ny, nx。）

标签： python numpy scipy sparse-matrix

【解决方案1】：

我用数组替换了你对列表的使用：

import numpy as np
from scipy import sparse

nx, ny = 3, 3
N  = nx*ny
main_diag = np.ones(N)*-4.0
side_diag = np.ones(N-1)
side_diag[np.arange(1,N)%4==0] = 0
up_down_diag = np.ones(N-3)
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = sparse.diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

生产

[[-4.  1.  0.  1.  0.  0.  0.  0.  0.]
 [ 1. -4.  1.  0.  1.  0.  0.  0.  0.]
 [ 0.  1. -4.  1.  0.  1.  0.  0.  0.]
 [ 1.  0.  1. -4.  0.  0.  1.  0.  0.]
 [ 0.  1.  0.  0. -4.  1.  0.  1.  0.]
 [ 0.  0.  1.  0.  1. -4.  1.  0.  1.]
 [ 0.  0.  0.  1.  0.  1. -4.  1.  0.]
 [ 0.  0.  0.  0.  1.  0.  1. -4.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0. -4.]]

侧对角线的 [1 1 1 0 1 1 1 0] 模式对吗？

对于像这样的小例子，它可能运行相同的速度，但使用数组而不是列表的大维度应该更快 - 而且它与稀疏底层的 numpy 代码更一致。

像diags 这样的统一对角线看起来很不错。

我只玩过另一个 SO 问题的块格式。 https://stackoverflow.com/a/34124377/901925

coo 适用于由重叠的较小矩阵组成的矩阵，例如有限元刚度。但是将你的对角线重铸成coo 是很乏味的。

对于它的价值，sparse.diags 使用 dia_matrix，已将对角线列表转换为 dia data 矩阵。你可以看一下，但它的布局并不那么明显。要制作csr 矩阵，diags 将此dia 格式转换为coo，然后再转换为csr。但通常你不应该担心所有这些转换。使用对您的问题最有意义的格式，并让sparse 处理转换细节。

如果您想更多地探索块格式，您需要概述如何将您的问题视为块。

【讨论】：

侧对角线的图案是错误的 - 如果是 3x3，它应该是 [1 1 0 1 1 0 1 1]
所以需要使用side_diag[np.arange(1,N)%nx==0] = 0。或者在你的版本中if i%nx == 0:

【解决方案2】：

N 维拉普拉斯算子可以表示为一维拉普拉斯算子的克罗内克积：

import scipy.sparse as sp
def laplacian2D(N):
    diag=np.ones([N*N])
    mat=sp.spdiags([diag,-2*diag,diag],[-1,0,1],N,N)
    I=sp.eye(N)
    return sp.kron(I,mat,format='csr')+sp.kron(mat,I)

Connectivity pattern of 2D Laplacian - img

【讨论】：