Numpy быстрее, чем Numba и Cython, как улучшить код Numba

Question

У меня есть простой пример, который поможет мне понять использование numba и cython. Я новичок как в numba, так и в cython. Я старался изо всех сил, чтобы включить все приемы, чтобы сделать Numba быстрой и в некоторой степени такой же для Cython, но мой код NumPy почти в 2 раза быстрее, чем Numba (для float64), более чем в 2 раза быстрее, если использовать float32. Не уверен, что мне здесь не хватает.

Я подумал, что, возможно, проблема не в кодировании, а в компиляторе и в том, с чем я не очень знаком.

Я прочитал много статей о переполнении стека о numpy, numba и cython и не нашел прямых ответов.

NumPy версия:

def py_expsum(x):
    return np.sum( np.exp(x) )

версия numba:

@numba.jit( nopython=True)    
def nb_expsum(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in range(ny):
            val += np.exp(x[ix, iy])
    return val

Версия Cython:

import numpy as np
import cython
from libc.math cimport exp

@cython.boundscheck(False) 
@cython.wraparound(False)
cpdef double cy_expsum2 ( double[:,:] x, int nx, int ny ):
    cdef: 
        double val = 0.0
        int ix, iy    
    for ix in range(nx):
        for iy in range(ny):
            val += exp(x[ix, iy])
    return val

играть с массивом размером 2000 x 1000 и повторять циклы более 100 раз. Для numba, первый раз, когда он активирован, не учитывается в цикле.

Использование python 3 (дистрибутив anaconda), окно 10

               float64       /   float32
    1. numpy : 0.56 sec      /   0.23 sec
    2. numba : 0.93 sec      /   0.74 sec      
    3. cython: 0.83 sec

Цитон близок к Нумбе. Поэтому большой вопрос для меня заключается в том, почему numba не может побить время выполнения numpy? Что я сделал не так или упустил здесь? Как могут влиять другие факторы и как мне это выяснить?

Answer 1

Как мы увидим, поведение зависит от того, какое numpy-распределение используется.

В этом ответе основное внимание будет уделено распространению Anacoda с использованием VML-библиотеки Intel (векторной математической библиотеки), в зависимости от другого оборудования и numpy-версии.

Также будет показано, как VML можно использовать через Cython или numexpr, если не используется Anacoda-дистрибутив, который подключает VML под капотом для некоторого numpy- операции.

---

Я могу воспроизвести ваши результаты, для следующих измерений

N,M=2*10**4, 10**3
a=np.random.rand(N, M)

Я получаю:

%timeit py_expsum(a)  #   87ms
%timeit nb_expsum(a)  #  672ms
%timeit nb_expsum2(a)  #  412ms

Львиная доля (около 90%) времени вычислений используется для оценки exp - функции, и, как мы увидим, это задача, интенсивно использующая ЦП.

Быстрый взгляд на top -статистику показывает, что версия numpy выполнена парализованной, но это не относится к numba. Однако на моей виртуальной машине, имеющей только два процессора, параллелизация сама по себе не может объяснить огромную разницу фактора 7 (как показывает версия DavidW nb_expsum2).

Профилирование кода с помощью perf для обеих версий показывает следующее:

nb_expsum

Overhead  Command  Shared Object                                      Symbol                                                             
  62,56%  python   libm-2.23.so                                       [.] __ieee754_exp_avx
  16,16%  python   libm-2.23.so                                       [.] __GI___exp
   5,25%  python   perf-28936.map                                     [.] 0x00007f1658d53213
   2,21%  python   mtrand.cpython-37m-x86_64-linux-gnu.so             [.] rk_random

py_expsum

  31,84%  python   libmkl_vml_avx.so                                  [.] mkl_vml_kernel_dExp_E9HAynn                                   ▒
   9,47%  python   libiomp5.so                                        [.] _INTERNAL_25_______src_kmp_barrier_cpp_38a91946::__kmp_wait_te▒
   6,21%  python   [unknown]                                          [k] 0xffffffff8140290c                                            ▒
   5,27%  python   mtrand.cpython-37m-x86_64-linux-gnu.so             [.] rk_random  

Как видно: numpy использует парализованную векторизованную mkl / vml-версию Intel под капотом, которая легко превосходит версию из библиотеки gnu-math (lm.so), используемой numba (или параллельной версией numba или на Cython в этом отношении). Можно немного выровнять почву, используя парализатор, но все же векторизованная версия mkl превзойдет numba и cython.

Тем не менее, видение производительности только для одного размера не очень показательно, и в случае exp (как и для других трансцендентных функций) необходимо учитывать 2 аспекта:

• количество элементов в массиве - эффекты кэширования и разные алгоритмы для разных размеров (не случайно в numpy) могут привести к разной производительности.
• в зависимости от значения x, для расчета exp(x) требуется разное время. Обычно существует три различных типа ввода, приводящих к разным временам расчета: очень маленький, нормальный и очень большой (с не конечными результатами)

Я использую perfplot для визуализации результата (см. Код в приложении). Для «нормального» диапазона мы получаем следующие показатели:

и хотя производительность для 0.0 аналогична, мы можем видеть, что VML от Intel получает весьма негативное влияние, как только результаты становятся бесконечными:

Однако есть и другие вещи, на которые стоит обратить внимание:

• Для векторных размеров <= 8192 = 2^13 numpy использует непараллельную glibc-версию exp (используются те же numba и cython).
• Anaconda-дистрибутив, который я использую, переопределяет функциональность numpy и подключает VML-библиотеку Intel для размеров> 8192, которая векторизована и распараллелена - это объясняет сокращение времени выполнения для размеров около 10 ^ 4 .
• numba легко превосходит обычную glibc-версию (слишком много накладных расходов для numpy) для меньших размеров, но будет (если numpy не переключится на VML) небольшая разница для большого массива.
• Кажется, это задача, связанная с процессором - мы нигде не видим границ кэша.
• Парализованная numba-версия имеет смысл только при наличии более 500 элементов.

Так каковы последствия?

1. Если элементов не более 8192, следует использовать numba-версию.
2. в противном случае версия NumPy (даже если VML-плагин не доступен, он не сильно потеряет).

Примечание: numba не может автоматически использовать vdExp из VML-кода Intel (как это частично предлагается в комментариях), поскольку он вычисляет exp(x) индивидуально, тогда как VML работает с целым массивом.

---

Можно уменьшить количество кешей при записи и загрузке данных, что выполняется в numpy-версии с использованием следующего алгоритма:

1. Выполнение VML vdExp для части данных, которая соответствует кешу, но также не слишком мала (накладные расходы).
2. Суммируйте полученный рабочий массив.
3. Выполнить 1. + 2. для следующей части данных, пока не будут обработаны все данные.

Однако я бы не ожидал, что выиграю больше, чем на 10% (но, возможно, я ошибаюсь) по сравнению с версией numpy, поскольку 90% времени вычислений в любом случае тратится на MVL.

Тем не менее, вот возможная быстрая и грязная реализация в Cython:

%%cython -L=<path_mkl_libs> --link-args=-Wl,-rpath=<path_mkl_libs> --link-args=-Wl,--no-as-needed -l=mkl_intel_ilp64 -l=mkl_core -l=mkl_gnu_thread -l=iomp5
# path to mkl can be found via np.show_config()
# which libraries needed: https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor

# another option would be to wrap mkl.h:
cdef extern from *:
    """
    // MKL_INT is 64bit integer for mkl-ilp64
    // see https://software.intel.com/en-us/mkl-developer-reference-c-c-datatypes-specific-to-intel-mkl
    #define MKL_INT long long int
    void  vdExp(MKL_INT n, const double *x, double *y);
    """
    void vdExp(long long int n, const double *x, double *y)

def cy_expsum(const double[:,:] v):
        cdef:
            double[1024] w;
            int n = v.size
            int current = 0;
            double res = 0.0
            int size = 0
            int i = 0
        while current<n:
            size = n-current
            if size>1024:
                size = 1024
            vdExp(size, &v[0,0]+current, w)
            for i in range(size):
                res+=w[i]
            current+=size
        return res

Однако это именно то, что сделал бы numexpr, который также использует Intel vml в качестве бэкэнда:

 import numexpr as ne
 def ne_expsum(x):
     return ne.evaluate("sum(exp(x))")

Что касается времени, мы можем видеть следующее:

со следующими примечательными деталями:

• numpy, Numberxpr и Cython версии имеют почти одинаковую производительность для больших массивов - что неудивительно, поскольку они используют одинаковую vml-функциональность.
• из этих трех, Cython-версия имеет наименьшие накладные расходы и цифру xpr больше всего
• Numberxpr-версия, вероятно, самая простая для написания (учитывая, что не все numpy дистрибутивы подключаются к mvl-функциональности).

---

Тэг:

Сюжеты:

import numpy as np
def py_expsum(x):
    return np.sum(np.exp(x))

import numba as nb
@nb.jit( nopython=True)    
def nb_expsum(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in range(ny):
            val += np.exp( x[ix, iy] )
    return val

@nb.jit( nopython=True, parallel=True)    
def nb_expsum2(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in nb.prange(ny):
            val += np.exp( x[ix, iy]   )
    return val

import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
    setup=lambda n: factor*np.random.rand(1,n),
    n_range=[2**k for k in range(0,27)],
    kernels=[
        py_expsum, 
        nb_expsum,
        nb_expsum2, 
        ],
    logx=True,
    logy=True,
    xlabel='len(x)'
    )

Answer 2

Добавить распараллеливание. В Numba это просто включает создание внешнего цикла prange и добавление parallel=True к параметрам jit:

@numba.jit( nopython=True,parallel=True)    
def nb_expsum2(x):
    nx, ny = x.shape
    val = 0.0
    for ix in numba.prange(nx):
        for iy in range(ny):
            val += np.exp( x[ix, iy]   )
    return val

На моем ПК это дает ускорение в 3,2 раза по сравнению с непараллельной версией. Тем не менее, на моем ПК и Numba, и Cython победили Numpy, как написано.

Вы также можете выполнить распараллеливание в Cython - я не тестировал здесь, но я бы ожидал, что он будет похож на Numba по производительности. (Обратите внимание также, что для Cython вы можете получить nx и ny от x.shape[0] и x.shape[1], поэтому вам не нужно отключать проверку границ, а затем полностью полагаться на пользовательский ввод, чтобы не выходить за границы).

Answer 3

Это зависит от реализации exp и распараллеливания

Если вы используете Intel SVML в Numpy, используйте его и в других пакетах, таких как Numba, Numexpr или Cython. Советы по производительности Numba

Если команды Numpy распараллелены, попробуйте распараллелить их в Numba или Cython.

Код

import os
#Have to be before importing numpy
#Test with 1 Thread against a single thread Numba/Cython Version and
#at least with number of physical cores against parallel versions
os.environ["MKL_NUM_THREADS"] = "1" 

import numpy as np
import numba as nb

def py_expsum(x):
    return np.sum( np.exp(x) )

@nb.njit(parallel=False) #set it to True for a parallel version  
def nb_expsum(x):
    val = 0.
    for ix in nb.prange(x.shape[0]):
        for iy in range(x.shape[1]):
            val += np.exp(x[ix,iy])
    return val

N,M=2000, 1000
#a=np.random.rand(N*M).reshape((N,M)).astype(np.float32)
a=np.random.rand(N*M).reshape((N,M))

Тесты

#float64
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1" 
8.79 ms ± 85.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6" 
5.78 ms ± 131 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
9.51 ms ± 89.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=True
1.77 ms ± 136 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

#float32
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1" 
4.08 ms ± 86 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6" 
3.09 ms ± 335 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
7.35 ms ± 38 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=True
1.29 ms ± 31 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Очевидно, что у Numba, возможно, из-за автоматической системы набора текста, есть некоторые проблемы с float32, но параллельная версия в обоих случаях быстрее, чем версия Numpy.В обоих случаях использовался SVML.