如何避免稀疏矩阵到稠密矩阵的转换

Question

我有一个在Numpy/Scipy中处理稀疏矩阵的简单代码，代码如下所示：

import numpy as np
import scipy as sp
from scipy.sparse import csr_matrix as sparse
from scipy.sparse import vstack
from scipy.linalg import toeplitz

N = 100
d4 = np.array([-22/24, 17/24, 9/24 ,-5/24, 1/24])
s4 = np.array([1/24,-9/8,9/8,-1/24])
n = len(s4)
r4 = sparse((s4, (np.zeros(n), np.arange(n))), shape=[1, N+1])
c4 = sparse(([s4[0]], ([0], [0])), shape=[N-2, 1])
lnd = len(d4)
rd4 = sparse((d4, (np.zeros(lnd), np.arange(lnd))), shape=[1, N+1])
D = sparse(np.concatenate((rd4.todense(), toeplitz(c4.todense(),r4.todense()), np.fliplr(rd4.todense()))))

我想删除sparse到dense的转换，但不知道如何用转换来替换toeplitz函数和fliplr。现在我有这个：

D = vstack([rd4, sparse(toeplitz(c4.todense(),r4.todense())), sparse(np.fliplr(rd4.todense()))])

当然，我可以处理非稀疏矩阵，并在最后转换，但我希望总是处理稀疏矩阵。有更好的主意吗？你知道吗

Answer 1

下面是如何使用scipy.sparse.diags。diags与spdiags类似，但更方便一些。你知道吗

>>> import numpy as np
>>> from scipy import sparse, linalg
>>> 
>>> a = sparse.csr_matrix(np.random.randint(-50, 10, (1, 10)).clip(0, None))
>>> b = sparse.csr_matrix(np.random.randint(-50, 10, (1, 10)).clip(0, None))

这个例子有行向量，对于列向量可以强制转换为csc，然后以同样的方式继续。你知道吗

>>> # dense method for reference
>>> d_toepl = linalg.toeplitz(a.A, b.A)
>>> 
>>> idx = b.indices[0] == 0 # make sure first element of b is ignored
>>> vals, offs = np.r_[a.data, b.data[idx:]], np.r_[-a.indices, b.indices[idx:]]
>>> N = max(a.shape[1], b.shape[1])
>>> dtype = (a[0, ...] + b[0, ...]).dtype
>>> 
>>> s_toepl = sparse.diags(vals, offs, (N, N), dtype=dtype)
>>> 
>>> np.all(d_toepl == s_toepl)
True

fliplr可以通过索引来完成。小问题：不是所有稀疏矩阵类目前都支持索引，您可能必须强制转换。你知道吗

>>> np.all(np.fliplr(d_toepl) == s_toepl.tocsr()[:, ::-1])
True