向量化沒有scipy.sparse的稀疏總和

我試圖做空間衍生物，幾乎設法得到所有的循環出我的代碼，但是當我試圖總結一切在最後我有一個問題。我有一組N~=250k節點。我發現節點對i,j與i.size=j.size=~7.5M是在一定的搜索距離內，最初來自np.triu_indices(n,1)並通過一系列布爾面具來沖洗出不相互影響的節點。現在我想總結其他節點對每個節點的影響。向量化沒有scipy.sparse的稀疏總和

目前，我有這樣的：

def sparseSum(a,i,j,n): 
    return np.array([np.sum(a[np.logical_or(i==k,j==k)],axis=0) for k in range(n)])

這是非常緩慢的。我想要的是矢量化的東西。如果我有scipy我可以做

def sparseSum(a,i,j,n): 
    sp=scipy.sparse.csr_matrix((a,(i,j)),shape=(n,n))+ scipy.sparse.csr_matrix((a,(j,i)),shape=(n,n)) 
    return np.sum(sp, axis=0)

但我在一個不包括scipy的Abaqus實現中這樣做。有沒有辦法做到這種numpy只？

來源

2017-01-23 Daniel F

方法1：這裏有一個方法利用的matrix-multiplication和broadcasting -

K = np.arange(n)[:,None] 
mask = (i == K) | (j == K) 
out = np.dot(mask,a)

方法2：對於少數列的情況下，我們可以使用np.bincount這樣bin-如此 -

def sparseSum(a,i,j,n): 
    if len(a.shape)==1: 
     out=np.bincount(i,a,minlength=n)+np.bincount(j,a) 
    else: 
     ncols = a.shape[1] 
     out = np.empty((n,ncols)) 
     for k in range(ncols): 
      out[:,k] = np.bincount(i,a[:,k],minlength=n) + np.bincount(j,a[:,k]) 
    return out

來源

2017-01-23 10:08:44 Divakar

大概應該補充'n =〜250k'和'i.size =〜7.6M'。我認爲'mask'可能導致'memoryError'，但我會嘗試 –

是的，它掛起。 'mask'是一個238GB的布爾數組，我的計算 –

@DanielForsman我假設'a'是一個二維數組/矩陣。它有多少列？ – Divakar

這裏不是交鑰匙解決方案，而是添加一列稀疏矩陣。它本質上是計算和利用CSC代表

def sparse_col_sums(i, j, a, N): 
    order = np.lexsort(j, i) 
    io, jo, ao = i[order], j[order], a[order] 
    col_bnds = io.searchsorted(np.arange(N)) 
    return np.add.reduceat(ao, col_bnds)

來源

2017-01-23 10:56:52

'i，j'只是對稱影響矩陣的一半。我認爲這個實現需要完整的級別，我不知道如何解決這個問題。 –

@DanielForsman很明顯，你必須用i，j直接調用函數一次，一次翻轉，然後添加返回的向量。我不確定你的全名是什麼意思。代碼當然沒有創建一個密集的表示。事實上，如果scipy.sparse在底層做了類似的事情，我不會感到驚訝。另外，由於您不需要正確的csc rep（不使用'jo'，並且不需要對列進行排序），您可以使用'order = np.argsort（i）'而不是更多昂貴的lexsort。（不可否認，必須排序有點瑕疵。） –

Scipy稀疏使用矩陣乘法來執行行或列總和。 – hpaulj

向量化沒有scipy.sparse的稀疏總和

回答

相關問題