循环大型稀疏数组

Question

假设我有一个（稀疏）矩阵M 大小(N*N, N*N)。我想从该矩阵中选择元素，其中grid（(n,m) 数组，其中n*m=N）的外积为 True（它是布尔二维数组，na=grid.sum()）。这可以按如下方式完成

result = M[np.outer( grid.flatten(),grid.flatten() )].reshape (( N, N ) )

result 是一个(na,na) 稀疏数组（和na < N）。上一行是我想要实现的：从grid的乘积中获取M中为真的元素，并将不真的元素挤出数组。

随着n 和m（因此N）增长，而M 和result 是稀疏矩阵，我无法在内存或速度方面有效地做到这一点。我尝试过的最接近的是：

result = sp.lil_matrix ( (1, N*N), dtype=np.float32 )
# Calculate outer product
A = np.einsum("i,j", grid.flatten(), grid.flatten())  
cntr = 0
it = np.nditer ( A, flags=['multi_index'] )
while not it.finished:
    if it[0]:
        result[0,cntr] = M[it.multi_index[0], it.multi_index[1]]
        cntr += 1
# reshape result to be a N*N sparse matrix

最后一次整形可以由this approach 完成，但我还没有到那里，因为 while 循环将永远持续下去。

我也尝试过选择 A 的非零元素，然后循环，但这会占用所有内存：

A=np.einsum("i,j", grid.flatten(), grid.flatten())  
nzero = A.nonzero() # This eats lots of memory
cntr = 0
for (i,j) in zip (*nzero):
    temp_mat[0,cntr] = M[i,j]
    cnt += 1

上例中的“n”和“m”约为 300。

Answer 1

我不知道是拼写错误还是代码错误，但您的示例缺少iternext：

R=[]
it = np.nditer ( A, flags=['multi_index'] )
while not it.finished:
    if it[0]:
        R.append(M[it.multi_index])
    it.iternext()

我认为追加到列表比R[ctnr]=... 更简单、更快捷。如果R 是一个常规数组，它是有竞争力的，并且稀疏索引更慢（即使是最快的lil 格式）。

ndindex 将 nditer 的这种用法包装为：

R=[]
for index in np.ndindex(A.shape):
    if A[index]:
        R.append(M[index])

ndenumerate 也可以：

R = []
for index,a in np.ndenumerate(A):
   if a:
       R.append(M[index])

但我想知道您是否真的想推进cntr 的每个it 步骤，而不仅仅是True 案例。否则将result 重塑为(N,N) 没有多大意义。但是那样的话，不就是你的问题吗

M[:N, :N].multiply(A)

或者如果M 是一个密集数组：

M[:N, :N]*A

事实上，如果M 和A 都是稀疏的，那么multiply 的.data 属性将与R 列表相同。

In [76]: N=4
In [77]: M=np.arange(N*N*N*N).reshape(N*N,N*N)
In [80]: a=np.array([0,1,0,1])
In [81]: A=np.einsum('i,j',a,a)
In [82]: A
Out[82]: 
array([[0, 0, 0, 0],
       [0, 1, 0, 1],
       [0, 0, 0, 0],
       [0, 1, 0, 1]])

In [83]: M[:N, :N]*A
Out[83]: 
array([[ 0,  0,  0,  0],
       [ 0, 17,  0, 19],
       [ 0,  0,  0,  0],
       [ 0, 49,  0, 51]])

In [84]: c=sparse.csr_matrix(M)[:N,:N].multiply(sparse.csr_matrix(A))
In [85]: c.data
Out[85]: array([17, 19, 49, 51], dtype=int32)

In [89]: [M[index] for index, a in np.ndenumerate(A) if a]
Out[89]: [17, 19, 49, 51]