将稀疏矩阵转换为密集矩阵再转换回稀疏矩阵会降低密度。

Question

我正在使用scipy来生成一个稀疏的有限差分矩阵。我最初是从块矩阵构造它，然后编辑对角线以考虑边界条件。得到的稀疏矩阵属于BSR类型。我发现，如果我将该矩阵转换为稠密矩阵，然后再使用scipy.sparse.BSR_matrix函数将其转换回稀疏矩阵，则会得到比之前更稀疏的矩阵。以下是我用来生成矩阵的代码：

size = (4,4)

xDiff = np.zeros((size[0]+1,size[0]))
ix,jx = np.indices(xDiff.shape)
xDiff[ix==jx] = 1
xDiff[ix==jx+1] = -1

yDiff = np.zeros((size[1]+1,size[1]))
iy,jy = np.indices(yDiff.shape)
yDiff[iy==jy] = 1
yDiff[iy==jy+1] = -1

Ax = sp.sparse.dia_matrix(-np.matmul(np.transpose(xDiff),xDiff))
Ay = sp.sparse.dia_matrix(-np.matmul(np.transpose(yDiff),yDiff))

lap = sp.sparse.kron(sp.sparse.eye(size[1]),Ax) + sp.sparse.kron(Ay,sp.sparse.eye(size[0]))

#set up boundary conditions
BC_diag = np.array([2]+[1]*(size[0]-2)+[2]+([1]+[0]*(size[0]-2)+[1])*(size[1]-2)+[2]+[1]*(size[0]-2)+[2])

lap += sp.sparse.diags(BC_diag)

如果我检查这个矩阵的稀疏性，我会看到以下内容：

lap
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 160 stored elements (blocksize = 4x4) in Block Sparse Row format>

然而，如果我将其转换为密集矩阵，然后再转换回相同的稀疏格式，我会看到一个更稀疏的矩阵：

sp.sparse.bsr_matrix(lap.todense())
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 64 stored elements (blocksize = 1x1) in Block Sparse Row format>

我怀疑发生这种情况的原因是我使用了 sparse.kron 函数构建矩阵，但我的问题是是否有一种方法可以在不先转换为密集矩阵的情况下得到较小的稀疏矩阵，例如，如果我最终想要模拟一个非常大的域。

- algol

---

注意块大小的更改。在第二种情况下，大小为1x1。我想知道lap.tocsr()是否会做同样的事情。我没有太多使用BSR，但我认为它将块存储为密集数组。 - hpaulj

Answer 1

BSR 将数据存储在密集块中：

最初的回答

In [167]: lap.data.shape                                                        
Out[167]: (10, 4, 4)

在这种情况下，这些块有相当多的零。最初的回答。

In [168]: lap1 = lap.tocsr() 
In [170]: lap1                                                                  
Out[170]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 160 stored elements in Compressed Sparse Row format>
In [171]: lap1.data                                                             
Out[171]: 
array([-2.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,
        1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,
        1., -2.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0., -3.,  1.,  0.,
        0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  1., -4.,  1.,  0., 
        ...
        0.,  0.,  1., -2.])

现场清理：

In [172]: lap1.eliminate_zeros()                                                
In [173]: lap1                                                                  
Out[173]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>

如果我在使用kron时指定了csr格式：

最初的回答：

In [181]: lap2 = sparse.kron(np.eye(size[1]),Ax,format='csr') + sparse.kron(Ay,n
     ...: p.eye(size[0]), format='csr')                                         
In [182]: lap2                                                                  
Out[182]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>

- hpaulj

Answer 2

我已经被告知我的答案是不正确的。如果我理解正确，原因是Scipy没有使用Lapack来创建矩阵，而是使用自己的代码来实现这个目的。有趣。尽管意外，但这个信息具有权威性。我将遵从它！

我将保留答案供参考，但不再断言答案是正确的。

一般来说，当涉及到像稀疏矩阵这样的复杂数据结构时，有两种情况：

1. 构造函数提前知道结构的全部内容；或者
2. 该结构被设计成逐步构建，以便只有在结构完成后才知道其全部内容。

复杂数据结构的典型案例是二叉树。你可以在它完成后复制二叉树，使其更有效率。否则，树的标准红黑实现会使一些搜索路径长度最长达到其他路径长度的两倍——通常可以接受，但并不是最优的。

现在，你可能已经知道了所有这些，但我提到它是有原因的。Scipy依赖于Lapack。Lapack带来了几种不同的存储方案。其中两种是

• 一般稀疏和
• 带状

方案。看起来Scipy开始将您的矩阵存储为稀疏矩阵，其中每个非零元素的索引都被明确地存储；但是，在复制时，Scipy注意到带状表示更为适合——因为您的矩阵毕竟是带状的。

- thb

---

啊哈，我被踩了。我的回答有误吗？我想知道。如果回答是错误的，我很乐意删除它，但据我所知，这个答案是正确的。 - thb

scipy.sparse 不使用 Lapack 来创建矩阵。它有自己的代码，是 Python 和 cython 的混合。 - hpaulj