# Fastest way to iterate over Numpy array

I wrote a function to calculate the gamma coefficient of a clustering. The bottleneck is the comparison of values from dist_withing to dist_between. To speed this up, I tried to adapt and compile it using Cython (I dealt with C only few times). But I don't know, how to rapidly iterate over numpy arrays or if its possible at all to do it faster than

for i in range(len(arr)):
arr[i]


I thought I could use a pointer to the array data and indeed the code runs in only half of the time, but pointer1[i] and pointer2[j] in cdef unsigned int countlower won't give me the expected values from the arrays. So, how to properly and speedy iterate over an array? And where else can be made improvements, even if in this case it would not make such a difference concerning runtime-speed?

# cython: profile=True
import cython

import numpy as np
cimport numpy as np

from scipy.spatial.distance import squareform

DTYPE = np.float
DTYPEint = np.int

ctypedef np.float_t DTYPE_t
ctypedef np.int_t DTYPEint_t

@cython.profile(False)
cdef unsigned int countlower(np.ndarray[DTYPE_t, ndim=1] vec1,
np.ndarray[DTYPE_t, ndim=1] vec2,
int n1, int n2):
# Function output corresponds to np.bincount(v1 < v2)[1]

assert vec1.dtype == DTYPE and vec2.dtype == DTYPE

cdef unsigned int i, j
cdef unsigned int trues = 0
cdef unsigned int* pointer1 = <unsigned int*> vec1.data
cdef unsigned int* pointer2 = <unsigned int*> vec2.data

for i in range(n1):
for j in range(n2):
if pointer1[i] < pointer2[j]:
trues += 1

return trues

def gamma(np.ndarray[DTYPE_t, ndim=2] Y, np.ndarray[DTYPEint_t, ndim=1] part):
assert Y.dtype == DTYPE and part.dtype == DTYPEint

if len(Y) != len(part):
raise ValueError('Distance matrix and partition must have same shape')

# defined locals
cdef unsigned int K, c_label, n_, trues
cdef unsigned int s_plus = 0
cdef unsigned int s_minus = 0

# assigned locals
cdef np.ndarray n_in_ci = np.bincount(part)
cdef int num_clust = len(n_in_ci) - 1
cdef np.ndarray s = np.zeros(len(Y), dtype=DTYPE)

# Partition should have at least two clusters
K = len(set(part))
if K < 2:
return 0
# Loop through clusters
for c_label in range(1, K+1):
dist_within = squareform(Y[part == c_label][:, part == c_label])
dist_between = np.ravel(Y[part == c_label][:, part != c_label])
n1 = len(dist_within)
n2 = len(dist_between)

trues = countlower(dist_within, dist_between, n1, n2)
s_plus += trues
s_minus += n1 * n2 - trues

n_ =  s_plus + s_minus

return (<double>s_plus - <double>s_minus) / <double>n_ if n_ != 0 else 0


Edit1: Passing just the pointers, instead of the arrays to the time-critical function (>99% of time is spent there) made a ~ 10% speed-up. I guess some things just cannot be made faster

@cython.profile(False)
@cython.boundscheck(False)
@cython.nonecheck(False)
cdef unsigned int countlower(double* v1, double* v2, int n1, int n2):
''' Function output corresponds to np.bincount(v1 < v2)[1]'''
''' The upper is not correct. It rather corresponds to
sum([np.bincount(v1[i] < v2)[1] for i in range(len(v1))])'''
cdef unsigned int trues = 0

cdef Py_ssize_t i, j
with nogil, parallel():
for i in prange(n1):
for j in prange(n2):
if v1[i] < v2[j]:
trues += 1
return trues


### 1. Introduction

This question is difficult because:

1. It's not clear what the function countlower does. It's always a good idea to write a docstring for a function, specifying what it does, what arguments it takes, and what it returns. (And test cases are always appreciated.)

2. It's not clear what the role of the arguments n1 and n2 is. The code in the post only ever passes len(v1) for n1 and len(v2) for n2. So is that a requirement? Or is it sometimes possible to pass in other values?

I am going to assume in what follows that:

1. the specification of the countlower function is Return the number of pairs i, j such that v1[i] < v2[j];

2. n1 is always len(v1) and n2 is always len(v2);

3. the Cython details are not essential to the problem, and that it's OK to work in plain Python.

Here's my rewrite of the countlower function. Note the docstring, the doctest, and the simple implementation, which loops over the sequence elements rather than their indices:

def countlower1(v, w):
"""Return the number of pairs i, j such that v[i] < w[j].

>>> countlower1(list(range(0, 200, 2)), list(range(40, 140)))
4500

"""
return sum(x < y for x in v for y in w)


And here's a 1000-element test case, which I'll use in the rest of this answer to compare the performance of various implementations of this function:

>>> v = np.array(list(range(0, 2000, 2)))
>>> w = np.array(list(range(400, 1400)))
>>> from timeit import timeit
>>> timeit(lambda:countlower1(v, w), number=1)
8.449613849865273


### 2. Vectorize

The whole reason for using NumPy is that it enables you to vectorize operations on arrays of fixed-size numeric data types. If you can successfully vectorize an operation, then it executes mostly in C, avoiding the substantial overhead of the Python interpreter.

Whenever you find yourself iterating over the elements of an array, then you're not getting any benefit from NumPy, and this is a sign that it's time to rethink your approach.

So let's vectorize the countlower function. This is easy using a sparse numpy.meshgrid:

import numpy as np

