Solution

As it turns out, np.sum is a python function calling np.add.reduce. This latter ufunc call is reported by cProfile, I presume because this is still a python object. np.maximum and np.subtract are calls to pure C functions and cProfile thinks these are atomic expressions.

Question

I am trying to optimize a small piece of code that takes more time that I would expect. When running cProfile however, it does not specify which numpy function calls consume the most time, and just lists my function _H as consuming almost all time. Here is the code:

def _H(X, granularity, knots=None, out=None, buf=None, bias=0):
    '''This is the version that I am profiling'''
    np.subtract(knots, granularity*X[..., np.newaxis], out=buf)
    np.maximum(0, buf, out=buf)
    np.sum(buf, axis=1, out=out)
    np.subtract(1-bias, out, out=out)
    np.maximum(-bias, out, out=out)

I am using ufuncs to reduce the number of temporary allocations, which shaved off quite some time. A more pythonic version just for readability:

def _H(X, granularity, knots=None, out=None, buf=None, bias=0):
    '''slow but more readable version'''
    return np.maximum(
         0, 1 - np.maximum(
             0, knots - granularity*X[..., np.newaxis]
         ).sum(1),
     ) - bias

This function is called from a for loop (because the buf array does not fit in memory for large stepsizes):

def H(X, knots, granularity, step=1_000):
    t, d = X.shape
    _, k = knots.shape

    buf = np.empty((step, d, k))
    out = np.empty((t, k))

    for i in range(0, t, step):
        _H(X[i:i + step], granularity, knots=knots,
           out=out[i:i + step], buf=buf[:min(step, t-i)], bias=0)
    return out

The profiling

Here is my profiling snippet:

t = 1_000_000
d = 3
p = 2
X = np.random.uniform(size=(t, d))

cProfile.run('H(X, p, 10_000)', 'profiler')
pstats.Stats('profiler').strip_dirs().sort_stats('tottime').print_stats()

output:

         1178 function calls in 0.689 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      100    0.371    0.004    0.680    0.007 grid.py:333(_H)
      100    0.306    0.003    0.306    0.003 {method 'reduce' of 'numpy.ufunc' objects}
        1    0.008    0.008    0.689    0.689 <string>:1(<module>)
        1    0.001    0.001    0.681    0.681 grid.py:308(H)
      100    0.001    0.000    0.307    0.003 fromnumeric.py:70(_wrapreduction)
      100    0.001    0.000    0.308    0.003 fromnumeric.py:2105(sum)
...

From what I can see, almost half (0.306 seconds) of the time is spent in numpy ufuncs, i.e. np.subtract, np.maximum, and np.sum. However, more than half the time (0.371 seconds) is spent on "other things" in _H. But what exactly? What piece of code is not further specified by cProfile, and why?