Лекция 9. Оптимизация и параллельные вычисления

Оптимизация и параллельные вычисленияПоиск места, где программа теряет время, и ускорение расчёта на Python: профилирование, приёмы языка, память, Numba, Cython, NumPy, потоки, процессы и видеокарта
Вячеслав Федоров
Лекция 9
Разработка и применение ПО
в физических исследованиях
100x
Изложение. Девятая лекция посвящена скорости вычислений на Python. Язык удобен, но интерпретация, динамическая типизация и объектная модель обходятся дорого: наивный цикл на Python выполняется в десятки и сотни раз медленнее того же кода на C. Для скрипта, запускаемого раз в день, это несущественно, а для обработки данных с детектора или моделирования миллиона частиц разница составляет часы против недель. В рамках лекции рассматриваются причины медленности, измерение времени и профилирование, приёмы, не требующие ничего, кроме самого языка, экономия памяти, компиляция через Numba и Cython, векторизация NumPy и параллельные вычисления: потоки и глобальная блокировка интерпретатора, пул процессов, сборка Python без GIL и видеокарта. Надпись на титуле, 100x, обозначает порядок ускорения, которого позволяет достичь компиляция горячего цикла; параллельные вычисления на процессоре добавляют к нему множитель, не превышающий числа ядер.
Демонстрации. Слайды с розовым выделением слова «Демонстрация» отмечают места, где изложение прерывается для работы в терминале виртуальной машины lab. На каждом таком слайде проигрывается запись тех же команд; числа в записях получены на этой виртуальной машине с двумя ядрами и Python 3.12, а в седьмой части — также со свободнопоточным Python 3.14; далее машина называется стендом. На другой машине иными будут и сами числа, и соотношения между ними: выигрыш параллельных версий зависит от числа ядер, а выигрыш BLAS — от процессора и библиотеки. Видеокарты на стенде нет, поэтому ядро CUDA выполняется в симуляторе на процессоре и служит только для проверки индексов, а не для замера скорости.

СодержаниеОт замера времени до компиляции и параллельных вычислений
01
Причины медленности
02
Измерение
03
Приёмы языка
04
Память
05
Компиляция
06
NumPy и BLAS
07
Потоки, процессы и GPU
Изложение. Лекция состоит из семи частей. Первая объясняет, откуда берётся медленность Python: динамическая типизация, изменяемость всего во время выполнения, устройство интерпретатора. Вторая посвящена измерению: когда оптимизация оправдана, как замерить время фрагмента модулем timeit и найти горячее место профилировщиками cProfile, line_profiler и py-spy. Третья собирает приёмы самого языка с замерами на стенде. Четвёртая — память: как её измерить и как хранить данные компактнее. Пятая — компиляция горячих циклов через Numba и Cython на примерах умножения матриц и пространственного заряда пучка. Шестая — векторизация NumPy и библиотека BLAS. Седьмая — параллельные вычисления: потоки и глобальная блокировка интерпретатора, пул процессов, сборка Python без GIL и модель вычислений CUDA на видеокарте.

Цель лекцииНайти место, где программа теряет время, и ускорить именно его
Профилировщик показывает, на что уходит время. Далее выбирается самый дешёвый способ ускорения: другой алгоритм, приём языка, NumPy, компиляция горячего цикла или параллельные вычисления. Каждый шаг подтверждается замером и тестом: ускоренный неверный расчёт остаётся неверным.
Изложение. Python обладает многими достоинствами, кроме скорости. Медленный код на Python почти всегда можно ускорить, не переписывая его целиком на C++, но для этого необходимо знать, где именно он медленный. Интуитивные представления о том, где программа проводит время, как правило, ошибочны: функция, которая кажется наиболее затратной, часто занимает единицы процентов, а основное время уходит в незаметный цикл. Таким образом, порядок работы всегда один и тот же: измерить, найти горячее место, ускорить его наименее затратным способом и снова измерить.
Для физиков: расчёт динамики пучка или обработка данных эксперимента обычно проводят большую часть времени в одном-двух вложенных циклах. Ускорение именно этих циклов в сто раз превращает расчёт, не укладывающийся в ночь, в расчёт на несколько минут.

Порядок
ускорения
01
Измерение
профилировщик показывает горячее место
02
Алгоритм и структура данных
O(n²) → O(n log n), список → множество
03
Приёмы языка
встроенные функции, локальные имена, идиомы
04
NumPy
цикл уходит в C и BLAS
05
Numba или Cython
компиляция горячего цикла в машинный код
06
Потоки, процессы и GPU
все ядра процессора и видеокарта
Важно
Сначала правильно, потом быстро: каждая правка проверяется тестом и замером на тех же данных.
Для физика
Порядок выигрыша на стенде: приёмы языка — в несколько раз, NumPy и компиляция — в десятки и сотни раз, процессы — не больше числа ядер.
Пояснение
Шаги идут по возрастанию цены: дальше вниз — больше зависимостей и сложнее код.
Изложение. Способы ускорения удобно выстроить по возрастанию цены. Сначала измерение: без профиля неизвестно, что ускорять. Затем алгоритм и структура данных: замена квадратичного алгоритма на n log n или списка на множество меняет асимптотику, чего никакая компиляция не обеспечивает. Приёмы самого языка дают выигрыш от нескольких процентов до нескольких раз и не требуют дополнительных затрат. NumPy переносит цикл в скомпилированный код и библиотеку линейной алгебры BLAS. Если задача не векторизуется, горячий цикл компилируют в машинный код через Numba или Cython. Последняя ступень — параллельные вычисления: потоки, процессы и сборка без GIL на нескольких ядрах или видеокарта. Чем ниже ступень, тем больше зависимостей и тем сложнее код, поэтому переходить к следующему шагу имеет смысл, только пока расчёт не укладывается в требуемое время.

Причины медленностиДинамическая типизация, изменяемость, устройство интерпретатора CPython
01
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. В первой части рассмотрим, откуда берётся разница со скоростью C. Гибкость Python препятствует оптимизациям: всякая оптимизация опирается на предположения, а интерпретатору почти ничего не позволено считать заранее известным.

Динамическая
типизация
def add(a, b):
    return a + b
 
add(1, 2)          # 3
add(1.5, 2)        # 3.5
add('ab', 'c')     # 'abc'
add([1], [2])      # [1, 2]

PythonC
тип значенияизвестен при выполненииизвестен компилятору
числообъект в куче, int — 28 байт8 байт в регистре
a + bвыбор операции по типамодна инструкция
Пояснение
Интерпретатор проверяет типы на каждой операции; заранее вычисляются лишь выражения из констант, остальной код исполняется буквально.
Изложение. Тип значения, связанного с именем, становится известен только во время выполнения: одна и та же функция add складывает целые и вещественные числа, склеивает строки и списки. Поэтому на каждой операции интерпретатор выясняет типы операндов и выбирает подходящую реализацию сложения, а эти проверки, повторяющиеся миллионы раз, требуют времени. Каждое число в Python является объектом в куче со счётчиком ссылок и указателем на тип: целое число занимает 28 байт, тогда как в C переменная типа double — 8 байт, а сложение — одна инструкция процессора. Кроме того, интерпретатор не знает заранее, с какими данными работает, и обязан исполнять написанный код почти буквально: компилятор заранее вычисляет лишь выражения из литералов, например 60 * 60 * 24, и отбрасывает ветвь if False, а выражение с именами вычислить заранее нельзя, поскольку любое имя может быть переопределено.

	Python	C
тип значения	известен при выполнении	известен компилятору
число	объект в куче, int — 28 байт	8 байт в регистре
a + b	выбор операции по типам	одна инструкция

Изменяемость
всего
import builtins
 
len = lambda obj: 'mock!'
len('abc')         # 'mock!'
len = builtins.len
 
def f(a, b): return a + b
def g(a, b): return a * b
f.__code__ = g.__code__
f(1, 2)            # 2
Во время выполнения подменяются встроенные имена, тело функции, атрибуты классов. Интерпретатор не вправе считать, что len — это len, а range — это range.
for i in range(len(a)):
    if i == 0:         # на каждом шаге
        a[i] = 1
    else:
        a[i] = i
Пояснение
Проверку, которую компилятор C вынес бы из цикла, переносит только программист: 13,1 → 9,8 мкс.
Изложение. В Python почти всё может быть изменено во время выполнения: встроенное имя len переопределяется присваиванием, тело уже определённой функции заменяется через атрибут __code__, атрибуты классов и модулей меняются в любой момент. Поэтому интерпретатор не может заранее подставить тело функции в место вызова или отделить от цикла первую итерацию, на которой проверка i == 0 истинна: он не знает, не будут ли по ходу выполнения переопределены имена range и len. Оптимизацию, разрешённую любому компилятору C, приходится выполнять вручную. Покажем это на примере из главы книги: проверка i == 0 выполняется на каждой из тысячи итераций, а вынесенная из цикла — один раз, и на стенде функция ускоряется с 13,1 до 9,8 мкс, на четверть.

Интерпретатор
CPython
>>> dis.dis('a + b')
  0     0 RESUME        0
  1     2 LOAD_NAME     0 (a)
        4 LOAD_NAME     1 (b)
        6 BINARY_OP     0 (+)
       10 RETURN_VALUE
Компилятор переводит код в байт-код, а цикл интерпретатора выполняет инструкции по одной: каждая — ветвь кода на C с проверкой типов и счётчиками ссылок.

ВерсияЧто изменилось
3.11Faster CPython: в среднем в 1,25 раза быстрее 3.10
3.13экспериментальный JIT и отдельная сборка без GIL
PyPyдругая реализация языка с JIT-компилятором
Важно
Производительность не была главной целью CPython, а совместимость с C API ограничивает изменения внутреннего устройства.
Изложение. Третья причина — сама эталонная реализация языка. CPython компилирует исходный текст в байт-код, а затем цикл интерпретатора выполняет инструкции по одной: для выражения a + b это загрузка двух имён и инструкция BINARY_OP, внутри которой выясняются типы операндов и вызывается нужная реализация сложения. Каждая инструкция представляет собой ветвь кода на C с проверками, работой со счётчиками ссылок и созданием объектов-результатов. CPython начинали писать задолго до многоядерных процессоров, производительность никогда не была его главной целью, а совместимость с C API, на котором построены NumPy и тысячи расширений, ограничивает изменения внутреннего устройства. Начиная с версии 3.11 проект Faster CPython ускоряет интерпретатор: специализирующий адаптивный интерпретатор подменяет общие инструкции быстрыми вариантами для встреченных типов, и в среднем 3.11 работает в 1,25 раза быстрее 3.10. В 3.13 появились экспериментальный JIT-компилятор и отдельная сборка без глобальной блокировки интерпретатора. Альтернативная реализация PyPy с JIT-компилятором ускоряет чистый Python в несколько раз, но совместима не со всеми расширениями на C.

Версия	Что изменилось
3.11	Faster CPython: в среднем в 1,25 раза быстрее 3.10
3.13	экспериментальный JIT и отдельная сборка без GIL
PyPy	другая реализация языка с JIT-компилятором

Демонстрация:
цена интерпретатора
01
Сумма миллиона чисел циклом for
02
Встроенная sum: тот же счёт на C
03
NumPy: a.sum() над массивом целых
04
Размер одного int и списка из миллиона ссылок
Что наблюдается
Цикл — 12,7 мс, sum — 2,5 мс, NumPy — 0,23 мс: та же сумма в 55 раз быстрее без цикла интерпретатора.
root@lab — запись со стенда
Демонстрация. python3 -m timeit измеряет одну и ту же сумму миллиона целых чисел тремя способами; ключ -s задаёт подготовку, которая в замер не входит. Цикл for с накоплением выполняется за 12,7 мс: на каждой итерации интерпретатор достаёт объект из списка, выбирает операцию сложения по типам и создаёт новый объект-результат. Встроенная функция sum делает то же самое в цикле на C — 2,5 мс, в пять раз быстрее. NumPy хранит числа подряд как 8-байтовые целые и суммирует их без объектов — 0,23 мс, ещё в десять раз быстрее. Последняя команда объясняет разницу: sys.getsizeof(10**6) показывает 28 байт на одно целое число, а список из миллиона элементов занимает 8 000 056 байт только под ссылки на объекты.

