Fast way to generate large-scale random ndarray

Question

I want to generate a random matrix of shape (1e7, 800). But I find numpy.random.rand() becomes very slow at this scale. Is there a quicker way?

Answer 1

A simple way to do that is to write a multi-threaded implementation using Numba:

import numba as nb
import random

@nb.njit('float64[:,:](int_, int_)', parallel=True)
def genRandom(n, m):
    res = np.empty((n, m))

    # Parallel loop
    for i in nb.prange(n):
        for j in range(m):
            res[i, j] = np.random.rand()

    return res

This is 6.4 times faster than np.random.rand() on my 6-core machine.

Note that using 32-bit floats may help to speed up a bit the computation although the precision will be lower.

Answer 2

Numba is a good option, another option that might work well is dask.array, which will create lazy blocks of numpy arrays and perform parallel computations on blocks. On my machine I get a factor of 2 improvement in speed (for 1e6 x 1e3 matrix since I don't have enough memory on my machine).

rows = 10**6
cols = 10**3
import dask.array as da
x = da.random.random(size=(rows, cols)).compute() # takes about 5 seconds

# import numpy as np
# x = np.random.rand(rows, cols) # takes about 10 seconds

Note that .compute at the end is only to bring the computed array into memory, however in general you can continue to exploit the parallel computations with dask to get much faster calculations (that can also scale beyond a single machine), see docs.

Answer 3

An attempt to find an answer from answers given till now:

I just wrote a script which is compiled from already given (by SultanOrazbayev and Jérôme Richard) answers and contains 3 functions for each numba, dask and numpy approach and measure the time spent for n number of different sized arrays.

The code

import dask.array as da
import matplotlib.pyplot as plt
import numba as nb
import timeit
import numpy as np


@nb.njit('float64[:,:](int_, int_)', parallel=True)
def nmb(n, m):
    res = np.empty((n, m))

    # Parallel loop
    for i in nb.prange(n):
        for j in range(m):
            res[i, j] = np.random.rand()

    return res


def nmp(n, m):
    return np.random.random((n, m))


def dask(n, m):
    return da.random.random(size=(n, m)).compute()


if __name__ == '__main__':
    data = []
    for i in range(1, 16):
        dmm = 2 ** i
        s_nmb = timeit.default_timer()
        nmb(dmm, dmm)
        e_nmb = timeit.default_timer()

        s_nmp = timeit.default_timer()
        nmp(dmm, dmm)
        e_nmp = timeit.default_timer()

        s_dask = timeit.default_timer()
        dask(dmm, dmm)
        e_dask = timeit.default_timer()

        data.append([
            dmm,
            e_nmb - s_nmb,
            e_nmp - s_nmp,
            e_dask - s_dask
        ])

    data = np.array(data)

    plt.plot(data[:, 0], data[:, 1], "-r", label="Numba")
    plt.plot(data[:, 0], data[:, 2], "-g", label="Numpy")
    plt.plot(data[:, 0], data[:, 3], "-b", label="Dask")

    plt.xlabel("Number of Element on each axes")
    plt.ylabel("Time spent (s)")
    plt.legend()
    plt.show()

The result