def countlower2(v, w):
"""Return the number of pairs i, j such that v[i] < w[j].

>>> countlower2(np.arange(0, 2000, 2), np.arange(400, 1400))
450000

"""
grid = np.meshgrid(v, w, sparse=True)
return np.sum(grid[0] < grid[1])


Let's see how fast that is on the 1000-element test case:

>>> timeit(lambda:countlower2(v, w), number=1)
0.005706002004444599


That's about 1500 times faster than countlower1.

### 3. Improve the algorithm

The vectorized countlower2 still takes $O(n^2)$ time on arrays of length $O(n)$, because it has to compare every pair of elements. Is it possible to do better than that?

Suppose that I start by sorting the first array v. Then consider an element y from the second array w, and find the point where y would fit into the sorted first array, that is, find i such that v[i - 1] < y <= v[i]. Then y is greater than i elements from v. This position can be found in time $O(\log n)$ using bisect.bisect_left, and so the algorithm as a whole has a runtime of $O(n \log n)$.

Here's a straightforward implementation:

from bisect import bisect_left

def countlower3(v, w):
"""Return the number of pairs i, j such that v[i] < w[j].

>>> countlower3(list(range(0, 2000, 2)), list(range(400, 1400)))
450000

"""
v = sorted(v)
return sum(bisect_left(v, y) for y in w)


This implementation is about three times faster than countlower3 on the 1000-element test case:

>>> timeit(lambda:countlower3(v, w), number=1)
0.0021441911812871695


This shows the importance of finding the best algorithm, not just speeding up the algorithm you've got. Here an $O(n \log n)$ algorithm in plain Python beats a vectorized $O(n^2)$ algorithm in NumPy.

### 4. Vectorize again

Now we can vectorize the improved algorithm, using numpy.searchsorted:

import numpy as np

def countlower4(v, w):
"""Return the number of pairs i, j such that v[i] < w[j].

>>> countlower4(np.arange(0, 20000, 2), np.arange(4000, 14000))
45000000

"""
return np.sum(np.searchsorted(np.sort(v), w))


And this is six times faster still:

>>> timeit(lambda:countlower4(v, w), number=1)
0.0003434771206229925


1. "What does vectorizing mean?" Please read the "What is NumPy?" section of the NumPy documentation, in particular the section starting:

Vectorization describes the absence of any explicit looping, indexing, etc., in the code - these things are taking place, of course, just “behind the scenes” (in optimized, pre-compiled C code).

2. "What is meshgrid?" Please read the documentation for numpy.meshgrid.

I use meshgrid to create a NumPy array grid containing all pairs of elements x, y where x is an element of v and y is an element of w. Then I apply the < function to those pairs, getting an array of Booleans, which I sum. Try it out in the interactive interpreter and see for yourself:

>>> import numpy as np
>>> v = [2, 4, 6]
>>> w = [1, 3, 5]
>>> np.meshgrid(v, w)
[array([[2, 4, 6],
[2, 4, 6],
[2, 4, 6]]), array([[1, 1, 1],
[3, 3, 3],
[5, 5, 5]])]
>>> _[0] < _[1]
array([[False, False, False],
[ True, False, False],
[ True,  True, False]], dtype=bool)
>>> np.sum(_)
3

• fuu, since for me reading code is quite hard, I'll try to do better documentation in future. The assumptions you made are correct. The sorting stuff broadens my way of thinking. What does vectorizing mean in general and in the case of NumPy? What did you do in Section 2 (what is meshgrid useful here)? The main function got about 30 times faster in comparison to when using my edit1 version. No need for using cython anymore. Anyway, if I make a cython file out of that it will take twice as much time. Any idea why is that? – embert Jan 28 '14 at 14:33
• See revised answer. – Gareth Rees Jan 28 '14 at 15:49
• @embert "if I make a cython file out of that it will take twice as much time. Any idea why is that?" The Python→Cython→Python just adds indirection, where the Python type has to be cast into a "special form" for Cython, but that "special form" is never used. Plus, it just ends up in a fast C routine anyway. I've found that inlining these routines into Cython can help by removing overhead, but don't always help much. – Veedrac May 24 '14 at 11:35
• Also look at this if you want to inline it. My guess is not, even if I wasn't months late. – Veedrac May 24 '14 at 11:51

When you deal with performance in cython, I would suggest using the --annotate flag (or use IPython with cython magic that allow you quick iteration with anotate flag too), it will tell you which part of your code may be slow. It generates an Html report with highlighted lines. The more yellow, potentially the slower. You can also click on the line o see the generated C, and generally you just call out into Python world from C when things get slow, like checking array bounds, negative indexing, catching exceptions... So, you might want to use the following decorators on your functions if you know you won't have out of bounds errors, or negative indexing from the end :

@cython.boundscheck(False)
@cython.wraparound(False)


Keep in mind that it you do have out of bounds, you will segfault.

This memview bench might give you ideas.

Yo might also want to look at numpy view, if you like to avoid copy and know things won't be muted (but I think it's the default now)

• The memview bench is quite helpful. Do understand a little bit better now. The -a option also is helpful, however checking for zero division produces some c-code but effects speed insignificantly. If I work with numpy types, should I always use smth like np.float_t instead of double? I need a line profiler, which works with cython! – embert Jan 7 '14 at 16:19
• My knowledge of cython doesn't go that far. Though, it seem to me that division in C is generally much slower than multiplication. Quick search give me @cython.cdivision(True) decorator. looking at the def for n_ I would say that you can simplify logics(x+trues-trues=x) . – Matt Jan 8 '14 at 8:01