Scipy 稀疏矩阵求幂：a**16 比 a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a 慢？

Question

我正在使用 scipy-0.17 进行简单的稀疏矩阵求幂，a**16。 （注意，不是元素乘法）。然而，在我的机器（运行 Debian stable 和 Ubuntu LTS）上，这比使用 for 循环或做类似a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a 之类的傻事要慢十倍。这没有意义，所以我认为我做错了什么，但是什么？

import scipy.sparse
from time import time

a=scipy.sparse.rand(2049,2049,.002)

print ("Trying exponentiation (a**16)")
t=time()
x=a**16
print (repr(x))
print ("Exponentiation took %f seconds\n" % (time()-t))

print ("Trying expansion (a*a*a*...*a*a)")
t=time()
y=a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a
print (repr(y))
print ("Expansion took %f seconds\n" % (time()-t))

print ("Trying a for loop (z=z*a)")
t=time()
z=scipy.sparse.eye(2049)
for i in range(16):
    z=z*a
print (repr(z))
print ("Looping took %f seconds\n" % (time()-t))

# Sanity check, all approximately the same answer, right? 
assert (abs(x-z)>=1e-9).nnz==0
assert (abs(x-y)>=1e-9).nnz==0

Answer 1

@hpaulj 关于非零数的评论很重要。 当您计算 a 的更高幂时，非零元素的数量 增加。对于稀疏矩阵，计算矩阵的时间 乘积随着非零元素的数量而增加。

用于计算a**16 的算法实际上是：

a2 = a*a
a4 = a2*a2
a8 = a4*a4
a16 = a8*a8

现在看看这些矩阵中非零元素的数量 对于a = sparse.rand(2049, 2049, 0.002)：

matrix      nnz    fraction nnz
  a        8396       0.0020
  a2      34325       0.0082
  a4     521593       0.1240
  a8    4029741       0.9598

在最后一个产品a16 = a8*a8 中，因子为 96% 非零。计算 使用稀疏矩阵乘法的那个产品是慢。 最后一步占用了 97% 的时间来计算 a**16。

另一方面，当您计算 a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a 时， 稀疏矩阵乘法执行 15 次，但只有一次 每个产品中的因子总是有一小部分 (0.002) 非零值，因此每个产品都可以合理地执行 高效。

这表明可能存在计算乘积的最佳策略，平衡乘法次数与因子的稀疏性。比如计算a2 = a*a; a16 = a2*a2*a2*a2*a2*a2*a2*a2比a16 = a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a快：

In [232]: %timeit a2 = a*a; a4 = a2*a2; a8 = a4*a4; a16 = a8*a8
14.4 s ± 199 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [233]: %timeit a16 = a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a
1.77 s ± 4.78 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [234]: %timeit a2 = a*a; a16 = a2*a2*a2*a2*a2*a2*a2*a2
1.42 s ± 3.16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

或者，由于您知道最终结果将是密集的，因此从一开始或在密集矩阵乘法比稀疏矩阵乘法更有效的某个中间步骤切换到标准 numpy 数组。