python – Most efficient way to find mode in numpy array

python – Most efficient way to find mode in numpy array

Check `scipy.stats.mode()` (inspired by @tom10s comment):

``````import numpy as np
from scipy import stats

a = np.array([[1, 3, 4, 2, 2, 7],
[5, 2, 2, 1, 4, 1],
[3, 3, 2, 2, 1, 1]])

m = stats.mode(a)
print(m)
``````

Output:

``````ModeResult(mode=array([[1, 3, 2, 2, 1, 1]]), count=array([[1, 2, 2, 2, 1, 2]]))
``````

As you can see, it returns both the mode as well as the counts. You can select the modes directly via `m[0]`:

``````print(m[0])
``````

Output:

``````[[1 3 2 2 1 1]]
``````

Update

The `scipy.stats.mode` function has been significantly optimized since this post, and would be the recommended method

This is a tricky problem, since there is not much out there to calculate mode along an axis. The solution is straight forward for 1-D arrays, where `numpy.bincount` is handy, along with `numpy.unique` with the `return_counts` arg as `True`. The most common n-dimensional function I see is scipy.stats.mode, although it is prohibitively slow- especially for large arrays with many unique values. As a solution, Ive developed this function, and use it heavily:

``````import numpy

def mode(ndarray, axis=0):
# Check inputs
ndarray = numpy.asarray(ndarray)
ndim = ndarray.ndim
if ndarray.size == 1:
return (ndarray[0], 1)
elif ndarray.size == 0:
raise Exception(Cannot compute mode on empty array)
try:
axis = range(ndarray.ndim)[axis]
except:
raise Exception(Axis {} incompatible with the {}-dimension array.format(axis, ndim))

# If array is 1-D and numpy version is > 1.9 numpy.unique will suffice
if all([ndim == 1,
int(numpy.__version__.split(.)[0]) >= 1,
int(numpy.__version__.split(.)[1]) >= 9]):
modals, counts = numpy.unique(ndarray, return_counts=True)
index = numpy.argmax(counts)
return modals[index], counts[index]

# Sort array
sort = numpy.sort(ndarray, axis=axis)
# Create array to transpose along the axis and get padding shape
transpose = numpy.roll(numpy.arange(ndim)[::-1], axis)
shape = list(sort.shape)
shape[axis] = 1
# Create a boolean array along strides of unique values
strides = numpy.concatenate([numpy.zeros(shape=shape, dtype=bool),
numpy.diff(sort, axis=axis) == 0,
numpy.zeros(shape=shape, dtype=bool)],
axis=axis).transpose(transpose).ravel()
# Count the stride lengths
counts = numpy.cumsum(strides)
counts[~strides] = numpy.concatenate([[0], numpy.diff(counts[~strides])])
counts[strides] = 0
shape = numpy.array(sort.shape)
shape[axis] += 1
shape = shape[transpose]
slices = [slice(None)] * ndim
slices[axis] = slice(1, None)
# Reshape and compute final counts
counts = counts.reshape(shape).transpose(transpose)[slices] + 1

# Find maximum counts and return modals/counts
slices = [slice(None, i) for i in sort.shape]
del slices[axis]
index = numpy.ogrid[slices]
index.insert(axis, numpy.argmax(counts, axis=axis))
return sort[index], counts[index]
``````

Result:

``````In [2]: a = numpy.array([[1, 3, 4, 2, 2, 7],
[5, 2, 2, 1, 4, 1],
[3, 3, 2, 2, 1, 1]])

In [3]: mode(a)
Out[3]: (array([1, 3, 2, 2, 1, 1]), array([1, 2, 2, 2, 1, 2]))
``````

Some benchmarks:

``````In [4]: import scipy.stats

In [5]: a = numpy.random.randint(1,10,(1000,1000))

In [6]: %timeit scipy.stats.mode(a)
10 loops, best of 3: 41.6 ms per loop

In [7]: %timeit mode(a)
10 loops, best of 3: 46.7 ms per loop

In [8]: a = numpy.random.randint(1,500,(1000,1000))

In [9]: %timeit scipy.stats.mode(a)
1 loops, best of 3: 1.01 s per loop

In [10]: %timeit mode(a)
10 loops, best of 3: 80 ms per loop

In [11]: a = numpy.random.random((200,200))

In [12]: %timeit scipy.stats.mode(a)
1 loops, best of 3: 3.26 s per loop

In [13]: %timeit mode(a)
1000 loops, best of 3: 1.75 ms per loop
``````

EDIT: Provided more of a background and modified the approach to be more memory-efficient

python – Most efficient way to find mode in numpy array

If you want to use numpy only:

``````x = [-1, 2, 1, 3, 3]
vals,counts = np.unique(x, return_counts=True)
``````

gives

``````(array([-1,  1,  2,  3]), array([1, 1, 1, 2]))
``````

And extract it:

``````index = np.argmax(counts)
return vals[index]
``````