python - 寻找良好的索引和稀疏矩阵方法以从现有矩阵创建矩阵

Question

我有一个形状为 M*3 的大型 numpy 数组 A，其每一行的元素都是唯一的非负整数，范围从 0 到 N - 1。实际上，在我的有限元分析中，每一行对应一个三角形.

例如，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(每个的比例为 10^6)，我如何使用好的 Numpy 索引技术来加速这个过程？此外，我觉得还需要稀疏矩阵技术，因为1个字节的M * N个数据大约是10**12，即1 000 G。

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

Answer 1

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

正如您在问题中已经提到的，由 uint8 值组成的密集矩阵需要 1 TB。通过使用稀疏矩阵，这可以减少到大约。 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 矩阵

看看 scipy-doc或一般概述 wiki文章也可能有用。

#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))