寻找良好的索引和稀疏矩阵方法来从现有的矩阵创建一个矩阵

Question

我有一个大的M*3形状的numpy数组，每一行的元素都是唯一的，非负整数，从0到N1。实际上，在我的有限元分析中，每一行都对应一个三角形。

例如，M=4、N=5和矩阵A如下所示

array([[0, 1, 2],
       [0, 2, 3],
       [1, 2, 4],
       [3, 2, 4]])

现在，我需要构造另一个M*N大小的数组B，这样

B[m,n] = 1 if n is in A[m], or else 0 

上述示例性A的相应B将是

1 1 1 0 0
1 0 1 1 0
0 1 1 0 1
0 0 1 1 1

基于循环的代码将是

B = np.zeros((M,N))
for m in range(M):
  for n in B[m]:
    B[m,n]=1

但是由于我有较大的M和N(每个M和N的比例尺为10^6 )，我如何使用良好的Numpy索引技术来加速这个过程？此外，由于1字节的M*N数据约为10**12，即1 000 G，因此也需要稀疏矩阵技术。

总的来说，我觉得使用numpy的矢量化技术，比如索引和广播，看起来更像是一种临时的、容易出错的活动，依赖于相当多的街头智慧(如果你愿意的话，称为艺术)。是否有任何编程语言可以系统地将基于循环的代码转换为高性能的向量化版本？

Answer 1

您可以直接从数据创建稀疏的csr矩阵。

正如您在问题中已经提到的，由uint8值组成的密集矩阵需要1TB。通过使用稀疏矩阵，可以将其简化为大约。19 MB，如下面示例所示。

创建具有相关大小的输入

这应该包括在问题中，因为它暗示了矩阵的稀疏性。

from scipy import sparse
import numpy as np

M=int(1e6)
N=int(1e6)

A=np.random.randint(low=0,high=N,size=(M,3))

创建稀疏csr矩阵

看一看枕部-医生，或者对于一般的概述来说，维基的文章也可能是有用的。

#array of ones, the same size of non-zero values (3 MB if uint8)
data   =np.ones(A.size,dtype=np.uint8)

#you already have the indices, they are expected as an 1D-array (12 MB)
indices=A.reshape(-1)

#every A.shape[1], a new row beginns (4 MB)
indptr =np.arange(0,A.size+1,A.shape[1])

B=sparse.csr_matrix((data,indices,indptr),shape=(M, N))