Scipy稀疏矩阵求幂:a ** 16比a * a * a * a * a * a * a * a * a * a * a * a * a * a * a * a * a *
我正在做一个简单的稀疏矩阵求幂,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
这样的傻事慢了十倍。这没有意义,所以我认为我做错了什么,但是什么?Scipy稀疏矩阵求幂:a ** 16比a * 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
@ 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的阵列,然后从开始处或位于稠密矩阵乘法是一些中间步骤比稀疏矩阵乘法更有效。
这是'scipy.sparse'功能请求的坚实基础。 –
无法重现。指数在我的测试中速度更快。 – user2357112
(另外,你正在打印结果的'__repr__'方法,而不是'repr'esentation。) – user2357112
嗯,我想这很适合某人。你在使用scipy-0.17吗? – hackerb9