具有可变和不可变值混合的Numpy数组

5

我对在忽略数组中的某些值时执行数组操作(点、外部、加等)的最佳/最快方法感兴趣。我主要关注一些值(可能是50%-30%)被忽略且是有效地零,同时使用相当大的数组,可能有100,000到1,000,000个元素。我可以考虑许多解决方案,但似乎没有一个能真正从能够忽略某些值的潜在优势中受益。例如:

import numpy as np
A = np.ones((dim, dim)) # the array to modify
B = np.random.random_integers(0, 1, (dim, dim)) # the values to ignore are 0
C = np.array(B, dtype = np.bool)
D = np.random.random((dim, dim)) # the array which will be used to modify A

# Option 1: zero some values using multiplication.
# some initial tests show this is the fastest
A += B * D

# Option 2: use indexing
# this seems to be the slowest
A[C] += D[C]

# Option 3: use masked arrays
A = np.ma.array(np.ones((dim, dim)), mask = np.array(B - 1, dtype = np.bool))
A += D

编辑1:

正如cyborg所建议的,稀疏数组可能是另一种选择。不幸的是,我对该包不太熟悉,无法获得可能会有的速度优势。例如,如果我有一个由稀疏矩阵A定义的限制连通性的加权图,另一个定义连通性(1 = 连通,0 = 不连通)的稀疏矩阵B和一个密集的numpy矩阵C,我希望能够做类似于A = A + B.multiply(C)的操作,并利用AB的稀疏性。


在您的应用程序中,B是否保持不变,还是在处理过程中可能会经常改变?如果它不会改变,您可以采用类似于掩码数组到压缩数组的方法。我假设从索引跳转是需要很长时间的部分,这是根据选项2的结果判断的。 - yosukesabai
另外,你在掩码数组中发现了什么?我很好奇。 - yosukesabai
根据我的目的,B将保持不变。压缩数组看起来很合适,但似乎没有任何后续操作,并且由于它在幕后使用掩码数组进行所有算术运算,我怀疑它是否会更好,但谁知道呢。在我进行的简单测试中,掩码数组比仅使用乘法将值清零慢了约4倍。 - alto
A += B.multiply(C) 是正确的做法。你为什么说你没有优势? - cyborg
如果 B 是稀疏的并且 C 是稠密的,B.multiply(C) 将会是一个稠密矩阵。我想使用稀疏矩阵的任何好处都将被转换结果回到稀疏矩阵的成本所抵消,然后再加到 A 中。也许有更好的方法,但是我不太熟悉 scipy.sparse。另外我的语法是错误的。在scipy中,+=未定义为稀疏矩阵,应该是 A = A + B.multiply(C) - alto
请查看我的编辑,关于对 scipy.sparse.coo_matrix 扩展的说明。 - cyborg
1个回答

1

如果密度低于10%,使用稀疏矩阵可以获得性能提升。使用稀疏矩阵可能会更快,这取决于是否包括构建矩阵所需的时间。

import timeit

setup=\
'''
import numpy as np
dim=1000
A = np.ones((dim, dim)) # the array to modify
B = np.random.random_integers(0, 1, (dim, dim)) # the values to ignore are 0
C = np.array(B, dtype = np.bool)
D = np.random.random((dim, dim)) # the array which will be used to modify A
'''

print('mult    '+str(timeit.timeit('A += B * D', setup, number=3)))

print('index   '+str(timeit.timeit('A[C] += D[C]', setup, number=3)))

setup3 = setup+\
''' 
A = np.ma.array(np.ones((dim, dim)), mask = np.array(B - 1, dtype = np.bool))
'''
print('ma      ' + str(timeit.timeit('A += D', setup3, number=3)))

setup4 = setup+\
''' 
from scipy import sparse
S = sparse.csr_matrix(C)
DS = S.multiply(D)
'''
print('sparse- '+str(timeit.timeit('A += DS', setup4, number=3)))

setup5 = setup+\
''' 
from scipy import sparse
'''
print('sparse+ '+str(timeit.timeit('S = sparse.csr_matrix(C); DS = S.multiply(D); A += DS', setup4, number=3)))

setup6 = setup+\
'''
from scipy import sparse
class Sparsemat(sparse.coo_matrix):
    def __iadd__(self, other):
        self.data += other.data
        return self
A = Sparsemat(sparse.rand(dim, dim, 0.5, 'coo')) # the array to modify
D = np.random.random((dim, dim)) # the array which will be used to modify A
anz = A.nonzero()
'''
stmt6=\
'''
DS = Sparsemat((D[anz[0],anz[1]], anz), shape=A.shape) # new graph based on random weights
A += DS
'''
print('sparse2 '+str(timeit.timeit(stmt6, setup6, number=3)))

输出:

mult    0.0248420299535
index   0.32025789431
ma      0.1067024434
sparse- 0.00996273276303
sparse+ 0.228869672266
sparse2 0.105496183846

编辑:您可以使用上面的代码(setup6)来扩展scipy.sparse.coo_matrix。它保持了稀疏格式。


似乎稀疏数组可能是一个合理的选择。不幸的是,我对这个包不够熟悉,无法利用它,即我总是最终得到密集矩阵。我已经编辑了我的原始问题,提供了一个更具体的例子。 - alto
1
一个矩阵如果有50%到30%的零元素,就不是稀疏矩阵。 - Simon Bergot

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接