3

I have a large numpy array A of shape M*3, whose elements of each row are unique, non-negative integers ranging from 0 to N - 1. In fact, each row corresponds to a triangle in my finite element analysis.

For example, M=4, N=5, and a matrix A looks like the following

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

Now, I need to construct another array B of size M*N, such that

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

The corresponding B for the exemplary A above would be

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

A loop-based code would be

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

But since I have large M and N (of scale 10^6 for each), how can I use good Numpy indexing techniques to accelerate this process? Besides, I feel that sparse matrix techniques are also needed since M * N data of 1 byte is about 10**12, namely 1 000 G.

In general, I feel using numpy's vectorization techniques, such as indexing and broadcasting, look more like an ad-hoc, error-prone activity relying on quite a bit of street smarts (or called art, if you prefer). Are there any programming language efforts that can systematically convert your loop-based code to a high-performance vectorized version?

zell
  • 9,830
  • 10
  • 62
  • 115
  • 1
    Advanced indexing is quite expressive in `numpy` and needs to be learned thoroughly: `A[np.arange(M)[:,None], idx] = 1`, with *A* the result array, *idx* the indices array. [`numba`](https://numba.pydata.org/) can trivially JIT compile your solution. If the result array is filled with >50% non-zweo elements sparse matrices wouldn't be smaller. – Michael Szczesny Aug 18 '22 at 10:21
  • Modern libraries like `pytorch` or `jax` have additional functions for common operations that are simpler to use. In `pytorch` you have [`scatter`](https://pytorch.org/docs/stable/generated/torch.Tensor.scatter_.html): `result.scatter_(dim=1, index=idx, value=1)`, with *result*, *idx* `torch.tensor`. – Michael Szczesny Aug 18 '22 at 10:51
  • I forgot `np.put_along_axis(A, idx, 1, 1)` exists in `numpy`. – Michael Szczesny Aug 18 '22 at 11:07
  • "*A loop-based code would be*" is it ```B``` instead of ```A``` in your loop based code? –  Aug 18 '22 at 12:26
  • 2
    Just to be bit clear what @MichaelSzczesny wrote: ```B = np.zeros((M,N)) ; B[np.arange(M)[:,None],A] = 1``` –  Aug 18 '22 at 13:06
  • @MichaelSzczesny Could you elaborate on "numba can trivially JIT compile your solution. "? I have felt numba frequently gives me compilation errors. – zell Aug 18 '22 at 13:23
  • @MichaelSzczesny Also, do you mind explaining a bit your magic line, as an answer, maybe? – zell Aug 18 '22 at 13:27
  • 2
    Dense matrices are not a good idea here because of the memory space. 10**12 is not 100 G but 1000 G. Because each item takes 8 bytes, this require 8000 GB: far more RAM that what any (realistic) computing machine have so far. Even it it would fit, the matrix will be too big to be read/written quickly (RAM is slow, storage devices are slower). Sparse matrices are a good idea but there are many kind of them designed for specific purposes because there is no 1 kind to rule them all. THus, you should tell us more about the purpose of such a matrix so we can help you further (eg. usage&sparsity) – Jérôme Richard Aug 18 '22 at 20:05

1 Answers1

3

You can directly create a sparse csr-matrix from your data

As you already mentioned in your question a dense matrix consisting of uint8 values would need 1 TB. By using a sparse matrix, this can be reduced to approx. 19 MB as shown in the example bellow.

Creating Inputs with relevant size

This should be included in the question, as it gives a hint on the sparsity of the matrix.

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

Creating a sparse csr-matrix

Have a look at the scipy-doc or for a general overview the wiki article could also be useful.

#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))
max9111
  • 6,272
  • 1
  • 16
  • 33
  • 1
    Good. [Here](https://stackoverflow.com/questions/73415826/how-to-covert-a-large-106-106-numpy-sparse-matrix-to-a-scipy-sparse-matrix/73416136#73416136) is an answer posted simultaneously for a related question of the same OP. Both are completing each other :) . – Jérôme Richard Aug 19 '22 at 11:56