ИзмерениеМомент оптимизации, timeit, cProfile, line_profiler, py-spy
02
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. Вторая часть посвящена измерению. Прежде чем ускорять, необходимо выяснить, требуется ли ускорение вообще и какую часть программы ускорять: на эти вопросы отвечают замер времени и профилировщик.

Момент
оптимизации
«О мелкой эффективности следует забыть в 97 % случаев: преждевременная оптимизация — корень всех зол. Но нельзя упускать возможности в тех критических 3 %»
— Д. Кнут, Structured Programming with go to Statements, 1974
Три правила оптимизации
01
Не делать
02
Делать позже
03
Делать оптимально
Пояснение
О производительности думают с самого начала, выбирая алгоритм и структуры данных; переписывают ради скорости — после замера.
Изложение. Фраза Кнута о преждевременной оптимизации является одной из самых цитируемых в программировании, и обычно её приводят без второй половины. Полностью она утверждает, что о мелкой эффективности программисты должны забыть в 97 процентах случаев, но не должны упускать возможности в критических 3 процентах. С ней согласуются три правила оптимизации: не делать; делать позже; делать оптимально. Противоречия здесь нет: думать о производительности необходимо с самого начала, выбирая алгоритм и структуры данных, а переписывать код ради скорости — лишь тогда, когда замер показал, что это требуется. Наименее затратной оказывается оптимизация, которую не пришлось выполнять, поскольку программа и так укладывается в отведённое время. Обратная крайность — архитектура, не продуманная вовсе: медленную систему, построенную без расчёта на объём данных, поздняя оптимизация не исправляет.

Порядок
оптимизации
01
Профилирование
горячее место показывает профилировщик, а не интуиция
02
Сохранение корректности
тесты до оптимизации: быстрый неверный расчёт остаётся неверным
03
Замер после правки
выигрыш подтверждается на тех же данных
S = 1 / ((1 − p) + p / s)

Доля pУскорение sИтог S
1 %10 раз1,01 раза
50 %10 раз1,8 раза
90 %10 раз5,3 раза
90 %бесконечно10 раз
Пояснение
Закон Амдала: ускорение части программы ограничено долей её времени.
Изложение. Главная задача — найти место, к которому имеет смысл прикладывать усилия, а порядок работы задают три правила. Первое — профилирование: угадывать, какая часть программы работает дольше всего, бесполезно, профилировщик отвечает на этот вопрос точно. Количественно это выражает закон Амдала: если доля p времени программы ускорена в s раз, вся программа ускоряется в 1 / ((1 − p) + p / s) раз. Оценим по этой формуле два случая: функция, занимающая 1 процент времени и ускоренная в десять раз, даёт выигрыш в один процент; даже бесконечно быстрая часть, занимающая 90 процентов, ускоряет программу лишь в десять раз. Второе правило — сохранение корректности: оптимизация легко нарушает правильность незаметно, оставляя результат правдоподобным, но неверным, поэтому перед переписыванием фрагмент покрывают тестами. Третье правило — замер после правки на тех же данных: выигрыш, не подтверждённый замером, считать выигрышем нельзя.

Доля p	Ускорение s	Итог S
1 %	10 раз	1,01 раза
50 %	10 раз	1,8 раза
90 %	10 раз	5,3 раза
90 %	бесконечно	10 раз

timeit:
время фрагмента
$ python3 -m timeit -s 'setup' 'stmt'
500000 loops, best of 5: 527 nsec per loop
 
>>> import timeit
>>> timeit.timeit('f(s)', setup, number=10000)
0.0052                    # сумма запусков
>>> timeit.repeat('f(s)', setup, number=10000)
[0.0053, 0.0052, ...]     # пять прогонов
01
number и repeat
число запусков в прогоне и число прогонов; берётся минимум
02
Условия замера
сборщик мусора отключён, часы — time.perf_counter
03
IPython
%timeit печатает среднее и отклонение по 7 прогонам
Важно
timeit.timeit возвращает сумму всех запусков, а не время одного.
Изложение. Когда требуется измерить время одного фрагмента, а не снимать полный профиль, применяется модуль timeit из стандартной библиотеки. Из командной строки он запускается как python3 -m timeit: ключ -s задаёт подготовку, которая в замер не входит, а сам фрагмент передаётся строкой; модуль подбирает число запусков так, чтобы прогон длился не меньше 0,2 секунды, повторяет прогон пять раз и печатает лучшее время одного запуска. Из программы вызываются функции timeit.timeit и timeit.repeat: первая возвращает суммарное время number запусков, вторая — список таких сумм по нескольким прогонам. Берут минимум: более медленные прогоны замедлены не самим кодом, а другими процессами. На время замера timeit отключает сборщик мусора, а время измеряет по time.perf_counter; если код создаёт много объектов и сборка мусора входит в его цену, сборщик включают первой строкой подготовки: gc.enable(). Код передаётся строками: timeit подставляет его в шаблон, что позволяет исключить накладные расходы на вызов функции-обёртки. В IPython для той же цели предусмотрена команда %timeit, печатающая среднее время и стандартное отклонение по семи прогонам.

Демонстрация:
замер timeit
01
Две функции разворота строки в rev.py
02
Цикл со сложением строк: 527 нс
03
Срез s[::-1]: 42 нс
04
Ключи -n и -r: число запусков и прогонов
Что наблюдается
Срез быстрее цикла в 12 раз; timeit сам подбирает число запусков и печатает лучшее из пяти прогонов.
root@lab — запись со стенда
Демонстрация. Файл rev.py содержит две функции разворота строки из 26 букв: reverse_loop собирает результат сложением строк в цикле, reverse_slice возвращает срез s[::-1]. python3 -m timeit с подготовкой from rev import … измеряет каждую: цикл выполняется за 527 нс, срез — за 42 нс, в 12 раз быстрее. Число запусков timeit подобрал сам — 500 000 и 5 000 000, — чтобы каждый прогон длился достаточно долго, и вывел лучшее из пяти прогонов. Ключи -n 1000 -r 3 задают число запусков и прогонов вручную; результат почти не меняется, потому что срез стабилен по времени.

Профилировщики

ИнструментЧто показываетКак запустить
cProfile + pstatsвремя и число вызовов функцийpython3 -m cProfile -o prof.out f.py
SnakeVizпрофиль cProfile диаграммойsnakeviz prof.out
line_profilerвремя каждой строки функции@profile и kernprof -l -v f.py
py-spyвыборки стека работающей программыpy-spy record --pid N
tracemallocвыделения памяти по строкамtracemalloc.start()
Пояснение
cProfile перехватывает каждый вызов и замедляет программу; py-spy опрашивает стек снаружи и почти не влияет на неё.
Для физика
Расчёт, идущий третьи сутки, профилируют py-spy: без перезапуска и правки кода.
Изложение. Профилировщик показывает, на что программа тратит время. Базовый набор входит в стандартную поставку: cProfile собирает профиль по функциям, pstats разбирает и сортирует собранное; сторонний пакет SnakeViz показывает тот же профиль диаграммой в браузере. cProfile представляет собой детерминированный профилировщик: он перехватывает каждый вызов и возврат, что обеспечивает точный подсчёт вызовов, но заметно замедляет программу. line_profiler измеряет время каждой строки отмеченной декоратором @profile функции и запускается командой kernprof. py-spy — сэмплирующий профилировщик: он по умолчанию сто раз в секунду считывает стек работающего процесса снаружи, из его памяти, и почти не влияет на скорость. tracemalloc из стандартной библиотеки отслеживает выделения памяти и рассматривается в части о памяти.
Для физиков: py-spy подключается к уже запущенному расчёту по номеру процесса, что позволяет исследовать расчёт, идущий третьи сутки, без перезапуска и правки кода; в Linux подключение к процессу, запущенному не из py-spy, требует прав root (sudo) или значения kernel.yama.ptrace_scope = 0.

Инструмент	Что показывает	Как запустить
cProfile + pstats	время и число вызовов функций	python3 -m cProfile -o prof.out f.py
SnakeViz	профиль cProfile диаграммой	snakeviz prof.out
line_profiler	время каждой строки функции	@profile и kernprof -l -v f.py
py-spy	выборки стека работающей программы	py-spy record --pid N
tracemalloc	выделения памяти по строкам	tracemalloc.start()

Демонстрация:
профиль cProfile
01
matrix.py: 16 умножений матриц 64×64
02
cProfile с записью профиля в prof.out
03
pstats: пять функций с наибольшим временем
04
matrix_product — 0,295 с из 0,308 с
Что наблюдается
96 % времени уходит в matrix_product; генерация случайных матриц занимает проценты.
root@lab — запись со стенда
Демонстрация. Сквозным примером лекции является умножение матриц: класс Matrix хранит матрицу списком списков, а matrix_product перемножает матрицы тройным циклом по индексам. Функция bench создаёт две случайные матрицы 64×64 и перемножает их 16 раз. python3 -m cProfile -o prof.out matrix.py выполняет программу и сохраняет профиль в файл, а pstats.Stats сортирует его по собственному времени функций (tottime) и печатает пять первых строк. Всего выполнено 66 666 вызовов за 0,308 с, из них 0,295 с — собственное время matrix_product, 96 процентов. Функции модуля random, создающие матрицы, вызывались по 8192 раза, но заняли миллисекунды. Таким образом, ускорять имеет смысл только умножение: ускорение генерации матриц, сколь угодно большое, по закону Амдала даст не больше четырёх процентов.

Демонстрация:
построчный профиль
01
@profile над matrix_product
02
kernprof -l -v: время каждой строки
03
Умножение со сложением — 60,7 %
04
Внутренний цикл for — 38,7 %
Что наблюдается
Почти всё время — две строки внутреннего цикла, выполненные 4 млн раз; профилирование замедлило счёт вчетверо.
root@lab — запись со стенда
Демонстрация. Декоратор @profile отмечает функцию matrix_product, а kernprof -l -v запускает программу под line_profiler и печатает время каждой строки отмеченной функции; ключ -u 1e-3 выводит время в миллисекундах. Тело внутреннего цикла Z[i][k] += X[i][j] * Y[j][k] выполнено 4 194 304 раза и заняло 60,7 процента времени функции, сама строка внутреннего цикла for — 38,7 процента; остальные строки — доли процента. Профиль позволяет установить, что именно ускорять: обращения по индексам во внутреннем цикле. Общее время под line_profiler — 1,28 с вместо 0,3 с без него: построчное профилирование замедляет программу в несколько раз, поэтому абсолютные времена под профилировщиком не сравнивают с обычными замерами, сравнивают только доли.

Демонстрация:
py-spy без перезапуска
01
beam.py: долгий расчёт в фоне
02
py-spy dump: текущий стек процесса
03
py-spy record: три секунды выборок
04
Остановка расчёта
Что наблюдается
Без перезапуска и правки кода видно, что процесс занят space_charge, строка 11; за 3 с снято 299 выборок.
root@lab — запись со стенда
Демонстрация. Программа beam.py моделирует долгий расчёт: на каждом шаге она вычисляет силы пространственного заряда для 400 частиц, а тысяча шагов занимает десятки секунд — дольше, чем длится демонстрация. Она запускается в фоне, а py-spy dump --pid $! подключается к процессу по номеру и печатает текущий стек: функция space_charge, строка 11, вызванная из simulate. py-spy record в течение трёх секунд снимает 100 выборок стека в секунду и записывает в файл beam.svg flame graph, представляющий собой диаграмму, в которой ширина каждой функции пропорциональна доле выборок, пришедшихся на неё. Наблюдение в реальном времени обеспечивает py-spy top — непрерывно обновляемая сводка, похожая на top. После показа процесс останавливается командой kill.

