Мотивация
Совсем недавно вышла новая версия 0.34 библиотеки оптимизирующего JIT компилятора Numba для Python. И там ура! появилась долгожданная семантика аннотаций и набор методов для организации параллельных вычислений. За основу была взята технология Intel Parallel Accelerator.
В данной статье я хочу поделиться результатами первого тестирования скорости вычислений на основе этой библиотеки для некоторой современной машины с четырехядерным процессором.
Введение
В настоящее время Python очень активно используется в научных вычислениях, а в области Machine Learning вообще является практически одним из стандартов. Но если посмотреть чуть глубже, то почти везде Python используется как обертка над библиотеками более низкого уровня, написанных большей частью на C/C++. А можно ли на чистом Python писать на самом деле быстрый и параллельный код?
Рассмотрим совсем простую задачу. Пусть нам даны два набора N точек в трехмерном пространстве: p и q. Необходимо вычислить специальную матрицу на основе попарных расстояний между всеми точками:
Для всех тестов возьмем 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 мс |
Комментарии (35)
saw_tooth
29.08.2017 22:48-1Python интерпретируемый язык, Python использует внутри себя GIL, что делает невозможным параллелизацию на уровне самого кода
И по всей статье ни одного (даже плохого примера) с параллельностью. Есть и multiprocessing, и concurrent futures…
Да и тема «чистого» питона не раскрыта, можно генератор написать, можно map функцию использовать.
В общем посыл не понятен, что Вы хотели этой статьей показать, что питон медленный и не «параллельный»? Так мы это и так знаем)
alec_kalinin Автор
30.08.2017 12:29Хм… А вы точно читали статью, например, фразу:
"… заменить 'range' на 'prange' для параллельного выполнения внешнего цикла"?saw_tooth
30.08.2017 12:35Я читал статью, и имел ввиду «встроенную» параллельность, на встроенных модулях которые я перечислил.
saw_tooth
30.08.2017 12:44А можно ли на чистом Python писать на самом деле быстрый и параллельный код?
Вот с этого места.
не увидел «чистого» питона в статье, уж простите.
Возможно у нас просто разные взгляды на «чистый питон», тогда я уточню. Для меня это пайтон, который ставится из коробки и используется только функционал его встроенных модулей.alec_kalinin Автор
30.08.2017 12:59Ну хорошо. А давайте попробуем быть чуть ближе к практике? Вы пишите на чистом Python очень красивый web сервис. Но вот проблема, вам в рамках это сервиса нужно решить задачу, приведенную в статье. Это очень простая задача. И достаточно распространенная.
Что для вас наиболее чистое решение этой задачи?saw_tooth
30.08.2017 15:05+2Если касаться именно веб сервисов, то нужно разделять быстрые задачи и медленные. Вторые обычно кладутся с очередь задач и там исполняются, по окончанию пользователь может забрать результат.
Если же требуется решение в реальном времени, то задачи с недетерминированным временем в рантайм выносить нельзя. Сегодня у Вас 5000 точек, завтра будет 10к, после -1млн. Очевидно что человек, выполнивший такую операцию повесит сервер для всех.
Что касается реализации, о которой вы спросили, то лично я бы таки выносил это в очередь, и делал на concurrent futures, с очередью, что касается чистого пайтона. Если бы скорость меня не устроила, я бы написал расширение DLL/SO на С++/С и просто бы его использовал и пайтона. Возможно это не самый простой способ, ноя не люблю полагаться на сторонний код, если реализация подобного занимает немного времени.
Второй момент, я не придирался к вашей реализации, возможно в жизни она действительно хороша и эффективна, я даже не сомневаюсь в этом. Единственное что меня смутило, что Вы обозвали питон не параллельным из коробки, указали примеры на сторонних библиотеках, при этом не привели решения, который таки дает «коробочный» питон.
Uranix
30.08.2017 00:55+3В этой задаче немаловажным должен быть порядок хранения матриц. Почему для q в питоновском коде используется транспонированная форма?
iroln
30.08.2017 10:32Думаю, с транспонированием 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.)
alec_kalinin Автор
30.08.2017 12:23+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)
Сейчас я попробую сделать тест без транспонирования и погляжу на результаты.alec_kalinin Автор
30.08.2017 14:38+1После тестов время осталось примерно такое же. Но вообще, вы правы! Порядок очень важен, поскольку от него зависит попадание в кеш процессора.
Uranix
30.08.2017 15:44А еще я бы в с++ коде попробовал бы использовать Point4D вместо Point3D. В этом случае, мне кажется, выравнивание может еще ускорить код. Насколько я знаю, стандартные библиотеки для работы с матрицами стараются выравнивать строки.
Senpos
30.08.2017 11:38+1Не ожидал такого результата от Numba! Поразительно.
А не пробовали запустить тест на PyPy? Интересно посмотреть, что он может продемонстрировать в сравнении с другими участниками тестирования.
И если есть возможность, поделитесь кодом, который использовался для тестирования, я бы с удовольствием поигрался у себя на компьютере.
Спасибо за тест!alec_kalinin Автор
30.08.2017 14:37+1Да, я тоже не ожидал. Более того, Numba еще очень неплохо и прозрачно интегрируется с GPU и тоже, всего лишь за счет аннотаций "@cuda.jit". В общем, очень интересный проект.
Код я выложил:
github.com/alkalinin/open-nuance/tree/master/src-cpp/tests
github.com/alkalinin/open-nuance/tree/master/src-py/open_nuance/tests
Evgen52
30.08.2017 11:38Как насчёт того, чтобы попробовать сравнить еще и Intel Distribution for Python?
ogoNEKto
30.08.2017 11:51+1В таблице не хватает только теста на фортране с параллелизацией )))
x67
30.08.2017 16:14+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(*,*)ogoNEKto
30.08.2017 17:03Не придраться ради, а просто из спортивного интереса:
1) real(8) будет много шустрее работать
2) print обязательно выносить из-под цикла
Можно еще поиграть с опциями, я как-то «оптимизировал» код методом подбора параметров компиляции с 80+ секунд до 29…x67
30.08.2017 18:09- real(8) — плавающая точка с двойной точностью, такие операции выполняются медленнее. У меня вышло 96 мс.
- принт в данном случае выносить из цикла не надо, так как затраты на вывод информации не учитываются (считается время исполнения кода между двумя вызовами cpu_time)
А оптимизацией я заниматься не буду, так как оптимизируя опции компилятора вкупе с использованием старого кода, досконально не разбираясь в каждом привнесенном или унесенном эффекте, легко что-нибудь сломать (например, изменив опции qsave и qzero можно получить совсем неадекватные результаты, если в коде все переменные «инициализируются» implicit'ом), а в моем случае скорость работы кода не критична, но все должно считаться правильно)ogoNEKto
30.08.2017 18:47Про print принимается, в бенчмарк он не попадает, это точно)
А вот про замедление на real(8) это неожиданно, у меня он самый быстрый…x67
30.08.2017 19:36Он не может быть самым быстрым, так как всегда медленнее real 4. А если взять real(16), то разница в скорости будет уже на порядки по сравнению с real(8).
Надо смотреть условия, код, настройки компилятора, при которых real(8) самый быстрый, там точно что-то не так и это значит, что не реал8 быстрый, а реал4 слишком медленныйerwins22
31.08.2017 09:23Если SSE задействовано то реал(8) будет в 2 раза медленнее реал(4) так как обрабатывается в 2 раза больше данных за такт + еще больше вероятность попадания в кеш, так что больше 2х раз.
Psychopompe
30.08.2017 20:25Чистый С на одном ядре выдаёт 127мс без трюков с оптимизацией и -О3 флагом. Ирония в том, что количество расчётов можно уменьшить в ~два раза (N(N-1)/2), но время вычисления только возрастает, где-то процентов на 80 =\
Upd. все расчёты велись для double, float выдаёт 041мсPsychopompe
30.08.2017 20:34Чёрт, я про корень забыл, с ним время сразу скачет до Плеяд. Тогда получаются совсем другие результаты: float — 494мс без оптимизации, 250 — с ней самой, как и ожидалось. Смена float на double увеличивает эти значения на 10-15мс.
x67
30.08.2017 23:05новые, озвученные вами цифры характерны для одноядерных расчетов, вы видимо использовали чистый чистый си на одном ядре) У фортрана в таких условиях аналогичная производительность была.
Psychopompe
31.08.2017 02:17Да, чистый С, одно ядро :) Распараллелить такую задачу — не проблема, а вот если что-то более изощрённое, вроде barnes-hut или cell-linked list, то там беда.
evseev
30.08.2017 18:18На мой взгляд не хватает теста на Matlab. Все-таки его можно в данной области считать стандартом. К тому-же он достаточно быстрый. Ну вообще было-бы интересно.
Psychopompe
31.08.2017 02:18А чем статья отличается от этой?
Senpos
31.08.2017 07:21Тут обновилась Numba, реализация на которой свелась к красивому декоратору, а производительность значительно выросла.
evocatus
Я не знаю C++ вообще, но честное слово, код на C++ выглядит проще всех этих извращений.
А если вместо векторов использовать массивы, то думаю, что C++ будет работать ещё быстрее.
alec_kalinin Автор
Тут вот в чем дело. Python очень красивый язык с большим количеством красивых решений. Но вот иногда в рамках красивого и чистого Python проекта возникает вычислительная задача, где нужно порой нужно пройтись по 2D или 3D циклу и что-то вычислить. Кроме цикла сложно что-то придумать, и чистый Python для этого как раз работает плохо.
И тут как раз возникает ряд возможных решений с той иной степенью красоты:
И я как раз пробую насколько хорошо работает каждое из этих решений.
encyclopedist
А ещё есть PyPy, где JIT работает без доп информации о типах.
evocatus
А ещё можно использовать CFFI