循環大型稀疏陣列

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]