Приёмы языкаВстроенные функции, структуры данных, идиомы и проверка замером
03
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. Третья часть объединяет приёмы, не требующие ничего, кроме самого языка. Выигрыш обычно скромнее, чем у компиляции, но и затраты отсутствуют: код остаётся обычным Python, понятным любому коллеге. Все числа в этой части получены на стенде с Python 3.12.

Уровни
оптимизации

УровеньЧто меняетсяПотенциал
архитектуракакие данные, в каком объёме и где хранятсянаибольший
алгоритмы и структуры данныхO(n²) → O(n log n), список → множествобольшой
реализациязапись выбранного алгоритма на языкепроценты и единицы раз
компиляцияJIT и AOT: Numba, Cythonдесятки и сотни раз
исполнениекеши, ленивые вычислениязависит от данных
Важно
Оптимизация имеет цену: время без гарантии результата, усложнение кода, а иногда и память.
Пояснение
Наибольший выигрыш дают архитектура и алгоритмы; ошибка на этих уровнях, найденная поздно, обходится дороже всего.
Изложение. Оптимизация не сводится к правке кода: ускорять можно на нескольких уровнях. Архитектура определяет, какие данные система обрабатывает, каким способом и где их хранит. Алгоритм и структура данных определяют асимптотику: замена квадратичного алгоритма на n log n или списка на множество при проверке вхождения даёт выигрыш, недостижимый никакой компиляцией. Реализация представляет собой способ записи выбранного алгоритма на языке; ей посвящена эта часть. Ниже лежат оптимизации во время компиляции, которые выполняют JIT- и AOT-компиляторы, и во время исполнения: кеши, ленивые вычисления, специализация под найденный профилировщиком горячий путь. За оптимизацию всегда приходится платить: она отнимает время без гарантии результата, код становится сложнее, а выигрыш в скорости легко оборачивается проигрышем в памяти.

Уровень	Что меняется	Потенциал
архитектура	какие данные, в каком объёме и где хранятся	наибольший
алгоритмы и структуры данных	O(n²) → O(n log n), список → множество	большой
реализация	запись выбранного алгоритма на языке	проценты и единицы раз
компиляция	JIT и AOT: Numba, Cython	десятки и сотни раз
исполнение	кеши, ленивые вычисления	зависит от данных

Встроенные
функции
def count(elements):
    total = 0
    for item in elements:
        total += 1
    return total
 
count(xs)              # 9,6 мс
len(xs)                # 12 нс
[x for x in xs if x % 2]           # 10,4 мс
list(filter(lambda x: x % 2, xs))  # 22,7 мс
 
[s for s in strs if s.isdigit()]   # 13,3 мс
list(filter(str.isdigit, strs))    # 12,4 мс
Пояснение
filter с lambda вызывает функцию Python на каждом элементе и проигрывает включению вдвое; с методом на C filter выигрывает.
Изложение. Если требуемая операция уже реализована встроенной функцией, почти всегда быстрее воспользоваться ею: встроенные функции написаны на C и обходятся без цикла на уровне интерпретатора. Подсчёт элементов списка из миллиона чисел циклом занимает 9,6 мс, а len возвращает длину, хранящуюся в самом списке, за 12 нс. Однако встроенные функции не ускоряют код сами по себе. Отбор нечётных чисел функцией filter с lambda вдвое медленнее спискового включения: filter вызывает функцию Python на каждом элементе, а включение выполняет проверку прямо в цикле байт-кода. Достаточно заменить lambda методом str.isdigit, написанным на C, и filter становится быстрее включения. Поэтому каждый конкретный случай необходимо проверять замером.

Структуры
данных


Операция, миллион элементовСложностьВремя на стенде
x in списокO(n)2,0 мс
x in множествоO(1)11 нс
set(список)O(n)10,9 мс
sorted(data) и data.sort()O(n log n)126 мс

Для физика
Множество выигрывает после пяти-шести поисков: построение длится как пять проверок по списку.
Пояснение
На Python 3.12 сортировка на месте и sorted со своей копией не различаются по времени: копия занимает 1,9 мс из 126.
Важно
Множество занимает больше памяти, чем список: хеш-таблица резервирует запас.
Изложение. Самый дешёвый выигрыш обычно даёт не микрооптимизация, а смена структуры данных. Проверка вхождения в список просматривает элементы по порядку, пока не найдёт нужный, — O(n): поиск числа в середине списка из миллиона элементов занимает 2 мс. Множество находит элемент по хешу за O(1) — 11 нс, в 180 тысяч раз быстрее. Платой за это являются время на построение множества, 10,9 мс, и память: хеш-таблица резервирует запас пустых ячеек. Строить множество ради одной проверки бессмысленно, а ради тысячи проверок необходимо; сравнив время построения со временем одной проверки по списку, получим, что затраты возмещаются после пяти-шести поисков. Сортировка на месте методом list.sort избегает копии, но на стенде с Python 3.12 разница с sorted незаметна: копия миллиона ссылок занимает 1,9 мс из 126. Подробно сложность операций с коллекциями разобрана в главе книги «Сложность операций с коллекциями».

Операция, миллион элементов	Сложность	Время на стенде
x in список	O(n)	2,0 мс
x in множество	O(1)	11 нс
set(список)	O(n)	10,9 мс
sorted(data) и data.sort()	O(n log n)	126 мс

Идиомы
и их цена


БыстрееМедленнее100 000 повторений
if not xif x == False0,77 и 1,12 мс
if not xsif len(xs) == 01,02 и 1,77 мс
[]list()1,44 и 2,20 мс
try при редком исключенииhasattr1,47 и 2,55 мс
hasattr при частом отсутствииtry1,84 и 15,9 мс
локальная переменнаяглобальная1,46 и 2,08 мс

>>> dis.dis('[]')
  1     2 BUILD_LIST     0
>>> dis.dis('list()')
  1     2 PUSH_NULL
        4 LOAD_NAME      0 (list)
        6 CALL           0
Пояснение
Литерал — одна инструкция, а list() требует поиска имени и вызова.
Изложение. Идиоматичная запись часто оказывается и самой быстрой. Проверка if not x быстрее сравнения с False, а if not xs — вызова len и сравнения с нулём. Литерал [] быстрее вызова list(): байт-код показывает, что литерал компилируется в одну инструкцию BUILD_LIST, тогда как вызов требует найти имя list, которое могло быть переопределено, и вызвать его. Выбор между проверкой hasattr и обработкой исключения зависит от частоты: если атрибут обычно есть, try дешевле, а если исключение возникает на каждой итерации, try оказывается почти в девять раз медленнее. Локальная переменная быстрее глобальной: локальные имена хранятся в массиве кадра и адресуются по номеру, глобальные ищутся в словаре модуля. Вызов функции также требует времени: на стенде квадраты десяти тысяч чисел через функцию square вычисляются за 304 мкс, а выражением прямо во включении — за 217 мкс.

Быстрее	Медленнее	100 000 повторений
if not x	if x == False	0,77 и 1,12 мс
if not xs	if len(xs) == 0	1,02 и 1,77 мс
[]	list()	1,44 и 2,20 мс
try при редком исключении	hasattr	1,47 и 2,55 мс
hasattr при частом отсутствии	try	1,84 и 15,9 мс
локальная переменная	глобальная	1,46 и 2,08 мс

Проверка
замером
def plain(xs):
    result = []
    for x in xs:
        result.append(x)     # 6,1 мс
 
def cached(xs):
    result = []
    append = result.append
    for x in xs:
        append(x)            # 7,7 мс
Совет сохранять метод в локальной переменной годами переходит из подборки в подборку, но с Python 3.11 он проигрывает: интерпретатор специализирует вызов метода в байт-коде.
Важно
Множественное присваивание и подобные приёмы дают проценты ценой читаемости; их применяют, только если профилировщик показал необходимость.
Изложение. Рекомендации по оптимизации устаревают с выходом новых версий интерпретатора. Поиск атрибута в Python имеет свою цену, и естественным решением кажется найти метод append один раз и сохранить его в локальную переменную. Совет годами переходил из одной подборки в другую, однако на Python 3.12 сохранённая ссылка оказывается на четверть медленнее обычного вызова result.append(x): начиная с версии 3.11 интерпретатор специализирует вызов метода списка в байт-коде, и прямой вызов обходится дешевле вызова через локальную переменную. Таким образом, каждую рекомендацию проверяют замером на используемой версии интерпретатора. Отдельную группу составляют приёмы, ухудшающие читаемость ради нескольких процентов: например, присваивание десяти переменных одной строкой через распаковку кортежа на стенде быстрее десяти отдельных строк на 19 процентов, но читается хуже. Такие приёмы применяют, только если профилировщик показал, что эти проценты необходимы.

Демонстрация:
приёмы и замер
01
Проверка вхождения: список и множество
02
Стоимость построения множества
03
Глобальная и локальная переменная
04
result.append и сохранённый метод
Что наблюдается
Множество в 180 000 раз быстрее списка, но строится 11 мс; локальная переменная быстрее на 30 %, а сохранённый метод append замедляет заполнение на 26 %.
root@lab — запись со стенда
Демонстрация. Скрипт tips.py измеряет приёмы из главы книги функцией timeit.repeat: число запусков подбирает Timer.autorange, из пяти прогонов берётся лучший, время печатается в микросекундах на один запуск. Поиск числа 500 000 в списке из миллиона элементов занимает около 2000 мкс, в множестве — 0,011 мкс, но построение множества обходится в 10 900 мкс. Сто тысяч увеличений глобальной переменной занимают 2,08 мс, локальной — 1,46 мс, на 30 процентов меньше. Последние две строки опровергают популярный совет: заполнение списка из миллиона элементов через сохранённую ссылку append = result.append на 26 процентов медленнее обычного result.append(x).

ПамятьИзмерение памяти, компактные структуры, __slots__, утечки
04
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. Четвёртая часть посвящена памяти. Расчёт, не помещающийся в оперативную память, не выполняется вовсе или многократно замедляется: выделение памяти завершается ошибкой MemoryError, а при подкачке данные читаются с диска. Кроме того, компактные данные быстрее обрабатываются: процессор читает их из кеша, а не из памяти.

Измерение
памяти
import sys, tracemalloc
 
xs = [float(i) for i in range(10**6)]
sys.getsizeof(xs)          # 8 МиБ: ссылки
 
tracemalloc.start()
xs = [float(i) for i in range(10**6)]
tracemalloc.get_traced_memory()
# (текущая, пиковая) — 31 МиБ
tracemalloc.stop()

СредствоЧто измеряет
sys.getsizeof(x)размер самого объекта без вложенных
tracemallocвыделения памяти интерпретатором, по строкам кода
массив NumPy: nbytesточный размер данных массива
RSS процесса: ps, topвся память процесса вместе с запасом аллокатора
Важно
getsizeof не учитывает объекты, на которые ссылается контейнер: по его оценке список частиц равен по размеру списку чисел.
Изложение. Измерение памяти в Python затруднено, и первый способ, подсказанный документацией, вводит в заблуждение. Функция sys.getsizeof возвращает размер самого объекта: для списка это массив ссылок, но не объекты, на которые эти ссылки ведут. Поэтому по оценке getsizeof список из миллиона экземпляров класса занимает столько же, сколько список из миллиона чисел. Модуль tracemalloc из стандартной библиотеки отслеживает выделения памяти интерпретатором: get_traced_memory возвращает текущий и пиковый объём, а снимки take_snapshot позволяют установить, какие строки кода выделили больше всего памяти. У массива NumPy точный размер данных возвращает атрибут nbytes. Память процесса в целом, RSS, показывают ps и top, но она включает запас, который аллокатор удерживает после освобождения объектов, поэтому для сравнения структур данных она плохо подходит.

