Measures branching in C++ from python

Links: notebook, html, PDF, python, slides, slides(2), GitHub

This notebooks looks into a couple of ways to write code, which one is efficient, which one is not when it comes to write fast and short loops. Both experiments are around branching. The notebook relies on C++ code implemented in cbenchmark.cpp and repeat_fct.h.

from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline

numpy is multithreaded. For an accurate comparison, this needs to be disabled. This can be done as follows or by setting environment variable MKL_NUM_THREADS=1.

try:
    import mkl
    mkl.set_num_threads(1)
except ModuleNotFoundError as e:
    print('mkl not found', e)
    import os
    os.environ['MKL_NUM_THREADS']='1'

First experiment: comparison C++ syntax

This all started with article Why is it faster to process a sorted array than an unsorted array?. It compares different implementation fo the following function for which we try different implementations for the third line in next cell. The last option is taken Checking whether a number is positive or negative using bitwise operators which avoids branching.

# int nb = 0;
# for(auto it = values.begin(); it != values.end(); ++it)
#     if (*it >= th) nb++; // this line changes
#     if (*it >= th) nb++; // and is repeated 10 times inside the loop.
#     // ... 10 times
# return nb;

The third line is also repeated 10 times to avoid the loop being too significant.

from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_A, measure_scenario_B
from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_C, measure_scenario_D
from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_E, measure_scenario_F
from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_G, measure_scenario_H
from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_I, measure_scenario_J
import pandas

def test_benchmark(label, values, th, repeat=10, number=20):
    funcs = [(k, v) for k, v in globals().copy().items() if k.startswith("measure_scenario")]
    rows = []
    for k, v in funcs:
        exe = v(values, th, repeat, number)
        d = exe.todict()
        d['doc'] = v.__doc__.split('``')[1]
        d['label'] = label
        d['name'] = k
        rows.append(d)
    df = pandas.DataFrame(rows)
    return df

