python - 为什么添加 scipy.sparse.dia_matrix 实例这么慢？

Question

我正在编写一个数字代码，其中使用scipy.sparse.dia_matrix。我的矩阵非常大(高达约 1000000 x 1000000)，但非常稀疏。有时是三对角线，有时还有更多对角线。

由于各种原因，从编码的角度来看，将几个这样的矩阵(当然大小相同)加在一起是非常方便和清晰的。然而，我发现添加这些稀疏矩阵非常慢。下面的例子说明了我的意思:

import numpy as np
from scipy.sparse import diags, dia_matrix

N = 100000

M1 = diags(diagonals = [np.random.random(N-1), np.random.random(N), np.random.random(N-1)], offsets = [-1, 0, 1])
M2 = diags(diagonals = [np.random.random(N-1), np.random.random(N), np.random.random(N-1)], offsets = [-1, 0, 1])
M3 = diags(diagonals = [np.random.random(N-1), np.random.random(N), np.random.random(N-1)], offsets = [-1, 0, 1])

def simple_add():
    M = M1 + M2 + M3
    
def complicated_add():
    M_ = dia_matrix((N, N))
    for d in [-1, 0, 1]:
        M_.setdiag(M1.diagonal(d) + M2.diagonal(d) + M3.diagonal(d), d)

%timeit simple_add()

%timeit complicated_add()

时序输出为:

16.9 ms ± 730 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
959 µs ± 39.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

我不明白为什么使用 + 运算符将矩阵相加比创建空的二边形矩阵并显式设置对角线慢 17 倍。我可以做些什么来加快速度吗？我更愿意保留带有 + 运算符的更简单的表达式，因为它更具可读性，但不会以计算时间增加一个数量级为代价。

更新:

我提议对 Scipy 进行一项更改，以便更快地添加两个 dia_matrix 实例，经过一番讨论后，我向 Scipy 提交了一个拉取请求，该请求现已合并。因此，将来添加两个 dia_matrix 实例将不再转换为 csr_matrix。

https://github.com/scipy/scipy/pull/14004

Answer 1

diags 从输入列表中生成 dia_matrix:

In [84]: M=sparse.diags([np.arange(1,4),np.arange(1,5),np.arange(1,4)], offsets=[-1,0,1])
In [85]: M
Out[85]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 10 stored elements (3 diagonals) in DIAgonal format>
In [86]: M.offsets
Out[86]: array([-1,  0,  1], dtype=int32)
In [87]: M.data
Out[87]: 
array([[1., 2., 3., 0.],
       [1., 2., 3., 4.],
       [0., 1., 2., 3.]])

对角线列表(不同长度)已转换为带有偏移量的 2 数组。这主要用作输入格式。大多数(如果不是全部)数学都是以 csr 格式实现的。即使在那里，矩阵乘法也是相对的强项。按元素计算明显不如 numpy 数组等价物。

In [89]: Mr=M.tocsr()
In [90]: Mr
Out[90]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 10 stored elements in Compressed Sparse Row format>
In [91]: Mr.data
Out[91]: array([1., 1., 1., 2., 2., 2., 3., 3., 3., 4.])
In [92]: Mr.indices
Out[92]: array([0, 1, 0, 1, 2, 1, 2, 3, 2, 3], dtype=int32)
In [93]: Mr.indptr
Out[93]: array([ 0,  2,  5,  8, 10], dtype=int32)

如果偏移量和形状都相同，dia 格式建议更快的添加速度。

In [94]: M.data += M.data + M.data
In [95]: M.data
Out[95]: 
array([[ 3.,  6.,  9.,  0.],
       [ 3.,  6.,  9., 12.],
       [ 0.,  3.,  6.,  9.]])
In [96]: M.A
Out[96]: 
array([[ 3.,  3.,  0.,  0.],
       [ 3.,  6.,  6.,  0.],
       [ 0.,  6.,  9.,  9.],
       [ 0.,  0.,  9., 12.]])

对于任何稀疏格式，如果所有参数和输出的稀疏度都相同，您通常可以直接在 data 属性上进行数学计算，保持隐含的 0 不变。