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关于nonzeros数量的评论很重要。 当您计算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的阵列，然后从开始处或位于稠密矩阵乘法是一些中间步骤比稀疏矩阵乘法更有效。