向量化没有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的稀疏总和

回答

相关问题