maximum() / minimum() is slow

[ZeppLu]

julia> using BenchmarkTools; const X = rand(1000, 1000);

julia> maximum(X); @benchmark maximum(X)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.502 ms (0.00% GC)
  median time:      2.591 ms (0.00% GC)
  mean time:        2.764 ms (0.00% GC)
  maximum time:     6.094 ms (0.00% GC)
  --------------
  samples:          1806
  evals/sample:     1

julia> minimum(X); @benchmark minimum(X)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.581 ms (0.00% GC)
  median time:      2.591 ms (0.00% GC)
  mean time:        2.620 ms (0.00% GC)
  maximum time:     3.688 ms (0.00% GC)
  --------------
  samples:          1906
  evals/sample:     1

julia> versioninfo()
Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)

Compared to python:

In [1]: import numpy as np; X = np.random.random([1000, 1000])

In [8]: %timeit np.max(X)
558 µs ± 45.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [9]: %timeit np.min(X)
467 µs ± 4.99 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

[KristofferC]

#26003

[ZeppLu]

But python/numpy handles NaN correctly while not losing speed:

In [5]: X[np.random.randint(1000,size=50), np.random.randint(1000,size=50)] = np.nan

In [6]: np.isnan(np.max(X))
/usr/lib/python3.7/site-packages/numpy/core/fromnumeric.py:83: RuntimeWarning: invalid value encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
Out[6]: True

In [7]: %timeit np.max(X)
/usr/lib/python3.7/site-packages/numpy/core/fromnumeric.py:83: RuntimeWarning: invalid value encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
515 µs ± 54.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

And even though I use mymax2() from #26003 (comment), it seems julia is still slower than numpy:

julia> @benchmark mymax2(X)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.193 ms (0.00% GC)
  median time:      1.294 ms (0.00% GC)
  mean time:        1.467 ms (0.00% GC)
  maximum time:     4.872 ms (0.00% GC)
  --------------
  samples:          3390
  evals/sample:     1

[KristofferC]

Yeah but cf the comment

Since maximum already does its own NaN handling, maybe it could simply pass a faster max to mapreduce?

[ZeppLu]

OK, I just wrote a C function, and found something interesting:

shell> cat maximum.c
double maximum(int length, double* X) {
        double m = X[0];
        for (int i=1; i<length; i++) {
                if (X[i] > m) {
                        m = X[i];
                }
        }
        return m;
}

shell> clang -O2 -shared maximum.c -o libmaximum.so

julia> @benchmark ccall((:maximum, "./libmaximum.so"), Cdouble, (Cint, Ptr{Cdouble}), length(X), Ref(X,1))
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.192 ms (0.00% GC)
  median time:      1.198 ms (0.00% GC)
  mean time:        1.206 ms (0.00% GC)
  maximum time:     1.696 ms (0.00% GC)
  --------------
  samples:          4134
  evals/sample:     1

And a julia implementation plus the ability to check NaN:

julia> function mymax(X)
           m = X[1]
           @inbounds for i = 2:length(X)
               if isnan(X[i])
                   return NaN
               end
               if X[i] > m
                   m = X[i]
               end
           end
           m
       end
mymax (generic function with 1 method)

julia> mymax(X) == maximum(X)
true

julia> isnan(mymax([NaN, 1, Inf]))
true

julia> isnan(mymax([Inf, NaN]))
true

julia> @benchmark mymax(X)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.194 ms (0.00% GC)
  median time:      1.207 ms (0.00% GC)
  mean time:        1.228 ms (0.00% GC)
  maximum time:     2.993 ms (0.00% GC)
  --------------
  samples:          4059
  evals/sample:     1

So basically speaking, a hand written mymax() in julia is almost as fast as C, but current implementation of maximum() does unnecessary NaN-checking, which makes it slower.

Now, my question becomes: Why is np.max() so fast? Can we borrow any idea?

[stevengj]

Now, my question becomes: Why is np.max() so fast? Can we borrow any idea?

np.max uses hand-coded SIMD loops on x86, see here and here, using _mm_min_pd etcetera. Apparently, this instruction forwards NaNs automatically as long as you arrange it so that the possible NaN (i.e. X[i]) is the second argument.

[zbrasseaux]

Hi, can someone point me to a place where I could get involved with this project? I would like to contribute to stuff like this and don't know where to start

[nalimilan]

@zbrasseaux Welcome! I'm not sure this is the easiest issue to get started. Unless you feel like you master the problem (I don't), have a look at issues with the "good first issue" label, or just find an issue which sounds simple enough. You can also hang around on Slack's, #my-first-pr channel.