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

Question

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