如何有效地对稀疏矩阵进行lambdify？

Question

我正在使用渐近来构建一个稀疏(N^2 x N^2)矩阵，并尝试将其转换为稀疏scipy矩阵。它被用于求解二维网格上依赖于k向量k= (kx，ky)的Schrödinger方程的有限差分方法。

SciPy的稠密矩阵

我目前的方法是使用lambdify从渐近矩阵A创建函数B，该矩阵为某个k向量提供有限差分矩阵，我将其转换为稀疏矩阵M。

N = 80
kx, ky = sp.symbols('kx ky')

A = dill.load(open("./A_{0}".format(N), "rb"))

B = sp.lambdify([kx, ky], A)

M = sparse.csc_matrix(B(1,0))

dill.dump(B, open("./B_{0}".format(N), "wb"))

我这里的问题是，一旦exceeding达到80，我就会收到>=的内存溢出错误(超过32 GB)，lambdify进程就会被终止。

稀疏Sympy矩阵的模型化

或者，我想直接为lambdify提供一个稀疏矩阵，但我不知道是哪种类型的稀疏矩阵出现……

B = sp.lambdify([kx, ky], sp.SparseMatrix(A))
M = B(1,0)

..。尝试计算M = B(1,0)行中的B会导致错误

File "<lambdifygenerated-1>", line 2, in _lambdifygenerated
TypeError: __init__() got multiple values for argument 'shape'

我使用的是Python 3.7.6、sympy 1.6和scipy 1.4.1

Answer 1

虽然我使用过很多scipy.sparse，但我以前还没有使用过sy.SparseMatrix。查看文档，我可以创建一个简单的文档(在isympy会话中)：

In [4]: A = SparseMatrix(4,4, {(1,1): x, (3,3):y})                              

In [5]: A                                                                       
Out[5]: 
⎡0  0  0  0⎤
⎢          ⎥
⎢0  x  0  0⎥
⎢          ⎥
⎢0  0  0  0⎥
⎢          ⎥
⎣0  0  0  y⎦

In [6]: B=lambdify([x,y], A)                                                    

In [7]: B                                                                       
Out[7]: <function _lambdifygenerated(x, y)>

In [10]: M=B(10,20)                                                             
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-442f4a9e4340> in <module>
----> 1 M=B(10,20)

<lambdifygenerated-1> in _lambdifygenerated(x, y)
      1 def _lambdifygenerated(x, y):
----> 2     return (coo_matrix([x, y], ([1, 3], [1, 3]), shape=(4, 4)))

TypeError: __init__() got multiple values for argument 'shape'

所以我收到了你的错误信息。也清楚地表明了它试图创建什么样的稀疏矩阵。你没看到吗？我很想在这一点上停下来，因为您显然持有有价值的信息！

===

B.__doc__为：

Created with lambdify. Signature:

func(x, y)

Expression:

Matrix([[0, 0, 0, 0], [0, x, 0, 0], [0, 0, 0, 0], [0, 0, 0, y]])

Source code:

def _lambdifygenerated(x, y):
    return (coo_matrix([x, y], ([1, 3], [1, 3]), shape=(4, 4)))

这就是回溯所显示的。很明显，它是在尝试制作一款sparse.coo_matrix。但是有一个错误，缺少一层()。设置矩阵的正确方法是：

In [11]: from scipy import sparse                                               
In [14]: sparse.coo_matrix(([x,y],([1,3], [1,3])), shape=(4,4))                 
Out[14]: 
<4x4 sparse matrix of type '<class 'numpy.object_'>'
    with 2 stored elements in COOrdinate format>

In [15]: print(_)                                                               
  (1, 1)    x
  (3, 3)    y

我不能用symbolic x和y对B做很多事情(甚至不能把它显示为密集的)。

但是如果我把它放在一个函数中，我可以为x和y提供数值。

In [16]: def foo(x,y): 
    ...:     return sparse.coo_matrix(([x,y],([1,3], [1,3])), shape=(4,4)) 
    ...:                                                                        

In [17]: foo(10,20).A                                                           
Out[17]: 
array([[ 0,  0,  0,  0],
       [ 0, 10,  0,  0],
       [ 0,  0,  0,  0],
       [ 0,  0,  0, 20]])

所以SparseMatrix的lambdify是错误的。我不会尝试提出解决办法。