循环大型稀疏阵列

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是稀疏的，则是.data的multiply属性将与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]