Средство	Что измеряет
sys.getsizeof(x)	размер самого объекта без вложенных
tracemalloc	выделения памяти интерпретатором, по строкам кода
массив NumPy: nbytes	точный размер данных массива
RSS процесса: ps, top	вся память процесса вместе с запасом аллокатора

Компактные
структуры

СтруктураНа элементПочему
list из float32 байтассылка 8 байт и объект 24 байта
array('d')8 байтчисла подряд, без объектов
np.ndarray float648 байтчисла подряд и операции на C
tuple вместо listбез запасане растёт, не резервирует место под append
__slots__−32 байта на объектнет словаря атрибутов
range, генератор0 байтзначение вычисляется по запросу
class Particle:
    __slots__ = ('x', 'y', 'z')
    def __init__(self, x, y, z):
        self.x, self.y, self.z = x, y, z
Пояснение
__slots__ указывают в каждом классе иерархии: иначе словарь атрибутов возвращается.
Изложение. Список чисел хранит не числа, а ссылки на объекты: на стенде каждый элемент списка вещественных чисел обходится в 8 байт ссылки и 24 байта объекта float. Модуль array хранит числа примитивного типа подряд, по 8 байт на double, а массив NumPy устроен так же, но вдобавок поддерживает вычисления. Кортеж немного компактнее списка: ему не требуется запас под рост, который резервирует список, созданный включением или append. Атрибут __slots__ отменяет у экземпляров словарь атрибутов и размещает объявленные поля в фиксированных ячейках; на стенде это экономит четверть памяти на частицу из трёх полей. С наследованием __slots__ неудобен: его указывают в каждом классе иерархии, иначе словарь атрибутов возвращается и экономия исчезает. Иногда последовательность не требуется хранить вовсе: range вычисляет элемент по индексу, а генератор выдаёт значения по одному.

Структура	На элемент	Почему
list из float	32 байта	ссылка 8 байт и объект 24 байта
array('d')	8 байт	числа подряд, без объектов
np.ndarray float64	8 байт	числа подряд и операции на C
tuple вместо list	без запаса	не растёт, не резервирует место под append
__slots__	−32 байта на объект	нет словаря атрибутов
range, генератор	0 байт	значение вычисляется по запросу

Демонстрация:
память под данные
01
tracemalloc: память под структуру целиком
02
Список из миллиона float — 30,9 МиБ
03
array и NumPy — 7,8 и 7,6 МиБ
04
Сто тысяч частиц: со словарём и со __slots__
Что наблюдается
Список чисел вчетверо объёмнее массива, а getsizeof учитывает лишь 8 МиБ ссылок; __slots__ экономит четверть памяти.
root@lab — запись со стенда
Демонстрация. Скрипт mem.py создаёт структуры данных под включённым tracemalloc и печатает объём выделенной памяти рядом со значением sys.getsizeof. Список из миллиона вещественных чисел занимает 30,9 МиБ, а getsizeof показывает 8,1 МиБ — только массив ссылок. array('d') и массив NumPy хранят те же числа в 7,8 и 7,6 МиБ, вчетверо компактнее, и здесь getsizeof с tracemalloc совпадают. Сто тысяч частиц, у которых все три координаты ссылаются на одно целое число, занимают 12,2 МиБ с обычными экземплярами и 9,1 МиБ с __slots__, на 32 байта меньше на экземпляр; getsizeof списка в обоих случаях — 0,8 МиБ ссылок.

Утечки
памяти
def append_to(x, acc=[]):     # список один
    acc.append(x)            # на все вызовы
    return acc
 
class Model:
    cache = {}               # не очищается
    def calc(self, arg):
        if arg not in self.cache:
            self.cache[arg] = solve(arg)
        return self.cache[arg]
Объект без ссылок освобождается сразу, недостижимые циклы находит сборщик мусора. Утечкой в Python называют долгоживущие ненужные объекты: на них ссылаются кеш, глобальный список, аргумент по умолчанию.
Пояснение
Кеш ограничивают: functools.lru_cache(maxsize=…). Рост памяти показывает сравнение двух снимков tracemalloc.
Изложение. В смысле C++ утечек в Python почти нет: объект освобождается, когда на него не остаётся ссылок, а циклические ссылки находит поколенческий сборщик мусора, разобранный в главе книги «Объекты и память». Утечкой в Python обычно называют долгоживущие ненужные объекты, на которые по-прежнему существуют ссылки. Классических источников три. Изменяемый аргумент по умолчанию создаётся один раз при определении функции и накапливает данные всех вызовов. Кеш, заведённый на атрибуте класса, является общим для всех экземпляров и растёт неограниченно: записи из него не удаляются. Большой массив в локальной переменной длинной функции живёт до её конца, даже если больше не нужен. Кеш ограничивают по размеру, например декоратором functools.lru_cache с параметром maxsize; применённый к методу, такой кеш хранит ссылку на self и удерживает экземпляры, пока запись не вытеснена. Место роста памяти находят сравнением двух снимков tracemalloc, сделанных до и после подозрительного участка.

КомпиляцияNumba и Cython на примерах умножения матриц и пространственного заряда
05
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. В пятой части рассматривается компиляция горячего цикла в машинный код. Приёмы языка ускоряют умножение матриц лишь в 2,8 раза; дальнейшему ускорению препятствует сам интерпретатор, и выигрыш в десятки и сотни раз обеспечивает только компиляция — при первом вызове, через Numba, или заранее, через Cython.

Пример:
умножение матриц
def matrix_product(X, Y):
    n_xrows, n_xcols = X.shape
    n_yrows, n_ycols = Y.shape
    Z = Matrix.zeros((n_xrows, n_ycols))
    for i in range(n_xrows):
        for j in range(n_xcols):
            for k in range(n_ycols):
                Z[i][k] += X[i][j] * Y[j][k]
    return Z

Чистый Python, 64×64Время
тройной цикл по индексам17,7 мс
перестановка циклов, накопитель11,3 мс
транспонирование и zip8,9 мс
sum(map(operator.mul, row, col))6,3 мс
Пояснение
Приёмы языка ускоряют в 2,8 раза; дальше каждое умножение по-прежнему проходит через интерпретатор.
Изложение. Сквозной пример — умножение матриц, операция, встречающаяся в любом физическом расчёте. Класс Matrix хранит матрицу списком списков, а matrix_product перемножает матрицы тройным циклом по индексам: на стенде умножение двух матриц 64×64 занимает 17,7 мс. Профиль показал, что время уходит на обращения по индексам во внутреннем цикле. Приёмами языка их можно сократить: перестановка циклов с накоплением суммы в локальной переменной даёт 11,3 мс, транспонирование второй матрицы и проход по парам элементов через zip — 8,9 мс, а сумма sum(map(operator.mul, row, col)) переносит внутренний цикл в C — 6,3 мс. В итоге ускорение составляет 2,8 раза, однако каждое умножение по-прежнему выполняется интерпретатором над объектами, и большего выигрыша этим путём не достичь.

Чистый Python, 64×64	Время
тройной цикл по индексам	17,7 мс
перестановка циклов, накопитель	11,3 мс
транспонирование и zip	8,9 мс
sum(map(operator.mul, row, col))	6,3 мс

Numba:
JIT-компиляция
import numba
import numpy as np
 
@numba.njit
def jit_product(X, Y):
    n, m = X.shape
    p = Y.shape[1]
    Z = np.zeros((n, p))
    for i in range(n):
        for k in range(p):
            acc = 0.0
            for j in range(m):
                acc += X[i, j] * Y[j, k]
            Z[i, k] = acc
    return Z
01
При первом вызове
LLVM компилирует функцию в машинный код под типы аргументов
02
Числа и массивы
работает с NumPy, а не со списками объектов
03
Прогрев
замер — после первого вызова: он включает компиляцию
04
Параметры
cache=True — код на диск, parallel=True и prange — по ядрам
Важно
Компиляция не меняет асимптотику: тройной цикл остаётся O(n³).
Изложение. Numba представляет собой JIT-компилятор для числового кода на Python. Декоратор @numba.njit заменяет функцию обёрткой, которая при первом вызове выводит типы аргументов, переводит функцию в промежуточное представление LLVM и компилирует его в машинный код; при следующих вызовах с теми же типами сразу выполняется готовый код. Numba работает с числами и массивами NumPy, а не со списками произвольных объектов, поэтому матрицы передаются как массивы float64. Лучше всего Numba ускоряет явные циклы по элементам: выражение над срезами внутри цикла, например (X[i, :] * Y[:, k]).sum(), создаёт временные массивы и замедляет функцию в несколько раз. Первый вызов включает компиляцию и длится в тысячи раз дольше самого счёта, поэтому замер проводят после прогрева; параметр cache=True позволяет сохранить скомпилированный код на диске между запусками, а parallel=True вместе с prange — распределить итерации цикла по ядрам. Компиляция ускоряет константу, но не меняет асимптотику: квадратичный алгоритм остаётся квадратичным.

Демонстрация:
Numba
01
nb.py: умножение 64×64 на списках и массивах
02
Чистый Python — 17,8 мс
03
Первый вызов @njit — 156 мс: компиляция
04
Повторный вызов — 0,08 мс
05
X @ Y в NumPy — 0,06 мс
Что наблюдается
После прогрева Numba быстрее чистого Python в 217 раз; первый вызов тратит 156 мс на компиляцию.
root@lab — запись со стенда
Демонстрация. Скрипт nb.py перемножает матрицы 64×64 тремя способами и печатает время одного вызова по time.perf_counter. Чистый Python на списках — 17,8 мс. Первый вызов функции с декоратором @numba.njit занимает 156 мс: в это время входит компиляция под типы аргументов. Второй вызов — 0,082 мс, в 217 раз быстрее чистого Python. Умножение средствами NumPy, X @ Y, — 0,057 мс за единственный вызов; повторные замеры timeit в демонстрациях «Cython» и «BLAS и тип данных» дают около 0,02 мс, вчетверо меньше, чем у скомпилированного тройного цикла. Последняя строка проверяет, что результаты Numba и NumPy совпадают.

Пространственный
заряд
F[i] = Σ (r[i] − r[j]) / |r[i] − r[j]|³,  j ≠ i
from numba import njit, prange
 
def forces_par(x, y, z, fx):
    for i in prange(len(x)):       # по ядрам
        for j in range(len(x)):
            if i != j:
                r3 = ((x[j]-x[i])**2 + (y[j]-y[i])**2
                      + (z[j]-z[i])**2)**1.5
                fx[i] += (x[i]-x[j]) / r3
 
par = njit(parallel=True)(forces_par)
Для физика
Сильноточный пучок расталкивает сам себя; в коде REDPIC суммирование кулоновских сил между макрочастицами — самая затратная функция расчёта.
Пояснение
Каждая частица взаимодействует с каждой: сложность O(N²), и слагаемые для разных i независимы — внешний цикл распараллеливается.
Изложение. Умножение матриц — учебный пример; перейдём к реальной задаче, на которой выигрыш оказывается таким же. В коде REDPIC, которому в книге посвящена отдельная глава, пучок представлен ансамблем макрочастиц, а самой затратной функцией расчёта является суммирование кулоновских сил между частицами — эффект пространственного заряда, из-за которого сильноточный пучок расталкивает сам себя. Каждая частица взаимодействует с каждой, поэтому сложность — O(N²). Цикл, записанный напрямую, допускает и JIT-компиляцию, и распараллеливание: слагаемые для разных i вычисляются независимо. Для параллельной версии внешний range заменяется на prange, что указывает Numba, какой цикл распределить по ядрам. Запись njit(parallel=True)(forces_par) вместо декоратора удобна, когда одну функцию измеряют в нескольких режимах.

