what is the fastest way to get the mode of a numpy array

The implementation of scipy.stats.mode has a Python loop for handling the axis argument with multidimensional arrays. The following simple implementation, for one-dimensional arrays only, is faster:

def mode1(x):
    values, counts = np.unique(x, return_counts=True)
    m = counts.argmax()
    return values[m], counts[m]

Here's an example. First, make an array of integers with length 1000000.

In [40]: x = np.random.randint(0, 1000, size=(2, 1000000)).sum(axis=0)

In [41]: x.shape
Out[41]: (1000000,)

Check that scipy.stats.mode and mode1 give the same result.

In [42]: from scipy.stats import mode

In [43]: mode(x)
Out[43]: ModeResult(mode=array([1009]), count=array([1066]))

In [44]: mode1(x)
Out[44]: (1009, 1066)

Now check the performance.

In [45]: %timeit mode(x)
2.91 s ± 18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [46]: %timeit mode1(x)
39.6 ms ± 83.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

2.91 seconds for mode(x) and only 39.6 milliseconds for mode1(x).

Here's one approach based on sorting -

def mode1d(ar_sorted):
    ar_sorted.sort()
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1])
    count = np.empty(idx.size+1,dtype=int)
    count[1:-1] = idx[1:] - idx[:-1]
    count[0] = idx[0] + 1
    count[-1] = ar_sorted.size - idx[-1] - 1
    argmax_idx = count.argmax()

    if argmax_idx==len(idx):
        modeval = ar_sorted[-1]
    else:
        modeval = ar_sorted[idx[argmax_idx]]
    modecount = count[argmax_idx]
    return modeval, modecount

Note that this mutates/changes the input array as it sorts it. So, if you want to keep the input array un-mutated or do mind the input array being sorted, pass a copy.

Sample run on 1M elements -

In [65]: x = np.random.randint(0, 1000, size=(1000000)).astype(float)

In [66]: from scipy.stats import mode

In [67]: mode(x)
Out[67]: ModeResult(mode=array([ 295.]), count=array([1098]))

In [68]: mode1d(x)
Out[68]: (295.0, 1098)

Runtime test

In [75]: x = np.random.randint(0, 1000, size=(1000000)).astype(float)

# Scipy's mode
In [76]: %timeit mode(x)
1 loop, best of 3: 1.64 s per loop

# @Warren Weckesser's soln
In [77]: %timeit mode1(x)
10 loops, best of 3: 52.7 ms per loop

# Proposed in this post
In [78]: %timeit mode1d(x)
100 loops, best of 3: 12.8 ms per loop

With a copy, the timings for mode1d would be comparable to mode1.

I added the two functions mode1 and mode1d from replies above to my script and tried to compare with the scipy.stats.mode.

dir_name="C:/Users/test_mode"
file_name="myfile2.h5"
ds_name="myds"
f_in=os.path.join(dir_name,file_name)

def mode1(x):
    values, counts = np.unique(x, return_counts=True)
    m = counts.argmax()
    return values[m], counts[m]

def mode1d(ar_sorted):
    ar_sorted.sort()
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1])
    count = np.empty(idx.size+1,dtype=int)
    count[1:-1] = idx[1:] - idx[:-1]
    count[0] = idx[0] + 1
    count[-1] = ar_sorted.size - idx[-1] - 1
    argmax_idx = count.argmax()

    if argmax_idx==len(idx):
        modeval = ar_sorted[-1]
    else:
        modeval = ar_sorted[idx[argmax_idx]]
    modecount = count[argmax_idx]
    return modeval, modecount


startTime=time.time()
with h5py.File(f_in, "a") as f1:

        myds=f1[ds_name].value
        time1=time.time()
        file_read_time=time1-startTime
        print(str(file_read_time)+"\t"+"s"+"\t"+str((file_read_time)/60)+"\t"+"min")

        print("mode_scipy=")
        mode_scipy=scipy.stats.mode(myds)
        print(mode_scipy)
        time2=time.time()
        mode_scipy_time=time2-time1
        print(str(mode_scipy_time)+"\t"+"s"+"\t"+str((mode_scipy_time)/60)+"\t"+"min")

        print("mode1=")
        mode1=mode1(myds)
        print(mode1)
        time3=time.time()
        mode1_time=time3-time2
        print(str(mode1_time)+"\t"+"s"+"\t"+str((mode1_time)/60)+"\t"+"min")

        print("mode1d=")
        mode1d=mode1d(myds)
        print(mode1d)
        time4=time.time()
        mode1d_time=time4-time3
        print(str(mode1d_time)+"\t"+"s"+"\t"+str((mode1d_time)/60)+"\t"+"min")

The result from running the script for a numpy array of around 1M is :

mode_scipy= ModeResult(mode=array([ 1.11903353e-06], dtype=float32), count=array([304909])) 938.8368742465973 s
15.647281237443288 min

mode1=(1.1190335e-06, 304909)

0.06500649452209473 s
0.0010834415753682455 min

mode1d=(1.1190335e-06, 304909)

0.06200599670410156 s
0.0010334332784016928 min