femto 0.0.1.dev0

What's the fastest way to sum a NumPy array?

What’s the fastest way to sum a NumPy array? I don’t know—that’s why I created femto.

femto, written in C, contains several implementations of a sum function. To keep things simple the input array must be at least 2d and the axis of summation cannot be None. Limiting ourselves to the 1d case would have be even simpler. But I am interested in both summing along an axis where the array elements are closely packed in memory (e.g. axis=-1 of a C contiguous array) and where they are widely spaced (axis=0). Both cases require different optimizations.

My goal is to find fast ways to implement reduction functions (sum, mean, std, max, nansum, etc.) that are bound by memory I/O. I chose summation as a test case because very little time is spent with arithmetic which makes it easier to measure improvements from things like manual loop unrolling, SIMD, and OpenMP.

femto, based on code from bottleneck, comes with several benchmark suites:

>>> ss.bench_axis0()
femto performance benchmark
    femto 0.0.1.dev0; Numpy 1.11.2
    Speed is NumPy time divided by femto time
    Score is harmonic mean of speeds

      (1000,1000)(1000,1000)(1000,1000)(1000,1000)
        float64    float32     int64      int32
         axis=0     axis=0     axis=0     axis=0     score
sum00     0.34       0.22       0.63       1.19       0.40
sum01     0.35       0.22       0.65       1.20       0.41
sum02     0.47       0.27       0.69       1.19       0.49
sum03     0.56       0.45       0.87       1.54       0.69
sum04     0.81       0.45       0.56       1.55       0.68
sum10     0.80       0.45       1.20       1.90       0.83
sum11     1.19       1.34       1.20       1.88       1.35
sum12     1.19       0.45       1.21       1.89       0.91
p_sum01   1.30       0.22       2.02       1.36       0.62
p_sum02   1.38       0.23       2.41       1.36       0.63
p_sum03   1.90       1.13       2.73       4.28       1.99
p_sum04   3.22       1.07       2.71       4.33       2.16

>>> ss.bench_axis1()
femto performance benchmark
    femto 0.0.1.dev0; Numpy 1.11.2
    Speed is NumPy time divided by femto time
    Score is harmonic mean of speeds

      (1000,1000)(1000,1000)(1000,1000)(1000,1000)
        float64    float32     int64      int32
         axis=1     axis=1     axis=1     axis=1     score
sum00     0.48       0.44       0.84       1.32       0.63
sum01     0.47       0.45       1.05       1.85       0.68
sum02     1.13       1.29       1.38       3.07       1.47
sum03     1.13       1.27       1.37       3.07       1.47
sum04     1.12       1.29       1.36       3.11       1.47
sum10     1.12       1.28       1.37       3.06       1.47
sum11     1.42       3.04       1.38       3.07       1.92
sum12     1.45       1.29       1.39       3.03       1.59
p_sum01   3.11       3.04       4.16       6.80       3.85
p_sum02   4.37       4.37       5.48      10.65       5.45
p_sum03   4.20       4.35       5.52      10.48       5.37
p_sum04   4.43       4.30       5.49      10.36       5.43

I chose numpy.sum as a benchmark because it is fast and convenient. It should be possible to beat NumPy’s performance. That’s because femto has an unfair advantage. We will not duplicate the pairwise summation NumPy uses to deal with the accumulated round-off error in floating point arrays.

The overall fastest function is the one with the highest benchmark score. Let’s consider the case where we benchmark each function with two arrays (five are used by default in the benchmark) and the speeds are 0.5 (half as fast as NumPy) and 2.0 (twice as fast). What should the overall score be? Some possibilities are the mean (1.25, which is faster than NumPy), the geometric mean (1.0, same as NumPy), or the harmonic mean (0.8, slower). I chose the harmonic mean. If a NumPy program spends equal time summing the two benchmark arrays, each 1 unit of time, then it will take 1/2 + 2 units of time with femto, which is a speed of 2/2.5 = 0.8.

Let’s look at function call overhead by benchmarking with small input arrays:

>>> ss.bench_overhead()
femto performance benchmark
    femto 0.0.1.dev0; Numpy 1.11.2
    Speed is NumPy time divided by femto time
    Score is harmonic mean of speeds

        (10,10)    (10,10)   (100,100)  (100,100)
        float64    float64    float64    float64
         axis=0     axis=1     axis=0     axis=1     score
sum00    16.11      16.11       1.12       0.98       1.96
sum01    13.82      13.22       1.17       1.03       2.03
sum02    15.83      15.91       2.76       2.51       4.51
sum03    16.94      13.38       2.81       2.40       4.42
sum04    17.85      13.45       4.53       2.38       5.19
sum10    16.67      18.16       2.01       2.52       3.96
sum11    18.33      18.62       3.33       3.52       5.78
sum12    18.29      16.14       3.27       3.01       5.30
p_sum01   1.75       1.73       2.59       2.20       2.01
p_sum02   1.80       1.77       2.83       2.58       2.15
p_sum03   1.83       1.87       2.94       2.55       2.21
p_sum04   1.90       1.85       3.30       2.46       2.25

Please help me avoid over optimizing for my particular operating system, CPU, and compiler. Let me know the benchmark results on your system. If you have ideas on how to speed up the code then share them.

License

femto is distributed under the GPL v3+. See the LICENSE file for details.

Requirements

Currently femto only compiles on GNU/Linux. Please help us with getting it to compile on OSX and Windows.

• SSE3, AVX, x86intrin.h, OpenMP

• Python 2.7, 3.4, 3.5

• NumPy 1.11

• gcc

• nose