Демонстрация:
пространственный заряд
01
nproc: два ядра виртуальной машины
02
Цикл Python над списками и над массивами
03
@njit: в 26 раз быстрее цикла над списками
04
parallel=True и prange: ещё вдвое
05
NumPy с массивами N×N
Что наблюдается
При удвоении N время растёт вчетверо — O(N²); @njit быстрее цикла над списками в 26 раз, над массивами — в 120; prange — вдвое на двух ядрах.
root@lab — запись со стенда
Демонстрация. Скрипт space.py вычисляет силы пространственного заряда для N случайных частиц пятью способами; функции Numba компилируются заранее, вызовом на двух частицах, чтобы компиляция не попала в замер. Один и тот же цикл Python выполняется над списками и над массивами NumPy, и только до N = 1024: над списками при 256, 512 и 1024 частицах — 8,8, 36,5 и 154 мс, над массивами — в 4,5–5 раз дольше, потому что каждое обращение x[j] к массиву создаёт объект numpy.float64, арифметика над которым дороже, чем над float. @njit считает те же силы за 0,34, 1,45 и 5,7 мс — в 26 раз быстрее цикла над списками и в 120 раз быстрее цикла над массивами. Отметим, что при каждом удвоении N время растёт вчетверо: компиляция ускоряет константу, но сложность остаётся O(N²). Параллельная версия с prange на двух ядрах виртуальной машины быстрее последовательной вдвое: 12,0 мс против 23,6 при N = 2048. NumPy-версия без циклов строит массивы разностей N×N и при N = 2048 считает 56 мс: временные массивы размером N² многократно записываются в память и читаются из неё, и векторизация уступает скомпилированному циклу. Для каждого N скрипт сверяет результаты всех версий и в конце сообщает, что они совпадают.

Cython:
Python с типами C
@cython.boundscheck(False)
@cython.wraparound(False)
def unchecked(double[:, :] X, double[:, :] Y):
    cdef Py_ssize_t n = X.shape[0], m = X.shape[1]
    cdef Py_ssize_t p = Y.shape[1]
    cdef Py_ssize_t i, j, k
    cdef double acc
    Z = np.zeros((n, p))
    cdef double[:, :] Zv = Z
    for i in range(n):
        for k in range(p):
            acc = 0.0
            for j in range(m):
                acc += X[i, j] * Y[j, k]
            Zv[i, k] = acc
    return Z
01
Сборка
.pyx → cythonize -i → код на C → компилятор → модуль .so
02
Типы
cdef и представления памяти double[:, :] убирают вызовы Python
03
Без проверок
boundscheck и wraparound отключают проверку индексов
04
cython -a
HTML-отчёт: жёлтые строки обращаются к C API Python
Важно
Без проверки границ ошибка индекса не вызывает IndexError, а портит чужую память.
Изложение. Второй путь к машинному коду — Cython, язык, представляющий собой Python с объявлениями типов C. Файл .pyx транслируется в код на C, компилятор C собирает из него модуль расширения, и модуль импортируется как обычный; команда cythonize -i выполняет всю сборку на месте. Без объявлений типов Cython порождает вызовы C API Python для каждой операции, и выигрыша почти нет. Объявления cdef для локальных переменных и типизированные представления памяти double[:, :] для массивов позволяют обращаться к числам напрямую, без объектов. Декораторы boundscheck(False) и wraparound(False) отключают проверку выхода за границы и поддержку отрицательных индексов. Взамен ошибка в индексе перестаёт вызывать IndexError и приводит к обращению в чужую память без какого-либо сообщения, поэтому границы такого цикла выверяют вручную. Команда cython -a строит HTML-отчёт, в котором жёлтым отмечены строки, всё ещё обращающиеся к C API Python.

Демонстрация:
Cython
01
cythonize -i: mm.pyx → mm.c → модуль .so
02
untyped: тот же код без типов — 67 мс
03
typed: типы C и представления памяти — 0,11 мс
04
unchecked: без проверок границ — 0,095 мс
05
X @ Y — 0,021 мс
Что наблюдается
Без типов Cython медленнее чистого Python; объявленные типы ускоряют его в 620 раз; NumPy остаётся быстрее.
root@lab — запись со стенда
Демонстрация. Файл mm.pyx содержит три версии умножения матриц: untyped — тройной цикл без объявлений типов над массивами NumPy, typed — те же циклы с cdef-переменными и представлениями памяти double[:, :], unchecked — то же с отключёнными boundscheck и wraparound. cythonize -i транслирует файл в mm.c и собирает модуль mm.cpython-312-aarch64-linux-gnu.so компилятором gcc. Скрипт cy.py проверяет совпадение результатов с X @ Y и измеряет время: версия без типов работает 67 мс — медленнее чистого Python на списках, потому что каждое обращение к элементу массива создаёт объект числа. Объявленные типы дают 0,108 мс, в 620 раз быстрее версии без типов и в 165 раз быстрее чистого Python на списках; отключение проверок границ сокращает время ещё на 12 процентов, до 0,095 мс. Умножение средствами NumPy — 0,021 мс.

Numba
или Cython

NumbaCython
компиляцияпри первом вызовепри сборке
что требуетсядекоратор.pyx, компилятор C
типывыводятся по аргументамобъявляются вручную
данныечисла и массивы NumPyлюбые, включая структуры C
сильная стороначисленные циклы без правки кодаобёртки библиотек C и C++
Для физика
Для числового цикла над массивами в научном коде первый выбор — Numba: один декоратор, без сборки.
Пояснение
Cython выбирают, когда нужен готовый модуль без компилятора на стороне пользователя или обёртка библиотеки на C.
Изложение. Numba компилирует функцию во время выполнения, при первом вызове, и выводит типы по аргументам: для численного цикла над массивами достаточно добавить декоратор, и никакой сборки не требуется. Её ограничение — поддерживаемое подмножество Python: числа, массивы NumPy, кортежи, простые классы через jitclass; со строками, словарями произвольных объектов и сторонними библиотеками она работает плохо или не работает вовсе. Cython компилирует заранее: типы объявляются вручную, требуется компилятор C и шаг сборки, однако результат представляет собой обычный модуль расширения, который распространяется готовым пакетом формата wheel и не требует компилятора на стороне пользователя, а из кода на Cython можно напрямую вызывать функции библиотек на C и C++. Таким образом, Cython чаще используют для обёрток существующих библиотек, а для ускорения числовых циклов научного кода первым выбором служит Numba.

	Numba	Cython
компиляция	при первом вызове	при сборке
что требуется	декоратор	.pyx, компилятор C
типы	выводятся по аргументам	объявляются вручную
данные	числа и массивы NumPy	любые, включая структуры C
сильная сторона	численные циклы без правки кода	обёртки библиотек C и C++

NumPy и BLASВекторизация: цикл в коде на C, библиотека BLAS, тип данных и плата памятью
06
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. Шестая часть возвращается к векторизации: прежде чем компилировать Python, имеет смысл попытаться не писать циклы вообще. Рассмотрим, за счёт чего NumPy выполняет операцию над массивом быстрее цикла, какую роль играет библиотека BLAS и почему тип данных определяет скорость не меньше, чем запись выражения.

Векторизация
NumPy
X @ Y                 # BLAS, float64
X.dot(Y)              # то же самое
 
dx = x[:, None] - x[None, :]    # N×N
r3 = (dx**2 + dy**2 + dz**2)**1.5
np.fill_diagonal(r3, np.inf)
fx = (dx / r3).sum(axis=1)
01
Цикл уходит в C
операция выполняется над массивом целиком
02
BLAS
линейная алгебра на векторных инструкциях и с учётом кеша
03
Тип данных
int64 перемножается без BLAS: на 64×64 вчетверо медленнее, дальше разрыв растёт
04
Плата — память
массивы N×N для пространственного заряда
Изложение. Векторизация переносит цикл из интерпретатора в скомпилированный код NumPy: операция записывается над массивом целиком, а поэлементный цикл выполняется на C. Матричное умножение X @ Y, как и X.dot(Y), вызывает BLAS — библиотеку линейной алгебры, использующую векторные инструкции процессора и обеспечивающую эффективную работу с кешем. Тип данных определяет больше, чем запись: BLAS поддерживает только вещественные и комплексные типы, и целочисленные матрицы NumPy перемножает собственным циклом: на матрицах 64×64 на стенде вчетверо медленнее, а с ростом размера разрыв многократно увеличивается, поскольку этот цикл не использует ни блочной обработки, ни нескольких ядер. Векторизация расходует дополнительную память: пространственный заряд без циклов требует массивов разностей размером N×N, и при больших N они не помещаются в кеш, а затем и в память. Numba и Cython необходимы там, где задача не векторизуется или векторизация требует слишком много памяти.

Демонстрация:
BLAS и тип данных
01
X @ Y для float64 — 19 мкс
02
Те же размеры, int64 — 79 мкс
03
Версия NumPy и библиотека BLAS
Что наблюдается
Целочисленные матрицы NumPy перемножает без BLAS — вчетверо медленнее; чистый Python — в 900 раз дольше float64.
root@lab — запись со стенда
Демонстрация. python3 -m timeit измеряет умножение матриц 64×64 средствами NumPy. Для массивов float64 из np.random.rand оператор @ выполняется за 19 мкс: вызов уходит в BLAS. Для целочисленных массивов из np.random.randint те же размеры перемножаются за 79 мкс, вчетверо медленнее: BLAS не работает с целыми, и NumPy использует собственный цикл. Последняя команда печатает версию NumPy, 2.5.3, и имя библиотеки BLAS, с которой собран wheel-пакет NumPy, — scipy-openblas. Для сравнения: чистый Python тратит на то же умножение 17,7 мс, в 900 раз больше.

Потоки, процессы и GPUGIL, пул процессов, Python без GIL, модель CUDA
07
Причины медленности
Измерение
Приёмы языка
Память
Компиляция
NumPy и BLAS
Потоки, процессы и GPU
Изложение. Седьмая часть посвящена параллельным вычислениям. Рассмотренные выше способы, за исключением prange и многопоточной библиотеки BLAS, ускоряют работу одного ядра, тогда как виртуальной машине стенда доступны два ядра, рабочему компьютеру — от четырёх до шестнадцати, а видеокарта содержит тысячи простых вычислителей. Рассмотрим различие процессов и потоков, глобальную блокировку интерпретатора и её влияние на вычисления и ожидание, пул процессов, сборку CPython без GIL и модель вычислений CUDA.

Процессы
и потоки
import math
from concurrent.futures import ThreadPoolExecutor
 
def integrate(f, a, b, *, n_iter=1000):
    acc, step = 0, (b - a) / n_iter
    for i in range(n_iter):
        acc += f(a + i * step) * step
    return acc
 
