python - Scipy 稀疏矩阵求幂 : a**16 is slower than a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a?

Question

我正在使用 scipy-0.17 进行简单的稀疏矩阵求幂 a**16。 (注意，不是元素乘法)。但是，在我的机器上(运行 Debian 稳定版和 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 数组。