优化与numpy的sparsey阵列

我有一个缓慢的numpy的操作挣扎，使用Python 3优化与numpy的sparsey阵列

的操作，我有以下操作：

np.sum(np.log(X.T * b + a).T, 1)

其中

(30000,1000) = X.shape 
(1000,1) = b.shape 
(1000,1) = a.shape

我的问题这个操作很慢（大约1.5秒），并且它在一个循环内，所以它重复了大约100次，这使得我的代码的运行时间非常长。

我想知道是否有更快的实现这个功能。

也许有用的事实：X非常稀疏（只有0.08％的条目非零），但是一个NumPy数组。

来源

2017-10-14 Ulderique Demoitre

什么你用什么结构来存储X？你有没有看过scipy.sparse.csc_matrix？ –

这是一个现在的数组。不，我没有看它.. –

冒险将这部分信息添加到帖子中，这是至关重要的。此外，人们可能会混淆NumPy稀疏矩阵和scipy稀疏矩阵，所以也编辑过标题。这个“稀疏”是我可以想象的最好的，而不会与SciPy稀疏矩阵混淆。 – Divakar

我们可以优化这似乎是瓶颈对数运算和被超越的一个函数可以用numexpr module加速，然后sum-reduce用NumPy加速，因为NumPy做得好得多，因此给了我们一个混合版本，就像这样 -

import numexpr as ne 

def numexpr_app(X, a, b): 
    XT = X.T 
    return ne.evaluate('log(XT * b + a)').sum(0)

仔细观察广播业务：XT * b + a，我们看到广播有两个阶段，我们可以进一步优化。目的是要看看是否可以减少到一个阶段，这在某些部门似乎是可能的。这给我们一个稍微修改后的版本，如下图所示 -

def numexpr_app2(X, a, b): 
    ab = (a/b) 
    XT = X.T 
    return np.log(b).sum() + ne.evaluate('log(ab + XT)').sum(0)

运行测试和验证

原始的方法 -

def numpy_app(X, a, b): 
    return np.sum(np.log(X.T * b + a).T, 1)

计时 -

In [111]: # Setup inputs 
    ...: density = 0.08/100 # 0.08 % sparse 
    ...: m,n = 30000, 1000 
    ...: X = scipy.sparse.rand(m,n,density=density,format="csr").toarray() 
    ...: a = np.random.rand(n,1) 
    ...: b = np.random.rand(n,1) 
    ...: 

In [112]: out0 = numpy_app(X, a, b) 
    ...: out1 = numexpr_app(X, a, b) 
    ...: out2 = numexpr_app2(X, a, b) 
    ...: print np.allclose(out0, out1) 
    ...: print np.allclose(out0, out2) 
    ...: 
True 
True 

In [114]: %timeit numpy_app(X, a, b) 
1 loop, best of 3: 691 ms per loop 

In [115]: %timeit numexpr_app(X, a, b) 
10 loops, best of 3: 153 ms per loop 

In [116]: %timeit numexpr_app2(X, a, b) 
10 loops, best of 3: 149 ms per loop

只是为了证明在启动log部分与原NumPy的方法的瓶颈规定的观察，这里是它的时机 - 在其改善是显著

In [44]: %timeit np.log(X.T * b + a) 
1 loop, best of 3: 682 ms per loop

In [120]: XT = X.T 

In [121]: %timeit ne.evaluate('log(XT * b + a)') 
10 loops, best of 3: 142 ms per loop

来源

2017-10-14 15:55:18 Divakar

如果b有零条目，该怎么办？ –

@ P.Camilleri你的意思是'b'有零？是的，在那种情况下，'numexpr_app'就是要走的路。 – Divakar

@Divakar这非常有趣，谢谢。我只能使用numexpr_app，因为b有负值。但是，我可能会得到与输出略有不同的值，比较numpy_app和numexpr_app？ –

有点不清楚你为什么要做np.sum(your_array.T, axis=1)而不是np.sum(your_array, axis=0)。

您可以使用scipy sparse matrix：（对于X使用列压缩格式，使XT的行压缩，因为你用b具有XT的一个排的形状乘）

X_sparse = scipy.sparse.csc_matrx(X)

并更换XT * b通过：

X_sparse.T.multiply(b)

但是，如果一个不稀疏，它不会帮助你尽可能多。

这些都是速度提升我获得这个操作：

In [16]: %timeit X_sparse.T.multiply(b) 
The slowest run took 10.80 times longer than the fastest. This could mean that an intermediate result is being cached. 
1000 loops, best of 3: 374 µs per loop 

In [17]: %timeit X.T * b 
10 loops, best of 3: 44.5 ms per loop

有：

import numpy as np 
from scipy import sparse 

X = np.random.randn(30000, 1000) 
a = np.random.randn(1000, 1) 
b = np.random.randn(1000, 1) 
X[X < 3] = 0 
print(np.sum(X != 0)) 
X_sparse = sparse.csc_matrix(X)

来源

2017-10-14 14:42:47

从稠密矩阵创建稀疏矩阵需要时间。而且'a'的加入使得结果变得密集，并且可能很昂贵。 – hpaulj

优化与numpy的sparsey阵列

回答

相关问题