def integrate_async(ex, f, a, b, *, n_jobs, n_iter):
    step = (b - a) / n_jobs
    fs = [ex.submit(integrate, f, a + i * step,
                    a + (i + 1) * step,
                    n_iter=n_iter // n_jobs)
          for i in range(n_jobs)]
    return sum(f.result() for f in fs)
 
with ThreadPoolExecutor(max_workers=2) as ex:
    integrate_async(ex, math.cos, 0.0, math.pi / 2,
                    n_jobs=2, n_iter=10**7)
Процесс 1
общая память
поток 1
поток 2
Процесс 2
своя память
поток 1
01
Поток
threading.Thread: start() и join(); общие данные защищают Lock и очередь queue.Queue
02
Исполнитель
ThreadPoolExecutor и ProcessPoolExecutor: submit() возвращает Future, блок with закрывает пул
03
Future
result() ожидает результат, done() проверяет готовность, add_done_callback() — обработчик
Изложение. Процесс представляет собой запущенную программу, которой операционная система выделяет собственное адресное пространство и системные ресурсы, например открытые файлы. Процессы изолированы друг от друга, и данные между ними приходится копировать. Поток исполняется внутри процесса, имеет собственный стек вызовов и разделяет с другими потоками память процесса: все потоки видят одни и те же объекты, поэтому обмен данными не требует копирования, но одновременное изменение общих данных приводит к гонкам. Порядком исполнения процессов и потоков распоряжается операционная система. В Python поток является обычным системным потоком: класс threading.Thread получает функцию и аргументы, метод start запускает поток, а join ожидает его завершения. Общие данные защищает мьютекс threading.Lock, а обмен между потоками удобнее строить на потокобезопасной очереди queue.Queue по схеме «производитель — потребитель». Принудительно завершить поток нельзя: поток получает флаг, по которому завершается сам и освобождает захваченные ресурсы.
Модуль concurrent.futures скрывает эти детали за интерфейсом исполнителя. Метод submit ставит вызов функции в очередь пула и возвращает объект Future — обёртку над ещё не завершённым вычислением: result ожидает результат, done проверяет готовность, add_done_callback регистрирует обработчик завершения. Блок with закрывает пул и дожидается всех задач. Сквозным примером части служит интеграл методом прямоугольников из главы книги «Многопоточность и GIL»: integrate суммирует значения f на равномерном разбиении отрезка, а integrate_async делит отрезок на n_jobs частей и передаёт каждую исполнителю. У ThreadPoolExecutor и ProcessPoolExecutor один интерфейс, поэтому одна и та же функция распределяет работу и по потокам, и по процессам.

GIL:
глобальная блокировка
Вычисления: потоки исполняются по очереди
поток 1
GIL
GIL
GIL
поток 2
GIL
GIL
GIL
Ожидание: на время ввода-вывода GIL свободен
поток 1
GIL
ожидание
GIL
поток 2
вычисления, GIL
@njit(nogil=True)       # Numba: без GIL
def cos_sum(a, b, n): ...
 
with nogil:             # Cython: только типы C
    for i in range(n_iter):
        acc += cos(a + i * step) * step
01
Мьютекс интерпретатора
байт-код исполняет один поток; GIL защищает счётчики ссылок
02
Переключение
ожидающий поток запрашивает GIL через 5 мс: sys.getswitchinterval()
03
Освобождение
блокирующий ввод-вывод, time.sleep, NumPy, Cython с nogil, Numba с nogil=True
Пояснение
Потоки обычной сборки CPython обеспечивают конкурентность — чередование задач, но не параллельное исполнение кода на Python.
Изложение. Глобальная блокировка интерпретатора (global interpreter lock, GIL) представляет собой мьютекс, который в каждый момент времени позволяет исполнять байт-код только одному потоку. Внутреннее состояние CPython, прежде всего счётчики ссылок у каждого объекта, не защищено от одновременного изменения, и одна общая блокировка является самым простым способом обеспечить его целостность. Поток, ожидающий GIL, по истечении интервала переключения, по умолчанию 5 мс (sys.getswitchinterval()), требует освободить блокировку, и потоки исполняются по очереди. Таким образом, в обычной сборке CPython потоки обеспечивают конкурентность — чередование задач, но не параллелизм — одновременное выполнение на нескольких ядрах.
Влияние GIL зависит от характера задачи. Если программа вычисляет на чистом Python, потоки не дают выигрыша: к прежнему времени добавляются переключения. Если программа ожидает — читает файл, получает данные по сети, опрашивает прибор, — на время блокирующего вызова поток освобождает GIL, и остальные потоки продолжают работу; для такой задачи потоков достаточно. Освободить GIL может и скомпилированный код, не обращающийся к объектам Python. NumPy делает это внутри длинных операций над массивами; в Cython для этого предусмотрен блок with nogil, в котором допустимы только типы C; Numba освобождает блокировку на время выполнения функции, скомпилированной с параметром nogil=True. Параметр parallel=True с циклом prange, рассмотренный в части о компиляции, также распределяет итерации по ядрам без участия GIL.

Демонстрация:
потоки и GIL
01
nproc: два ядра виртуальной машины
02
Чистый Python: 0,31 с в одном потоке, 0,32 в двух
03
NumPy и Numba с nogil=True: почти вдвое быстрее
04
@njit без nogil: 0,30 с в обоих случаях
05
Восемь ожиданий по 0,25 с: 2,02 → 0,25 с
Что наблюдается
Потоки ускоряют только код, освобождающий GIL: np.sin — в 1,97 раза, Numba с nogil=True — в 1,77, ожидание — в 8 раз; чистый Python — 0,99.
root@lab — запись со стенда
Демонстрация. Команда nproc подтверждает, что виртуальной машине стенда доступны два ядра, поэтому ускорение вычислений двумя потоками не может превысить двух раз. Скрипт threads.py выполняет одну и ту же работу двумя способами: обе половины подряд в одном потоке и по половине в каждом из двух потоков ThreadPoolExecutor; из трёх запусков берётся лучшее время. Интеграл на чистом Python занимает 0,31 с в одном потоке и 0,32 с в двух: байт-код исполняется по очереди, и выигрыша нет. Двадцатикратное вычисление np.sin над массивом из миллиона чисел ускоряется в 1,97 раза: цикл функции np.sin на C выполняется без GIL. Для Numba сравниваются две компиляции одной функции: с обычным @njit время в двух потоках не меняется, 0,30 с, поскольку скомпилированный код удерживает GIL, а с nogil=True сокращается до 0,17 с, в 1,77 раза. Скрипт wait.py моделирует опрос восьми приборов, каждый из которых отвечает через 0,25 с; ожидание ответа воспроизводит вызов time.sleep, освобождающий GIL так же, как блокирующий ввод-вывод. Опрос по очереди занимает 2,02 с, восемью потоками — 0,25 с: ожидания перекрываются, и ускорение в 8 раз не ограничено числом ядер.

Пул
процессов
import math
from concurrent.futures import ProcessPoolExecutor
 
from integ import integrate_async  # из модуля
 
if __name__ == "__main__":         # модуль импортируется
    a, b = 0.0, math.pi / 2
    with ProcessPoolExecutor(max_workers=2) as ex:
        r = integrate_async(ex, math.cos, a, b,
                            n_jobs=2, n_iter=10**7)
        f = ex.submit(lambda x: x * x, 3)
        f.result()                 # PicklingError
Важно
Функция и аргументы передаются через pickle, функция — по имени: lambda не передаётся, а функцию из ячейки блокнота при spawn и forkserver рабочий процесс не находит.

МетодНовый процессПо умолчанию
forkкопия родительского процессаLinux по 3.13
spawnновый интерпретатор, импорт главного модуляWindows, macOS с 3.8
forkserverкопия заранее запущенного сервераLinux с 3.14
Пояснение
Запуск процессов и пересылка данных требуют времени: выигрыш дают крупные независимые задачи, а пул создаётся один раз.
Изложение. Поскольку GIL один на интерпретатор, наиболее прямой способ его обхода состоит в запуске нескольких интерпретаторов: каждый процесс имеет собственный GIL, и процессы вычисляют одновременно на разных ядрах. Модуль multiprocessing устроен подобно threading, а ProcessPoolExecutor из concurrent.futures имеет тот же интерфейс, что и пул потоков, поэтому integrate_async переходит на процессы заменой одного класса. Платой является изоляция: процессы не разделяют память, функция и её аргументы передаются рабочему процессу сериализацией модулем pickle, а результат возвращается тем же путём. Функция передаётся по имени модуля и функции, поэтому в пуле выполняются функции, определённые на верхнем уровне импортируемого модуля: lambda и вложенная функция не сериализуются: для lambda pickle возбуждает ошибку PicklingError, для вложенной функции до Python 3.14 — AttributeError, а функцию, определённую в интерактивном сеансе или в ячейке блокнота, рабочий процесс, запущенный методом spawn или forkserver, не находит, поскольку главный модуль такого сеанса импортировать нельзя. Ошибка сериализации не возникает в момент вызова submit: она сохраняется в объекте Future и возбуждается вызовом result(), поэтому результат каждой задачи необходимо запросить.
Способ создания рабочего процесса задаёт метод запуска. Метод fork копирует родительский процесс целиком, spawn запускает новый интерпретатор и импортирует главный модуль, forkserver порождает процессы копированием заранее запущенного однопоточного сервера. По умолчанию в Windows и, начиная с Python 3.8, в macOS используется spawn, в Linux до версии 3.13 включительно — fork, а начиная с 3.14 — forkserver: копия процесса, в котором уже работают потоки, может унаследовать захваченные ими блокировки. При spawn и forkserver главный модуль импортируется в рабочем процессе заново, поэтому создание пула необходимо поместить под проверку if __name__ == "__main__". Запуск процессов и пересылка данных требуют времени, поэтому выигрыш дают крупные независимые задачи, а пул создаётся один раз и используется повторно. В Python 3.14 появился промежуточный вариант — InterpreterPoolExecutor: его рабочие потоки исполняют задачи в изолированных подынтерпретаторах одного процесса, у каждого из которых собственный GIL, а функция, аргументы и результат копируются между интерпретаторами: функции без замыканий и глобальных имён передаются напрямую, остальные — через pickle.

Метод	Новый процесс	По умолчанию
fork	копия родительского процесса	Linux по 3.13
spawn	новый интерпретатор, импорт главного модуля	Windows, macOS с 3.8
forkserver	копия заранее запущенного сервера	Linux с 3.14

Демонстрация:
пул процессов
01
Python 3.12, метод fork: 0,31 → 0,18 с
02
Python 3.14, метод forkserver: 0,39 → 0,21 с
03
Запуск пула: 2 мс при fork, 30 мс при forkserver
04
lambda в пуле процессов: PicklingError
Что наблюдается
Два процесса на двух ядрах ускоряют чистый Python в 1,7–1,9 раза; forkserver запускает пул на порядок дольше fork; lambda не сериализуется.
root@lab — запись со стенда
Демонстрация. Скрипт procs.py вычисляет тот же интеграл в одном процессе и в пуле ProcessPoolExecutor из двух процессов, а также печатает версию Python и метод запуска по умолчанию. На Python 3.12 в Linux пул создаёт процессы методом fork: один процесс считает 0,31 с, два — 0,18 с, в 1,7 раза быстрее, а запуск пула занимает 2 мс. Второй интерпретатор стенда, Python 3.14 (его особенность — отсутствие GIL — рассматривается на следующем слайде), использует новый метод по умолчанию, forkserver: расчёт ускоряется с 0,39 до 0,21 с, в 1,9 раза, а запуск пула занимает 30 мс, на порядок дольше, поскольку сначала запускается процесс-сервер. Ускорение на двух ядрах не превышает двух раз. Скрипт lam.py передаёт пулу встроенную функцию pow, которая сериализуется по имени и выполняется, а затем lambda: pickle не находит её по имени в модуле __main__, и вызов завершается ошибкой PicklingError.

Python
без GIL

ВерсияСборка без GIL
3.13экспериментальная, PEP 703: python3.13t
3.14поддерживаемая, PEP 779: python3.14t
$ python3.14t -VV
Python 3.14.8 free-threading build ...
>>> import sys; sys._is_gil_enabled()
False
$ python3.14t -X gil=1 app.py      # с GIL
$ PYTHON_GIL=1 python3.14t app.py  # то же
Пояснение
Обычная сборка с GIL остаётся сборкой по умолчанию; свободнопоточная устанавливается рядом.
01
Потоки на всех ядрах
байт-код исполняется параллельно; list, dict и set защищены своими блокировками
02
Гонки остаются
counter += 1 не атомарна: общие данные по-прежнему защищает Lock
03
Расширения на C
объявляют поддержку; модуль без объявления снова включает GIL
04
Цена
однопоточный код медленнее: в 3.13 — на 40 %, в 3.14 — на 1–8 %
Изложение. Начиная с Python 3.13 интерпретатор CPython собирается и без глобальной блокировки (PEP 703). Такая сборка устанавливается рядом с обычной, её исполняемый файл имеет суффикс t — python3.13t, python3.14t, — а потоки в ней исполняют байт-код одновременно на всех ядрах. В 3.13 сборка была экспериментальной; в 3.14 она официально поддерживается (PEP 779), но сборкой по умолчанию не является. Работает ли GIL, позволяет узнать функция sys._is_gil_enabled(); в свободнопоточной сборке блокировку можно включить ключом -X gil=1 или переменной окружения PYTHON_GIL=1.
Без GIL интерпретатор защищает своё состояние иначе: встроенные типы — список, словарь, множество — снабжены собственными блокировками, поэтому их одновременное изменение не нарушает целостности объектов. Атомарными операции программы от этого не становятся: составное действие над общими данными, например counter += 1, состоит из чтения, сложения и записи, и гонки в программе по-прежнему устраняет Lock. Расширения на C должны объявить, что не полагаются на GIL; при импорте модуля без такого объявления интерпретатор снова включает блокировку и выводит предупреждение. Однопоточный код в свободнопоточной сборке выполняется медленнее. В 3.13, где специализирующий адаптивный интерпретатор в этой сборке был отключён, замедление составляло в среднем около 40 процентов на наборе тестов pyperformance; в 3.14 специализация включена, и по документации замедление составляет от 1 процента на macOS с процессором aarch64 до 8 процентов на Linux с x86-64; расход памяти при этом выше. Пакеты wheel для свободнопоточной сборки NumPy выпускает с версии 2.1, Numba — с версии 0.65, и обе библиотеки называют эту поддержку экспериментальной. Таким образом, сборка без GIL позволяет распараллелить потоками код на чистом Python, но требует и потокобезопасной программы, и подготовленных расширений.

Версия	Сборка без GIL
3.13	экспериментальная, PEP 703: python3.13t
3.14	поддерживаемая, PEP 779: python3.14t

Демонстрация:
Python без GIL
01
python3.14t -VV: free-threading build
02
-X gil=1: два потока — 0,39 с, как один
03
-X gil=0: два потока — 0,20 с вместо 0,38
04
Счётчик без Lock: 1 118 818 из 2 000 000
Что наблюдается
Без GIL чистый Python в двух потоках на двух ядрах ускоряется в 1,9 раза; счётчик без Lock теряет 44 % увеличений, с Lock результат точный.
root@lab — запись со стенда
Демонстрация. На стенде установлен интерпретатор CPython 3.14.8 в свободнопоточной сборке, без сторонних пакетов; ключ -VV печатает пометку free-threading build. Скрипт ft.py вычисляет интеграл на чистом Python пулом из одного и из двух потоков и печатает значение sys._is_gil_enabled(). С ключом -X gil=1 блокировка включена, и два потока считают столько же, сколько один, — 0,39 с. С ключом -X gil=0 тот же код в двух потоках выполняется за 0,20 с вместо 0,38 с, в 1,9 раза быстрее: потоки исполняют байт-код одновременно на двух ядрах. Отметим, что один поток считает здесь 0,38 с, тогда как Python 3.12 в демонстрации «пул процессов» — 0,31 с; это сравнение смешивает версию интерпретатора и параметры его сборки и цену отказа от GIL не измеряет: оценка 1–8 % относится к обычной и свободнопоточной сборкам одной версии на наборе тестов pyperformance. Скрипт race.py показывает, что отсутствие GIL не делает программу потокобезопасной: два потока увеличивают общий счётчик по миллиону раз, и без Lock итог составляет 1 118 818 вместо 2 000 000. Чтение, сложение и запись одного потока перемежаются с теми же действиями другого, и 44 процента увеличений теряются; с Lock результат точный.

Видеокарта:
SIMT и warp

CPUGPU
ядра4–64 сложныхтысячи простых
оптимизирован подзадержку одной задачипропускную способность
задержку памяти скрываютбольшие кешитысячи нитей
полоса памятидесятки ГБ/ссотни ГБ/с — единицы ТБ/с
нагрузкаветвистая, последовательнаяоднотипная над массивами
Для физика
Частицы, узлы сетки, события детектора — одна формула над миллионами независимых элементов.
01
Warp и SIMT
32 нити выполняют одну инструкцию над разными данными
02
Дивергенция
ветви if внутри warp выполняются по очереди
03
Иерархия памяти
регистры нити, разделяемая память блока, глобальная память карты
04
Объединение обращений
запросы warp сливаются в 32-байтные транзакции: соседним нитям — соседние адреса
Изложение. Процессор оптимизирован под задержку: одна цепочка инструкций должна выполниться как можно быстрее, для чего ядро располагает предсказателем переходов, внеочередным исполнением и большими кешами. Видеокарта оптимизирована под пропускную способность: тысячи простых ядер выполняют максимальный объём однотипной работы в секунду, а задержку памяти скрывают не кешами, а числом нитей — пока одни ожидают данные, вычисляют другие. Нити исполняются группами по 32, называемыми warp; все нити warp в каждый момент выполняют одну инструкцию над разными данными, и эта модель называется SIMT (Single Instruction, Multiple Threads). Если внутри warp код ветвится, обе ветви выполняются по очереди, а нити другой ветви простаивают; это явление называется дивергенцией, и сильно ветвящийся код на видеокарте ускоряется плохо.
Память видеокарты устроена ступенями: регистры отдельной нити, разделяемая память блока нитей объёмом порядка ста килобайт на мультипроцессор и глобальная память карты объёмом в гигабайты, задержка которой достигает сотен тактов. Обращения нитей одного warp к глобальной памяти объединяются в минимальное число 32-байтных транзакций: если соседние нити читают соседние четырёхбайтные числа, 32 нитям достаточно четырёх транзакций, а при чтении с большим шагом каждой нити требуется собственная транзакция, и полоса памяти расходуется впустую. Поэтому данные на видеокарте располагают так, чтобы соседние нити обращались к соседним элементам массива. Для физика это означает, что видеокарта подходит для задач массового параллелизма данных: движение миллионов частиц, пересчёт узлов разностной сетки, обработка событий детектора.

	CPU	GPU
ядра	4–64 сложных	тысячи простых
оптимизирован под	задержку одной задачи	пропускную способность
задержку памяти скрывают	большие кеши	тысячи нитей
полоса памяти	десятки ГБ/с	сотни ГБ/с — единицы ТБ/с
нагрузка	ветвистая, последовательная	однотипная над массивами

Модель
CUDA
Сетка (grid): 3 блока по 4 нити, n = 10
блок 0: i = 0…3
0
1
2
3
блок 1: i = 4…7
0
1
2
3
блок 2: i = 8…11
0
1
2
3
i = blockIdx.x · blockDim.x + threadIdx.x
Нить 1 блока 2: i = 2 · 4 + 1 = 9; серые нити 2 и 3 блока 2 простаивают: i ≥ n
Host
CPU и оперативная
память
PCIe
Device
GPU и память
карты
from numba import cuda
@cuda.jit
def add_kernel(a, b, out):
    i = cuda.grid(1)        # глобальный индекс
    if i < out.size:        # лишние нити не работают
        out[i] = a[i] + b[i]
 
d_a, d_b = cuda.to_device(a), cuda.to_device(b)
d_out = cuda.device_array_like(d_a)
tpb = 256                   # нитей в блоке
bpg = (n + tpb - 1) // tpb  # блоков
add_kernel[bpg, tpb](d_a, d_b, d_out)
out = d_out.copy_to_host()  # device → host
Пояснение
Встроенный numba.cuda устарел с Numba 0.61; NVIDIA развивает numba-cuda-mlir, numba-cuda только сопровождается.
Изложение. CUDA представляет собой платформу NVIDIA для вычислений общего назначения на видеокарте, и её терминологией пользуются все инструменты для GPU. Host объединяет процессор и оперативную память, device — видеокарту и её память; адресные пространства у них разные, и данные копируются между ними явно. Ядром (kernel) называется функция, которую host запускает на device: её одновременно исполняют тысячи нитей, каждая со своим элементом данных. Нити организованы в сетку блоков. Блок целиком выполняется на одном мультипроцессоре, и его нити могут обмениваться данными через разделяемую память, а блоки друг от друга не зависят и распределяются по мультипроцессорам в произвольном порядке.
Глобальный индекс нити получим из номера блока, размера блока и номера нити в блоке: i = blockIdx.x · blockDim.x + threadIdx.x; например, нить 5 в блоке 3 при блоках по 256 нитей обрабатывает элемент 773. Привычный цикл for исчезает: вместо него запускается n нитей. В Numba ядро записывается функцией с декоратором @cuda.jit; cuda.grid(1) возвращает глобальный индекс, а проверка i < out.size необходима, поскольку число нитей кратно размеру блока и может превышать число элементов. Конфигурация запуска — число блоков и число нитей в блоке — указывается в квадратных скобках; размер блока выбирают кратным 32, обычно от 128 до 512 нитей. Ядро ничего не возвращает и записывает результат в переданный массив, а cuda.to_device и copy_to_host явно копируют данные между host и device. С версии Numba 0.61 встроенный модуль numba.cuda объявлен устаревшим, и поддержка CUDA перешла в пакет NVIDIA numba-cuda. В 2026 году numba-cuda переведён в режим сопровождения: новые возможности NVIDIA развивает в пакете numba-cuda-mlir, совместимом с ядрами numba-cuda, и для большинства программ переход на него сводится к замене импорта from numba import cuda на from numba_cuda_mlir import cuda; код, использующий API расширения Numba, требует переработки, а сам пакет поддерживает видеокарты вычислительной способности не ниже 7.0.

Данные
на видеокарте
import cupy as cp       # numpy — тот же код на CPU
 
h = q * dt / (2 * m)              # схема Бориса
t = h * B0
s = 2 * t / (1 + t * t)
rng = cp.random.default_rng(42)   # сразу на GPU
x = rng.standard_normal(N, dtype=cp.float32) * 1e-3
y = rng.standard_normal(N, dtype=cp.float32) * 1e-3
vx = rng.standard_normal(N, dtype=cp.float32) * 1e5
vy = rng.standard_normal(N, dtype=cp.float32) * 1e5
y0 = float(y.mean())              # центр облака
for step in range(n_steps):
    vx += h * E0                  # полшага поля E
    px, py = vx + vy * t, vy - vx * t
    vx, vy = vx + py * s, vy - px * s  # поворот
    vx += h * E0
    x += vx * dt
    y += vy * dt
v = (float(y.mean()) - y0) / (n_steps * dt)  # на host
01
Шина и память карты
PCIe x16 — 16–64 ГБ/с; память GPU — до нескольких ТБ/с
02
Оценка для 1 ГБ
пересылка при 16 ГБ/с — около 60 мс, y = 2x на карте — около 2 мс
03
Асинхронный запуск
вызов ставит ядро в очередь; замер — после прогрева и synchronize()
04
CuPy
интерфейс NumPy на GPU; cp.asnumpy и cp.asarray копируют данные
Важно
Данные создаются на карте и остаются там до конца расчёта; на host возвращается результат.
Изложение. Видеокарта соединена с процессором шиной PCIe, полоса которой составляет десятки гигабайт в секунду, тогда как память самой карты передаёт от сотен гигабайт до единиц терабайт в секунду: в примере из руководства NVIDIA по оптимизации — 16 ГБ/с у шины PCIe x16 третьего поколения против 898 ГБ/с у памяти Tesla V100. Оценим цену пересылки на примере гигабайта данных: копирование через PCIe третьего поколения занимает около 60 мс, а поэлементная операция y = 2x на карте прочитает и запишет два гигабайта примерно за 2 мс. Таким образом, при малом числе операций на байт пересылка обходится дороже вычисления, и главное правило вычислений на GPU состоит в том, чтобы создать данные на карте один раз, выполнить там всю цепочку шагов и вернуть только результат.
Самый короткий путь к видеокарте из Python обеспечивает CuPy — реализация интерфейса NumPy поверх CUDA: массивы размещаются в памяти карты, операции над ними выполняет GPU, а cp.asnumpy и cp.asarray копируют данные в одну и в другую сторону. Пример на слайде, взятый из главы книги «CUDA и вычисления на GPU», моделирует дрейф электронов в скрещенных полях E = (E0, 0, 0) и B = (0, 0, B0) схемой Бориса. Массивы создаются генератором cp.random сразу на карте, каждая арифметическая операция над массивами выполняется отдельным ядром над всеми частицами, а на host возвращается одно число — скорость дрейфа, которая по теории равна по модулю E0/B0, направлена против оси y и не зависит ни от заряда, ни от массы. Запуск ядер асинхронный: вызов ставит работу в очередь и возвращается сразу, поэтому при замере времени необходимо дождаться завершения вызовом synchronize или воспользоваться функцией cupyx.profiler.benchmark, а первый вызов, компилирующий ядра, в замер не включают. Замена первой строки на import numpy as cp даёт версию для процессора, пригодную для отладки на машине без видеокарты.

Демонстрация:
ядро CUDA на CPU
01
cuda.is_available(): видеокарты на стенде нет
02
Симулятор CUDA: ядро на CPU, 3 блока по 4 нити
03
Индекс i = блок · 4 + нить; нити с i = 10 и 11 простаивают
04
Дрейф 10⁵ электронов: NumPy вместо CuPy
Что наблюдается
Симулятор проверяет индексы, а не скорость; скорость дрейфа −9,998·10⁴ м/с совпадает с −E0/B0 с точностью 0,02 %.
root@lab — запись со стенда
Демонстрация. На стенде видеокарты нет: cuda.is_available() возвращает False, и запуск ядра без симулятора завершается ошибкой поиска драйвера CUDA. Переменная окружения NUMBA_ENABLE_CUDASIM=1 включает симулятор CUDA, входящий в Numba: ядро выполняется на процессоре интерпретатором Python, блоки — по одному, а нити блока — потоками Python. Симулятор предназначен для отладки и проверки индексов, операций уровня warp не поддерживает и о скорости вычислений на видеокарте ничего не говорит. Скрипт kernel.py запускает ядро сложения векторов для 10 элементов с блоками по 4 нити: блоков получается 3, нитей — 12, и каждая нить записывает номер своего блока и свой номер в блоке. Таблица подтверждает формулу индекса: элемент 9 обрабатывает нить 1 блока 2, 2 · 4 + 1 = 9, а нити 2 и 3 последнего блока не проходят проверку i < out.size и простаивают. Скрипт drift.py повторяет пример с дрейфом, в котором вместо CuPy импортирован NumPy: 10⁵ электронов за 2144 шага рассчитываются на процессоре за 0,3 с, и скорость дрейфа −9,9983·10⁴ м/с отличается от теоретической −E0/B0 = −10⁵ м/с на 0,02 %. Число шагов выбрано равным шести ларморовским периодам: облако начинает движение со средней скоростью, равной нулю, его центр движется по циклоиде, и средняя скорость совпадает со скоростью дрейфа только за целое число периодов.

Итог
замеров

ЗадачаСредствоВыигрыш на стенде
цикл Python над спискамиприёмы языкав 2,8 раза
числовой циклNumba и Cythonв 217 и 165 раз
линейная алгебраNumPy и BLASв 900 раз
цикл Numba, операции NumPyprange, nogil=True, потокив 1,8–2 раза
чистый Python на всех ядрахпул процессов, python3.14tв 1,7–1,9 раза
ожидание ввода-выводапотоки, asyncioв 8 раз при восьми ожиданиях
массовый параллелизм данныхGPU: CuPy, numba-cuda-mlirвидеокарты на стенде нет
Важно
На двух ядрах стенда параллелизм дал не больше двух раз, компиляция и векторизация — сотни.
Для физика
Перебор параметров установки — естественная задача для пула процессов: расчёты независимы.
Изложение. Итог измерений лекции сведён в таблицу по задачам. Приёмы языка ускорили умножение матриц на списках в 2,8 раза; Numba после прогрева — в 217 раз, Cython с объявленными типами — в 165 раз по сравнению с чистым Python на списках, NumPy с BLAS — в 900 раз. Параллельные средства ограничены числом ядер: prange в Numba ускорил расчёт пространственного заряда почти вдвое, потоки с NumPy и с Numba nogil=True — в 1,8–2 раза, а пул из двух процессов и свободнопоточная сборка ускорили интеграл на чистом Python в 1,7–1,9 раза. Для ожидания ввода-вывода достаточно потоков, и число ядер их не ограничивает: восемь ожиданий перекрылись, и опрос ускорился в 8 раз. Видеокарты на стенде нет, поэтому её выигрыш не измерялся; он определяется размером задачи, числом операций на байт данных и тем, остаются ли данные на карте. Таким образом, к параллельным вычислениям переходят после того, как горячий цикл ускорен на одном ядре: два ядра дают не больше двух раз, тогда как компиляция и векторизация — сотни. Подробно эти темы разобраны в главах книги «Многопоточность и GIL», «Асинхронность» и «CUDA и вычисления на GPU».
Для физиков: перебор параметров установки или ансамбль независимых расчётов — естественная задача для пула процессов: каждый расчёт выполняется в своём процессе, а результаты собираются в конце.

Задача	Средство	Выигрыш на стенде
цикл Python над списками	приёмы языка	в 2,8 раза
числовой цикл	Numba и Cython	в 217 и 165 раз
линейная алгебра	NumPy и BLAS	в 900 раз
цикл Numba, операции NumPy	prange, nogil=True, потоки	в 1,8–2 раза
чистый Python на всех ядрах	пул процессов, python3.14t	в 1,7–1,9 раза
ожидание ввода-вывода	потоки, asyncio	в 8 раз при восьми ожиданиях
массовый параллелизм данных	GPU: CuPy, numba-cuda-mlir	видеокарты на стенде нет

Итоги
лекции
01
Медленность Python обусловлена динамической типизацией, изменяемостью и циклом интерпретатора над объектами.
02
Ускоряют только измеренное: timeit для фрагмента, cProfile и line_profiler для программы, py-spy для работающего процесса.
03
Структура данных и алгоритм важнее микрооптимизаций; любой совет проверяется замером на своей версии Python.
04
Память измеряют tracemalloc; числа хранят в массивах, объекты — со __slots__.
05
NumPy, Numba и Cython ускоряют горячий цикл в десятки и сотни раз.
06
Потоки ускоряют ожидание и код без GIL; Python на всех ядрах — процессы или сборка без GIL; работа над массивами — видеокарта.
Изложение. Подведём итог. Python работает медленно не потому, что плохо написан, а вследствие своей гибкости: тип каждого значения выясняется при выполнении, почти всё может быть изменено по ходу работы программы, а цикл интерпретатора работает с объектами в куче. Ускорять имеет смысл только измеренное: timeit отвечает на вопрос, какой из двух фрагментов быстрее, cProfile и line_profiler — где программа проводит время, py-spy — чем занят уже работающий расчёт. Самый дешёвый выигрыш даёт смена структуры данных или алгоритма; микрооптимизации дают проценты, а рекомендации устаревают с выходом новых версий интерпретатора, поэтому каждая проверяется замером. Память измеряют tracemalloc, а не getsizeof; числа хранят в массивах, а многочисленные объекты — со __slots__. NumPy, Numba и Cython переносят горячий цикл в машинный код и ускоряют его в десятки и сотни раз. Потоки в обычной сборке CPython ускоряют ожидание ввода-вывода и код, освобождающий GIL, — NumPy, Numba с nogil=True, Cython с блоком with nogil; вычисления на чистом Python распределяют по ядрам пул процессов или свободнопоточная сборка Python 3.14, а однотипную работу над миллионами элементов выполняет видеокарта, если данные остаются в её памяти. Таким образом, общее правило остаётся прежним: сначала правильно, потом быстро, и только затем параллельно.

Упражнения
для самопроверки
01
Профиль своего кода
cProfile и line_profiler на программе из проекта: какая функция и строка горячие
02
Пространственный заряд
версия со сглаживанием r² + ε² на Numba; проверка закона N² и выигрыша prange
03
Память частиц
список объектов, __slots__ и массив структур NumPy для миллиона частиц
04
Векторизация блоками
силы без массивов N×N: разбиение на блоки и сравнение с Numba
05
Параллельный перебор
ансамбль расчётов в пуле процессов и в потоках python3.14t: ускорение и число ядер
06
Дрейф на видеокарте
пример с CuPy в Colab: время с пересылкой массивов на каждом шаге и без неё
Пояснение
Каждое ускорение сравнивается с исходной версией на тех же данных и проверяется тестом.
Для физика
Проект курса: горячий цикл расчёта — первый кандидат на профиль и Numba.
Изложение. Отдельного задания к лекции нет: упражнения служат самопроверке. Первое — снять профиль программы из собственного проекта cProfile и line_profiler и найти горячие функцию и строку. Второе — добавить к расчёту пространственного заряда сглаживание, заменив r² на r² + ε², скомпилировать функцию Numba и проверить, что время растёт как N², а prange ускоряет расчёт не больше числа ядер. Третье — сравнить память под миллион частиц в виде списка объектов, объектов со __slots__ и массива структур NumPy. Четвёртое — векторизовать расчёт сил блоками, без массивов N×N целиком, и сравнить время и память с версией Numba. Пятое — распределить ансамбль независимых расчётов, например интегралов с разными параметрами, по пулу процессов и по потокам свободнопоточного интерпретатора python3.14t и сопоставить ускорение с числом ядер. Шестое — выполнить пример с дрейфом частиц на CuPy в Google Colab и сравнить время расчёта, при котором массивы остаются на видеокарте, со временем расчёта, при котором они на каждом шаге копируются на host и обратно. Каждое ускорение сравнивается с исходной версией на тех же данных и проверяется тестом.

Источники1
Книга курса: раздел «Ускорение расчётов», главы об оптимизации, скорости, GIL и CUDA
phys-dev.github.io/soft-dev-book
2
Лебедев С. Быстрее, Python, ещё быстрее. – Computer Science Center, 2015
youtube.com/watch?v=-lMiAKKyLFI
3
Лебедев С. Многопоточность и GIL. – Computer Science Center, 2015
youtube.com/watch?v=nR8WhdcRJwM
4
timeit, cProfile, tracemalloc, concurrent.futures, multiprocessing, free threading — документация Python
docs.python.org/3/howto/free-threading-python.html
5
py-spy — документация; line_profiler — kernprof.readthedocs.io
github.com/benfred/py-spy
6
Numba — документация; numba-cuda и numba-cuda-mlir — nvidia.github.io; CuPy — docs.cupy.dev
numba.readthedocs.io
7
Cython — документация
cython.readthedocs.io
8
Knuth D. Structured Programming with go to Statements. – ACM Computing Surveys, 1974
dl.acm.org/doi/10.1145/356635.356640
9
Gorelick M., Ozsvald I. High Performance Python. – 3rd ed. – O'Reilly, 2025
Изложение. Основой лекции служат главы раздела книги курса «Ускорение расчётов»; сквозной пример с умножением матриц и порядок изложения инструментов следуют лекции «Быстрее, Python, ещё быстрее» курса Python Computer Science Center 2015 года, а часть о потоках и процессах — лекции «Многопоточность и GIL» того же курса. Для самостоятельного изучения рекомендуются документация модулей timeit, cProfile, tracemalloc, concurrent.futures и multiprocessing, руководство по свободнопоточной сборке Python, документация line_profiler, py-spy, Numba, numba-cuda и numba-cuda-mlir, Cython и CuPy, статья Д. Кнута 1974 года и книга «High Performance Python».

Спасибо за вниманиеВопросы
phys-dev.github.io/soft-dev-book
@njit
Изложение. Надпись на последнем слайде — @njit, декоратор, который одной строкой позволяет ускорить горячий цикл в десятки и сотни раз, если до этого профилировщик показал, что ускорения требует именно этот цикл. Лекция завершается ответами на вопросы.