Быстрый тест производительности Python для вычислительных задач

    Мотивация


    Совсем недавно вышла новая версия 0.34 библиотеки оптимизирующего JIT компилятора Numba для Python. И там ура! появилась долгожданная семантика аннотаций и набор методов для организации параллельных вычислений. За основу была взята технология Intel Parallel Accelerator.

    В данной статье я хочу поделиться результатами первого тестирования скорости вычислений на основе этой библиотеки для некоторой современной машины с четырехядерным процессором.

    Введение


    В настоящее время Python очень активно используется в научных вычислениях, а в области Machine Learning вообще является практически одним из стандартов. Но если посмотреть чуть глубже, то почти везде Python используется как обертка над библиотеками более низкого уровня, написанных большей частью на C/C++. А можно ли на чистом Python писать на самом деле быстрый и параллельный код?

    Рассмотрим совсем простую задачу. Пусть нам даны два набора N точек в трехмерном пространстве: p и q. Необходимо вычислить специальную матрицу на основе попарных расстояний между всеми точками:

    $R(i, j) = \frac{1}{1+\|p(i) - q(j)\|_2}$


    Для всех тестов возьмем N = 5000. Время вычисления усредняется для 10 запусков.

    Реализация на C++


    Как точку отсчета возьмем следующую реализацию на C++:

    void getR(std::vector<Point3D> & p, std::vector<Point3D> & q, Matrix & R)
    {
        double rx, ry, rz;
        
        #pragma omp parallel for
        for (int i = 0; i < p.size(); i++) {
            for (int j = 0; j < q.size(); j++) {
    	    rx = p[i].x - q[j].x;
                ry = p[i].y - q[j].y;
                rz = p[i].z - q[j].z;
    
                R(i, j) = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz));
            }
        }
    }
    

    Внешний цикл по точкам p выполняется параллельно с использованием технологии OpenMP.

    Время выполнения: 44 мс.

    Чистый Python


    Начнем тест скорости с кода на чистом Python:

    def get_R(p, q):
        R = np.empty((p.shape[0], q.shape[1]))
    
        for i in range(p.shape[0]):
            for j in range(q.shape[1]):
                rx = p[i, 0] - q[0, j]
                ry = p[i, 1] - q[1, j]
                rz = p[i, 2] - q[2, j]
                R[i, j] = 1 / (1 + math.sqrt(rx * rx + ry * ry + rz * rz))
    
        return R
    

    Время выполнения: 52 861 мс, медленнее базовой реализации больше, чем в 1000 раз.

    Python интерпретируемый язык, Python использует внутри себя GIL, что делает невозможным параллелизацию на уровне самого кода. Это все очень медленно. Теперь начнем это все ускорять.

    Python + NumPy + SciPy


    Проблему медленности Python для численных задач осознали очень давно. И ответом на эту проблему была библиотека NumPy. Идеология NumPy во многом близка MatLab, который является общепризнанным инструментом научных расчетов.

    Мы прекращаем мыслить итерационно, мы начинаем мыслить матрицами и векторами как атомарными объектами для вычисления. А все операции с матрицам и векторами на нижнем уровне уже выполняются высокопроизводительными библиотеками линейной алгебры Intel MKL или ATLAS.

    Реализация на NumPy выглядит так:

    def get_R_numpy(p, q):
        Rx = p[:, 0:1] - q[0:1]
        Ry = p[:, 1:2] - q[1:2]
        Rz = p[:, 2:3] - q[2:3]
    
        R = 1 / (1 + np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz))
    
        return R
    

    В этой реализации вообще нет ни одного цикла!

    Время выполнения: 839 мс, что медленнее базовой реализации где-то в 19 раз.

    Более того, в NumPy и SciPy есть огромное количество встроенных функций. Реализация данной задачи на SciPy выглядит так:

    def get_R_scipy(p, q):
        D = spd.cdist(p, q.T)
        R = 1 / (1 + D)
        return R
    

    Время выполнения: 353 мс, что медленнее базовой реализации в 8 раз.

    Для большинства задач это уже вполне приемлемое время работы. Ценой этого является переключение способа мышления, теперь необходимо собирать код из базовых операций линейной алгебры. Иногда это выглядит очень красиво, но порой приходится придумывать разные трюки.

    Но а как быть с параллельностью? Здесь она неявная. Мы надеемся, что на низком уровне все операции с матрицами и векторами реализованы эффективно и параллельно.

    А что делать, если наш код не вписывается в линейную алгебру, или мы хотим явную параллелизацию?

    Python + Cython


    Тут приходит время Cython. Cython это специальный язык, который позволяет внутри обычного Python кода вставлять код на C-образом языке. Далее Cython преобразует это код в .c файлы, которые компилируется в python модули. Эти модули достаточно прозрачно можно вызывать в других частях Python кода. Реализация на Cython выглядит так:

    @cython.wraparound(False)
    @cython.nonecheck(False)
    @cython.boundscheck(False)
    def get_R_cython_mp(py_p, py_q):
        py_R = np.empty((py_p.shape[0], py_q.shape[1]))
    
        cdef int nP = py_p.shape[0]
        cdef int nQ = py_q.shape[1]
        cdef double[:, :] p = py_p
        cdef double[:, :] q = py_q
        cdef double[:, :] R = py_R
    
        cdef int i, j
        cdef double rx, ry, rz
    
        with nogil:
            for i in prange(nP):
                for j in xrange(nQ):
                    rx = p[i, 0] - q[0, j]
                    ry = p[i, 1] - q[1, j]
                    rz = p[i, 2] - q[2, j]
    
                    R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
    
        return py_R
    

    Что тут происходит? На вход функция принимает python numpy объекты, далее они преобразуются в типизированные Cython С-структуры, а далее отключается gil и при помощи специальной конструкции 'prange' внешний цикл выполняется параллельно.

    Время выполнения: 76 мс, что в 1.75 раз медленнее, чем базовая реализация.

    Ну что, мы почти приблизились к базовой реализации, мы начали писать явный параллельный код. Но ценой этого стал менее читаемый код, мы ушли от чистого Python.

    В целом, большинство численных расчетов так и пишутся. Большая часть на NumPy, а некоторые места критичные по быстродействию выносятся в отдельные модули и реализуются на cython.

    Python + Numba


    Мы проделали длинный путь. Мы начали с чистого Python, потом пошли по пути магии матричных вычислений, потом погрузились в специальных язык расширений. Пора вернуться обратно к тому, с чего мы начали. Итак, реализация на Python + Numba:

    @jit(float64[:, :](float64[:, :], float64[:, :]), nopython=True, parallel=True)
    def get_R_numba_mp(p, q):
        R = np.empty((p.shape[0], q.shape[1]))
    
        for i in prange(p.shape[0]):
            for j in range(q.shape[1]):
                rx = p[i, 0] - q[0, j]
                ry = p[i, 1] - q[1, j]
                rz = p[i, 2] - q[2, j]
    
                R[i, j] = 1 / (1 + math.sqrt(rx * rx + ry * ry + rz * rz))
    
        return R
    

    Время выполнения: 46 мс, что практически совпадает с базовой реализацией!

    И все, что нужно было сделать для этого с исходным медленным Python кодом:

    • добавить аннотацию с указанием типов и параллельным типом выполнения
    • заменить 'range' на 'prange' для параллельного выполнения внешнего цикла

    Numba это очень интересный проект. Это оптимизирующий JIT компилятор для Python на основе LLVM. Он абсолютно прозрачен для использования, не нужно никаких отдельных модулей. Все что нужно, это добавить несколько аннотаций в критические по скорости методы.

    В итоге, времена исполнения следующие:
    C++ Python Numpy SciPy Cython Numba
    44 мс 52 861 мс 839 мс 353 мс 76 мс 46 мс
    Поделиться публикацией
    Реклама помогает поддерживать и развивать наши сервисы

    Подробнее
    Реклама
    Комментарии 38
    • 0
      Я не знаю C++ вообще, но честное слово, код на C++ выглядит проще всех этих извращений.
      А если вместо векторов использовать массивы, то думаю, что C++ будет работать ещё быстрее.
      • 0
        Тут вот в чем дело. Python очень красивый язык с большим количеством красивых решений. Но вот иногда в рамках красивого и чистого Python проекта возникает вычислительная задача, где нужно порой нужно пройтись по 2D или 3D циклу и что-то вычислить. Кроме цикла сложно что-то придумать, и чистый Python для этого как раз работает плохо.

        И тут как раз возникает ряд возможных решений с той иной степенью красоты:
        • отказаться от итераций и работать только с матричными операциями (NumPy)
        • свести задачу к набору стандартных базовых операций (SciPy)
        • спуститься к уровню C/C++ (и только здесь возникает С++ код)
        • добавить информацию о типах для JIT компилятора

        И я как раз пробую насколько хорошо работает каждое из этих решений.
        • 0

          А ещё есть PyPy, где JIT работает без доп информации о типах.

          • 0
            А ещё можно использовать CFFI
        • –1
          Python интерпретируемый язык, Python использует внутри себя GIL, что делает невозможным параллелизацию на уровне самого кода

          И по всей статье ни одного (даже плохого примера) с параллельностью. Есть и multiprocessing, и concurrent futures…
          Да и тема «чистого» питона не раскрыта, можно генератор написать, можно map функцию использовать.
          В общем посыл не понятен, что Вы хотели этой статьей показать, что питон медленный и не «параллельный»? Так мы это и так знаем)
          • 0
            Хм… А вы точно читали статью, например, фразу:
            "… заменить 'range' на 'prange' для параллельного выполнения внешнего цикла"?
            • 0
              Я читал статью, и имел ввиду «встроенную» параллельность, на встроенных модулях которые я перечислил.
              • 0
                А можно ли на чистом Python писать на самом деле быстрый и параллельный код?

                Вот с этого места.
                не увидел «чистого» питона в статье, уж простите.
                Возможно у нас просто разные взгляды на «чистый питон», тогда я уточню. Для меня это пайтон, который ставится из коробки и используется только функционал его встроенных модулей.
                • 0
                  Ну хорошо. А давайте попробуем быть чуть ближе к практике? Вы пишите на чистом Python очень красивый web сервис. Но вот проблема, вам в рамках это сервиса нужно решить задачу, приведенную в статье. Это очень простая задача. И достаточно распространенная.

                  Что для вас наиболее чистое решение этой задачи?
                  • +2
                    Если касаться именно веб сервисов, то нужно разделять быстрые задачи и медленные. Вторые обычно кладутся с очередь задач и там исполняются, по окончанию пользователь может забрать результат.
                    Если же требуется решение в реальном времени, то задачи с недетерминированным временем в рантайм выносить нельзя. Сегодня у Вас 5000 точек, завтра будет 10к, после -1млн. Очевидно что человек, выполнивший такую операцию повесит сервер для всех.
                    Что касается реализации, о которой вы спросили, то лично я бы таки выносил это в очередь, и делал на concurrent futures, с очередью, что касается чистого пайтона. Если бы скорость меня не устроила, я бы написал расширение DLL/SO на С++/С и просто бы его использовал и пайтона. Возможно это не самый простой способ, ноя не люблю полагаться на сторонний код, если реализация подобного занимает немного времени.
                    Второй момент, я не придирался к вашей реализации, возможно в жизни она действительно хороша и эффективна, я даже не сомневаюсь в этом. Единственное что меня смутило, что Вы обозвали питон не параллельным из коробки, указали примеры на сторонних библиотеках, при этом не привели решения, который таки дает «коробочный» питон.
          • +3
            В этой задаче немаловажным должен быть порядок хранения матриц. Почему для q в питоновском коде используется транспонированная форма?
            • 0

              Думаю, с транспонированием q что-то не так. Зачем там транспонирование — не понятно. cdist принимает на вход данные в виде такой матрицы:


              [(x1, y1, z1, ..., n1),
               (x2, y2, z2, ..., n2),
               ...
               (xK, yK, zK, ..., nK)
              ]

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


              ValueError: XA and XB must have the same number of columns (i.e. feature dimension.)
              • +1
                Да, тут вы правы, спасибо за замечание! Тест вырос из NumPy кода, а там я использую такое хранение для того, чтобы избежать итераций и сразу получить матрицу за счет использования numpy broadcasing rules, например:

                import numpy as np
                
                N = 5000
                
                p = np.random.rand(N, 1)
                q = np.random.rand(N, 1)
                
                # no broadcasting
                print((p - q).shape)
                
                >>> (5000, 1)
                
                # with broadcasting
                print((p - q.T).shape)
                
                >>> (5000, 5000)
                


                Сейчас я попробую сделать тест без транспонирования и погляжу на результаты.
                • +1
                  После тестов время осталось примерно такое же. Но вообще, вы правы! Порядок очень важен, поскольку от него зависит попадание в кеш процессора.
                  • 0
                    А еще я бы в с++ коде попробовал бы использовать Point4D вместо Point3D. В этом случае, мне кажется, выравнивание может еще ускорить код. Насколько я знаю, стандартные библиотеки для работы с матрицами стараются выравнивать строки.
              • +1
                Не ожидал такого результата от Numba! Поразительно.

                А не пробовали запустить тест на PyPy? Интересно посмотреть, что он может продемонстрировать в сравнении с другими участниками тестирования.

                И если есть возможность, поделитесь кодом, который использовался для тестирования, я бы с удовольствием поигрался у себя на компьютере.

                Спасибо за тест!
              • 0
                Жаль нет pypy.
                • 0

                  Как насчёт того, чтобы попробовать сравнить еще и Intel Distribution for Python?

                  • +1
                    В таблице не хватает только теста на фортране с параллелизацией )))
                    • +1
                      У меня на интеловском компиляторе код ниже(переписал исходники автора) с O3 оптимизацией выходит 4.43E-2 секунды в среднем. Навряд ли у кого-то с 8 ядрами выйдет сильно быстрее.
                          USE IFPORT
                          REAL(4) p(5000,3),q(5000,3),R(5000,5000)
                          integer i,o
                          CALL SRAND(100)
                          do j=1,1000
                             N = 5000
                             do i=1,N
                                 do o=1,3
                                     p(i,o)=RAND()
                                     q(i,o)=RAND()
                                 end do
                             end do
                              CALL CPU_time(t1)
                              DO CONCURRENT  (i=1:size(p,1))
                                  DO CONCURRENT  (o=1:size(q,1))
                                      rx=p(i,1)-q(o,1)
                                      ry=p(i,2)-q(o,2)
                                      rz=p(i,3)-q(o,3)
                                      R(i,o)=1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
                                  end do
                              end do
                              CALL CPU_time(t2)
                              print *, t1, t2, t2-t1
                              dt=dt+(t2-t1)
                              print *, R(1,1) !чтобы компилятор не халтурил и выполнял всю работу
                          end do
                          print *, 'Avg time:', dt/1000.
                          read(*,*)
                      • 0
                        Не придраться ради, а просто из спортивного интереса:
                        1) real(8) будет много шустрее работать
                        2) print обязательно выносить из-под цикла
                        Можно еще поиграть с опциями, я как-то «оптимизировал» код методом подбора параметров компиляции с 80+ секунд до 29…
                        • 0
                          • real(8) — плавающая точка с двойной точностью, такие операции выполняются медленнее. У меня вышло 96 мс.
                          • принт в данном случае выносить из цикла не надо, так как затраты на вывод информации не учитываются (считается время исполнения кода между двумя вызовами cpu_time)

                          А оптимизацией я заниматься не буду, так как оптимизируя опции компилятора вкупе с использованием старого кода, досконально не разбираясь в каждом привнесенном или унесенном эффекте, легко что-нибудь сломать (например, изменив опции qsave и qzero можно получить совсем неадекватные результаты, если в коде все переменные «инициализируются» implicit'ом), а в моем случае скорость работы кода не критична, но все должно считаться правильно)
                          • 0
                            Про print принимается, в бенчмарк он не попадает, это точно)
                            А вот про замедление на real(8) это неожиданно, у меня он самый быстрый…
                            • 0
                              Он не может быть самым быстрым, так как всегда медленнее real 4. А если взять real(16), то разница в скорости будет уже на порядки по сравнению с real(8).
                              Надо смотреть условия, код, настройки компилятора, при которых real(8) самый быстрый, там точно что-то не так и это значит, что не реал8 быстрый, а реал4 слишком медленный
                              • 0
                                Если SSE задействовано то реал(8) будет в 2 раза медленнее реал(4) так как обрабатывается в 2 раза больше данных за такт + еще больше вероятность попадания в кеш, так что больше 2х раз.
                        • 0
                          Чистый С на одном ядре выдаёт 127мс без трюков с оптимизацией и -О3 флагом. Ирония в том, что количество расчётов можно уменьшить в ~два раза (N(N-1)/2), но время вычисления только возрастает, где-то процентов на 80 =\
                          Upd. все расчёты велись для double, float выдаёт 041мс
                          • 0
                            Чёрт, я про корень забыл, с ним время сразу скачет до Плеяд. Тогда получаются совсем другие результаты: float — 494мс без оптимизации, 250 — с ней самой, как и ожидалось. Смена float на double увеличивает эти значения на 10-15мс.
                            • 0
                              новые, озвученные вами цифры характерны для одноядерных расчетов, вы видимо использовали чистый чистый си на одном ядре) У фортрана в таких условиях аналогичная производительность была.
                              • 0
                                Да, чистый С, одно ядро :) Распараллелить такую задачу — не проблема, а вот если что-то более изощрённое, вроде barnes-hut или cell-linked list, то там беда.
                          • 0
                            .
                            • 0
                              На мой взгляд не хватает теста на Matlab. Все-таки его можно в данной области считать стандартом. К тому-же он достаточно быстрый. Ну вообще было-бы интересно.
                            • 0
                              А чем статья отличается от этой?
                              • 0
                                Тут обновилась Numba, реализация на которой свелась к красивому декоратору, а производительность значительно выросла.
                              • 0

                                Я конечно тоже считаю numpy неотъемлемой частью питона и не буду городить списки списков для двумерных массивов чисел, но всё-таки код на чистом питоне будет выглядеть как-то так:


                                def get_R_pure_python(p, q):
                                    R = [[0 for x in range(len(q))] for y in range(len(p))]
                                
                                    for i in range(len(p)):
                                        for j in range(len(q)):
                                            rx = p[i][0] - q[j][0]
                                            ry = p[i][1] - q[j][1]
                                            rz = p[i][2] - q[j][2]
                                            R[i][j] = 1 / (1 + math.sqrt(rx * rx + ry * ry + rz * rz))
                                
                                    return R

                                Кроме того, этот код быстрее чем numpy без broadcasting и его уже намного проще подсунуть pypy. В итоге на i5-2500 я получил следующие результаты:


                                | numpy 1.13.1 (без broadcasting) | 42532 мс | 477.89x |
                                | python 3.5.2 | 22247 мс | 249.97x |
                                | python 2.7.12 | 17773 мс | 199.70x |
                                | pypy 5.1.2 | 621 мс | 6.98x |
                                | numpy 1.13.1 (с broadcasting) | 511 мс | 5.74x |
                                | scipy 0.19.1 | 280 мс | 3.15x |
                                | numba 0.35.0rc1 | 89 мс | 1x |


                                Для python 2.7.12 функцию range я заменил на xrange там где это надо.
                                P.S. У кого-нибудь есть возможность поправить комментарий, чтобы заработал Markdown, потому что у меня, к сожалению, не работает.

                                • 0
                                  Не знал, что Markdown начинает работать только через какое-то время.
                                • 0
                                  Начнем тест скорости с кода на чистом Python:

                                      R = np.empty((p.shape[0], q.shape[1]))

                                  это разве не Numpy методы?
                                  т.е. это уже не совсем чистый python?
                                  /занудамодеоф

                                  Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.