test_benchmark("sorted", list(range(10)), 5)
average deviation doc label max_exec min_exec name number repeat
0 6.123500e-07 6.180500e-07 if (values[i] >= th) ++nb; sorted 1.580000e-06 0.000001 measure_scenario_A 20.0 10.0
1 6.123500e-07 6.180500e-07 if (*it >= th) ++nb; sorted 1.580000e-06 0.000001 measure_scenario_B 20.0 10.0
2 5.926000e-07 5.926001e-07 if (*it >= th) nb++; sorted 1.186000e-06 0.000001 measure_scenario_C 20.0 10.0
3 7.905000e-08 1.581001e-07 nb += *it >= th ? 1 : 0; sorted 3.960000e-07 0.000000 measure_scenario_D 20.0 10.0
4 5.926000e-07 5.926001e-07 if (*it >= th) nb += 1; sorted 1.186000e-06 0.000001 measure_scenario_E 20.0 10.0
5 1.383000e-07 1.884712e-07 nb += (*it - th) >= 0 ? 1 : 0; sorted 3.960000e-07 0.000000 measure_scenario_F 20.0 10.0
6 1.185500e-07 1.810882e-07 nb += (*it - th) < 0 ? 1 : 0; sorted 3.960000e-07 0.000000 measure_scenario_G 20.0 10.0
7 9.875000e-08 1.710400e-07 nb += *it < th ? 1 : 0; sorted 3.950000e-07 0.000000 measure_scenario_H 20.0 10.0
8 9.875000e-08 1.710400e-07 nb += 1 ^ ((unsigned int)(*it) >> (sizeof(int)... sorted 3.950000e-07 0.000000 measure_scenario_I 20.0 10.0
9 9.875000e-08 1.710400e-07 nb += values[i] >= th ? 1 : 0; sorted 3.950000e-07 0.000000 measure_scenario_J 20.0 10.0

Times are not very conclusive on such small lists.

values = list(range(100000))
df_sorted = test_benchmark("sorted", values, len(values)//2, repeat=200)
df_sorted
average deviation doc label max_exec min_exec name number repeat
0 0.248125 0.0 if (values[i] >= th) ++nb; sorted 0.052481 0.015181 measure_scenario_A 20.0 200.0
1 0.258902 0.0 if (*it >= th) ++nb; sorted 0.045524 0.013611 measure_scenario_B 20.0 200.0
2 0.231941 0.0 if (*it >= th) nb++; sorted 0.044393 0.013159 measure_scenario_C 20.0 200.0
3 0.011715 0.0 nb += *it >= th ? 1 : 0; sorted 0.006699 0.000985 measure_scenario_D 20.0 200.0
4 0.169927 0.0 if (*it >= th) nb += 1; sorted 0.032840 0.013392 measure_scenario_E 20.0 200.0
5 0.019457 0.0 nb += (*it - th) >= 0 ? 1 : 0; sorted 0.003985 0.001550 measure_scenario_F 20.0 200.0
6 0.018889 0.0 nb += (*it - th) < 0 ? 1 : 0; sorted 0.003149 0.001426 measure_scenario_G 20.0 200.0
7 0.010478 0.0 nb += *it < th ? 1 : 0; sorted 0.001324 0.001007 measure_scenario_H 20.0 200.0
8 0.017460 0.0 nb += 1 ^ ((unsigned int)(*it) >> (sizeof(int)... sorted 0.003487 0.001560 measure_scenario_I 20.0 200.0
9 0.014523 0.0 nb += values[i] >= th ? 1 : 0; sorted 0.002035 0.001389 measure_scenario_J 20.0 200.0

The article some implementations will be slower if the values are not sorted.

import random
random.shuffle(values)
values = values.copy()
values[:10]
[40500, 74171, 82572, 33524, 98173, 79727, 77980, 68334, 21586, 56113]
df_shuffled = test_benchmark("shuffled", values, len(values)//2, repeat=200)
df_shuffled
average deviation doc label max_exec min_exec name number repeat
0 0.159574 0.0 if (values[i] >= th) ++nb; shuffled 0.027933 0.013216 measure_scenario_A 20.0 200.0
1 0.237726 0.0 if (*it >= th) ++nb; shuffled 0.035669 0.017613 measure_scenario_B 20.0 200.0
2 0.242710 0.0 if (*it >= th) nb++; shuffled 0.035139 0.018303 measure_scenario_C 20.0 200.0
3 0.016829 0.0 nb += *it >= th ? 1 : 0; shuffled 0.002218 0.001059 measure_scenario_D 20.0 200.0
4 0.229266 0.0 if (*it >= th) nb += 1; shuffled 0.035332 0.017033 measure_scenario_E 20.0 200.0
5 0.027358 0.0 nb += (*it - th) >= 0 ? 1 : 0; shuffled 0.004077 0.001695 measure_scenario_F 20.0 200.0
6 0.025477 0.0 nb += (*it - th) < 0 ? 1 : 0; shuffled 0.003006 0.001531 measure_scenario_G 20.0 200.0
7 0.017814 0.0 nb += *it < th ? 1 : 0; shuffled 0.002583 0.001009 measure_scenario_H 20.0 200.0
8 0.030216 0.0 nb += 1 ^ ((unsigned int)(*it) >> (sizeof(int)... shuffled 0.005585 0.001666 measure_scenario_I 20.0 200.0
9 0.026263 0.0 nb += values[i] >= th ? 1 : 0; shuffled 0.003680 0.001476 measure_scenario_J 20.0 200.0
df = pandas.concat([df_sorted, df_shuffled])
dfg = df[["doc", "label", "average"]].pivot("doc", "label", "average")
ax = dfg.plot.bar(rot=30)
labels = [l.get_text() for l in ax.get_xticklabels()]
ax.set_xticklabels(labels, ha='right')
ax.set_title("Comparison of all implementations");
../_images/cbenchmark_branching_16_0.png

It seems that inline tests (cond ? value1 : value2) do not stop the branching and it should be used whenever possible.

sdf = df[["doc", "label", "average"]]
dfg2 = sdf[sdf.doc.str.contains('[?^]')].pivot("doc", "label", "average")
ax = dfg2.plot.bar(rot=30)
labels = [l.get_text() for l in ax.get_xticklabels()]
ax.set_xticklabels(labels, ha='right')
ax.set_title("Comparison of implementations using ? :");
../_images/cbenchmark_branching_18_0.png
sdf = df[["doc", "label", "average"]]
dfg2 = sdf[sdf.doc.str.contains('if')].pivot("doc", "label", "average")
ax = dfg2.plot.bar(rot=30)
labels = [l.get_text() for l in ax.get_xticklabels()]
ax.set_xticklabels(labels, ha='right')
ax.set_ylim([0.10, 0.20])
ax.set_title("Comparison of implementations using tests");
../_images/cbenchmark_branching_19_0.png

sorted, not sorted does not seem to have a real impact in this case. It shows branching really slows down the execution of a program. Branching happens whenever the program meets a loop condition or a test. Iterator *it are faster than accessing an array with notation [i] which adds a cost due to an extra addition.

Second experiment: dot product

The goal is to compare the dot product from numpy.dot and a couple of implementation in C++ which look like this:

# float vector_dot_product_pointer(const float *p1, const float *p2, size_t size)
# {
#     float sum = 0;
#     const float * end1 = p1 + size;
#     for(; p1 != end1; ++p1, ++p2)
#         sum += *p1 * *p2;
#     return sum;
# }
#
#
# float vector_dot_product(py::array_t<float> v1, py::array_t<float> v2)
# {
#     if (v1.ndim() != v2.ndim())
#         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
#     if (v1.ndim() != 1)
#         throw std::runtime_error("Vector v1 and v2 must be vectors.");
#     return vector_dot_product_pointer(v1.data(0), v2.data(0), v1.shape(0));
# }

numpy vs C++

%matplotlib inline
import numpy

def simple_dot(values):
    return numpy.dot(values, values)

values = list(range(10000000))
values = numpy.array(values, dtype=numpy.float32)
vect = values / numpy.max(values)
simple_dot(vect)
3333333.2
vect.dtype
dtype('float32')
from timeit import Timer

def measure_time(stmt, context, repeat=10, number=50):
    tim = Timer(stmt, globals=context)
    res = numpy.array(tim.repeat(repeat=repeat, number=number))
    mean = numpy.mean(res)
    dev = numpy.mean(res ** 2)
    dev = (dev - mean**2) ** 0.5
    return dict(average=mean, deviation=dev, min_exec=numpy.min(res),
                max_exec=numpy.max(res), repeat=repeat, number=number,
                size=context['values'].shape[0])

measure_time("simple_dot(values)", context=dict(simple_dot=simple_dot, values=vect))
{'average': 0.16509191320000055,
 'deviation': 0.0043107628936845475,
 'min_exec': 0.16100208099999946,
 'max_exec': 0.173292435999997,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("simple_dot(values)", repeat=100,
                     context=dict(simple_dot=simple_dot, values=vect[:i].copy()))
    res.append(t)

import pandas
dot = pandas.DataFrame(res)
dot.tail()
average deviation max_exec min_exec number repeat size
75 0.001635 0.000348 0.002779 0.000963 50 100 187510
76 0.001632 0.000309 0.002114 0.001150 50 100 190010
77 0.001635 0.000326 0.001984 0.000977 50 100 192510
78 0.001627 0.000322 0.002032 0.001035 50 100 195010
79 0.001720 0.000325 0.002444 0.001188 50 100 197510
res = []
for i in range(100000, 10000000, 1000000):
    t = measure_time("simple_dot(values)", repeat=10,
                     context=dict(simple_dot=simple_dot, values=vect[:i].copy()))
    res.append(t)

huge_dot = pandas.DataFrame(res)
huge_dot.head()
average deviation max_exec min_exec number repeat size
0 0.001156 0.000448 0.002394 0.000825 50 10 100000
1 0.014210 0.001607 0.016698 0.010715 50 10 1100000
2 0.038344 0.009988 0.057081 0.029680 50 10 2100000
3 0.047936 0.002603 0.053441 0.044621 50 10 3100000
4 0.063141 0.001562 0.065169 0.060751 50 10 4100000
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0])
huge_dot.plot(x='size', y="average", ax=ax[1])
ax[0].set_title("numpy dot product execution time");
ax[1].set_title("numpy dot product execution time");
../_images/cbenchmark_branching_30_0.png

Now the custom implementation. We start with an empty function to get a sense of the cost due to to pybind11.

from cpyquickhelper.numbers.cbenchmark_dot import empty_vector_dot_product
empty_vector_dot_product(vect, vect)
0.0
def empty_c11_dot(vect):
    return empty_vector_dot_product(vect, vect)

measure_time("empty_c11_dot(values)", context=dict(empty_c11_dot=empty_c11_dot, values=vect), repeat=10)
{'average': 0.0001688884999992979,
 'deviation': 6.0375089806571085e-05,
 'min_exec': 6.716000000039912e-05,
 'max_exec': 0.0002374319999915997,
 'repeat': 10,
 'number': 50,
 'size': 10000000}

Very small. It should not pollute our experiments.

from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product
vector_dot_product(vect, vect)
3334629.0
def c11_dot(vect):
    return vector_dot_product(vect, vect)

measure_time("c11_dot(values)", context=dict(c11_dot=c11_dot, values=vect), repeat=10)
{'average': 1.2002336787000005,
 'deviation': 0.03429446321745382,
 'min_exec': 1.1659099270000013,
 'max_exec': 1.2893128430000047,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("c11_dot(values)", repeat=10,
                     context=dict(c11_dot=c11_dot, values=vect[:i].copy()))
    res.append(t)

import pandas
cus_dot = pandas.DataFrame(res)
cus_dot.tail()
average deviation max_exec min_exec number repeat size
75 0.020400 0.001467 0.023072 0.018215 50 10 187510
76 0.021511 0.001281 0.023578 0.019078 50 10 190010
77 0.022221 0.001384 0.024188 0.019615 50 10 192510
78 0.023844 0.001819 0.026112 0.021054 50 10 195010
79 0.022494 0.001679 0.024817 0.020021 50 10 197510
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0], label="numpy")
cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
ax[0].set_title("numpy and custom dot product execution time");
ax[1].set_title("numpy and custom dot product execution time");
../_images/cbenchmark_branching_38_0.png

Pretty slow. Let’s see what it does to compute dot product 16 by 16.

Use of branching: 16 multplications in one row

The code looks like what follows. If there is more than 16 multiplications left, we use function vector_dot_product_pointer16, otherwise, there are done one by one like the previous function.

# float vector_dot_product_pointer16(const float *p1, const float *p2)
# {
#     float sum = 0;
#
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#     sum += *(p1++) * *(p2++);
#
#     return sum;
# }
#
# #define BYN 16
#
# float vector_dot_product_pointer16(const float *p1, const float *p2, size_t size)
# {
#     float sum = 0;
#     size_t i = 0;
#     if (size >= BYN) {
#         size_t size_ = size - BYN;
#         for(; i < size_; i += BYN, p1 += BYN, p2 += BYN)
#             sum += vector_dot_product_pointer16(p1, p2);
#     }
#     for(; i < size; ++p1, ++p2, ++i)
#         sum += *p1 * *p2;
#     return sum;
# }
#
# float vector_dot_product16(py::array_t<float> v1, py::array_t<float> v2)
# {
#     if (v1.ndim() != v2.ndim())
#         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
#     if (v1.ndim() != 1)
#         throw std::runtime_error("Vector v1 and v2 must be vectors.");
#     return vector_dot_product_pointer16(v1.data(0), v2.data(0), v1.shape(0));
# }
from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16
vector_dot_product16(vect, vect)
3333331.75
def c11_dot16(vect):
    return vector_dot_product16(vect, vect)

measure_time("c11_dot16(values)", context=dict(c11_dot16=c11_dot16, values=vect), repeat=10)
{'average': 0.545918839700002,
 'deviation': 0.02039451000541943,
 'min_exec': 0.5292382869999983,
 'max_exec': 0.6024312859999981,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("c11_dot16(values)", repeat=10,
                     context=dict(c11_dot16=c11_dot16, values=vect[:i].copy()))
    res.append(t)

cus_dot16 = pandas.DataFrame(res)
cus_dot16.tail()
average deviation max_exec min_exec number repeat size
75 0.009402 0.000889 0.010427 0.008099 50 10 187510
76 0.009694 0.000615 0.010664 0.008382 50 10 190010
77 0.011447 0.002160 0.015773 0.008372 50 10 192510
78 0.009392 0.000837 0.010663 0.008247 50 10 195010
79 0.010293 0.000732 0.011465 0.008879 50 10 197510
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0], label="numpy")
cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
ax[0].set_title("numpy and custom dot product execution time");
ax[1].set_title("numpy and custom dot product execution time");
../_images/cbenchmark_branching_45_0.png

We are far from numpy but the branching has clearly a huge impact and the fact the loop condition is evaluated only every 16 iterations does not explain this gain. Next experiment with SSE instructions.

Optimized to remove function call

We remove the function call to get the following version.

# float vector_dot_product_pointer16_nofcall(const float *p1, const float *p2, size_t size)
# {
#     float sum = 0;
#     const float * end = p1 + size;
#     if (size >= BYN) {
#         #if(BYN != 16)
#             #error "BYN must be equal to 16";
#         #endif
#         unsigned int size_ = (unsigned int) size;
#         size_ = size_ >> 4;  // division by 16=2^4
#         size_ = size_ << 4;  // multiplication by 16=2^4
#         const float * end_ = p1 + size_;
#         for(; p1 != end_;)
#         {
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#             sum += *p1 * *p2; ++p1, ++p2;
#         }
#     }
#     for(; p1 != end; ++p1, ++p2)
#         sum += *p1 * *p2;
#     return sum;
# }
#
# float vector_dot_product16_nofcall(py::array_t<float> v1, py::array_t<float> v2)
# {
#     if (v1.ndim() != v2.ndim())
#         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
#     if (v1.ndim() != 1)
#         throw std::runtime_error("Vector v1 and v2 must be vectors.");
#     return vector_dot_product_pointer16_nofcall(v1.data(0), v2.data(0), v1.shape(0));
# }
from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_nofcall
vector_dot_product16_nofcall(vect, vect)
3334629.0
def c11_dot16_nofcall(vect):
    return vector_dot_product16_nofcall(vect, vect)

measure_time("c11_dot16_nofcall(values)",
             context=dict(c11_dot16_nofcall=c11_dot16_nofcall, values=vect), repeat=10)
{'average': 1.1737164166,
 'deviation': 0.05164345150094166,
 'min_exec': 1.106978244000004,
 'max_exec': 1.26945232300001,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("c11_dot16_nofcall(values)", repeat=10,
                     context=dict(c11_dot16_nofcall=c11_dot16_nofcall, values=vect[:i].copy()))
    res.append(t)

cus_dot16_nofcall = pandas.DataFrame(res)
cus_dot16_nofcall.tail()
average deviation max_exec min_exec number repeat size
75 0.020536 0.001330 0.022946 0.018912 50 10 187510
76 0.020514 0.001460 0.022847 0.017647 50 10 190010
77 0.025768 0.003522 0.032687 0.020478 50 10 192510
78 0.022442 0.001405 0.024451 0.020280 50 10 195010
79 0.021800 0.001369 0.024049 0.019611 50 10 197510
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0], label="numpy")
cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
cus_dot16_nofcall.plot(x='size', y="average", ax=ax[0], label="pybind11x16_nofcall")
dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
cus_dot16_nofcall.plot(x='size', y="average", ax=ax[1], label="pybind11x16_nofcall")
ax[0].set_title("numpy and custom dot product execution time");
ax[1].set_title("numpy and custom dot product execution time");
../_images/cbenchmark_branching_52_0.png

Weird, branching did not happen when the code is not inside a separate function.

SSE instructions

We replace one function in the previous implementation.

# #include <xmmintrin.h>
#
# float vector_dot_product_pointer16_sse(const float *p1, const float *p2)
# {
#     __m128 c1 = _mm_load_ps(p1);
#     __m128 c2 = _mm_load_ps(p2);
#     __m128 r1 = _mm_mul_ps(c1, c2);
#
#     p1 += 4;
#     p2 += 4;
#
#     c1 = _mm_load_ps(p1);
#     c2 = _mm_load_ps(p2);
#     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
#
#     p1 += 4;
#     p2 += 4;
#
#     c1 = _mm_load_ps(p1);
#     c2 = _mm_load_ps(p2);
#     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
#
#     p1 += 4;
#     p2 += 4;
#
#     c1 = _mm_load_ps(p1);
#     c2 = _mm_load_ps(p2);
#     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
#
#     float r[4];
#     _mm_store_ps(r, r1);
#
#     return r[0] + r[1] + r[2] + r[3];
# }
from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_sse
vector_dot_product16_sse(vect, vect)
3333332.0
def c11_dot16_sse(vect):
    return vector_dot_product16_sse(vect, vect)

measure_time("c11_dot16_sse(values)", context=dict(c11_dot16_sse=c11_dot16_sse, values=vect), repeat=10)
{'average': 0.23317566440000093,
 'deviation': 0.0237807668204815,
 'min_exec': 0.21735401399999432,
 'max_exec': 0.29070227300000795,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("c11_dot16_sse(values)", repeat=10,
                     context=dict(c11_dot16_sse=c11_dot16_sse, values=vect[:i].copy()))
    res.append(t)

cus_dot16_sse = pandas.DataFrame(res)
cus_dot16_sse.tail()
average deviation max_exec min_exec number repeat size
75 0.003467 0.000436 0.003985 0.002719 50 10 187510
76 0.003366 0.000393 0.003866 0.002743 50 10 190010
77 0.003411 0.000406 0.003874 0.002591 50 10 192510
78 0.003413 0.000491 0.003938 0.002221 50 10 195010
79 0.003856 0.000463 0.004565 0.002685 50 10 197510
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0], label="numpy")
cus_dot16_sse.plot(x='size', y="average", ax=ax[0], label="pybind11x16_sse")
dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
cus_dot16_sse.plot(x='size', y="average", ax=ax[1], label="pybind11x16_sse")
cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
ax[0].set_title("numpy and custom dot product execution time");
ax[1].set_title("numpy and custom dot product execution time");
../_images/cbenchmark_branching_59_0.png

Better even though it is still slower than numpy. It is closer. Maybe the compilation option are not optimized, numpy was also compiled with the Intel compiler. To be accurate, multi-threading must be disabled on numpy side. That’s the purpose of the first two lines.

AVX 512

Last experiment with AVX 512 instructions but it does not work on all processor. I could not test it on my laptop as these instructions do not seem to be available. More can be found on wikipedia CPUs with AVX-512.

import platform
platform.processor()
'Intel64 Family 6 Model 78 Stepping 3, GenuineIntel'
import numpy
values = numpy.array(list(range(10000000)), dtype=numpy.float32)
vect = values / numpy.max(values)
from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_avx512
vector_dot_product16_avx512(vect, vect)
3333332.0
def c11_dot16_avx512(vect):
    return vector_dot_product16_avx512(vect, vect)

measure_time("c11_dot16_avx512(values)",
             context=dict(c11_dot16_avx512=c11_dot16_avx512, values=vect), repeat=10)
{'average': 0.2386972381999982,
 'deviation': 0.02440139257215092,
 'min_exec': 0.2189105550000079,
 'max_exec': 0.30276467900000625,
 'repeat': 10,
 'number': 50,
 'size': 10000000}
res = []
for i in range(10, 200000, 2500):
    t = measure_time("c11_dot16_avx512(values)", repeat=10,
                     context=dict(c11_dot16_avx512=c11_dot16_avx512, values=vect[:i].copy()))
    res.append(t)

cus_dot16_avx512 = pandas.DataFrame(res)
cus_dot16_avx512.tail()
average deviation max_exec min_exec number repeat size
75 0.003174 0.000449 0.003718 0.002469 50 10 187510
76 0.003114 0.000469 0.003789 0.002473 50 10 190010
77 0.003417 0.000593 0.003835 0.002502 50 10 192510
78 0.003976 0.001096 0.005927 0.002887 50 10 195010
79 0.003413 0.000465 0.003981 0.002700 50 10 197510
fig, ax = plt.subplots(1, 2, figsize=(14,4))
dot.plot(x='size', y="average", ax=ax[0], label="numpy")
cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
cus_dot16_sse.plot(x='size', y="average", ax=ax[0], label="pybind11x16_sse")
cus_dot16_avx512.plot(x='size', y="average", ax=ax[0], label="pybind11x16_avx512")
dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
cus_dot16_sse.plot(x='size', y="average", ax=ax[1], label="pybind11x16_sse")
cus_dot16_avx512.plot(x='size', y="average", ax=ax[1], label="pybind11x16_avx512")
cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
ax[0].set_title("numpy and custom dot product execution time");
ax[1].set_title("numpy and custom dot product execution time");
../_images/cbenchmark_branching_67_0.png

If the time is the same, it means that options AVX512 are not available.

from cpyquickhelper.numbers.cbenchmark_dot import get_simd_available_option
get_simd_available_option()
'Available options:  __SSE__ __SSE2__ __SSE3__ __SSE4_1__'

Back to numpy

This article Why is matrix multiplication faster with numpy than with ctypes in Python? gives some kints on why numpy is still faster. By looking at the code of the dot product in numpy: arraytypes.c.src, it seems that numpy does a simple dot product without using branching or uses the library BLAS which is the case in this benchmark (code for dot product: sdot.c). And it does use branching. See also function blas_stride. These libraries then play with compilation options and optimize for speed.