Причудливые алгоритмы, пошаговое руководство по исходному коду и возможные улучшения

В части I мы рассмотрели основы выборки с обратным преобразованием (ITS) и создали собственную реализацию ITS на чистом python для выборки чисел из стандартного нормального распределения. Затем мы сравнили скорость нашей несколько оптимизированной функции со скоростью встроенной функции SciPy и обнаружили, что нам чего-то не хватает — например, 40x медленнее.

В этой части цель состоит в том, чтобы объяснить, почему это так, копаясь в соответствующих битах кодовой базы SciPy и NumPy, чтобы увидеть, где эти улучшения скорости проявляются. В целом мы обнаружим, что он состоит из комбинации:

  • более быстрые функции либо из-за того, что они написаны на Cython, либо прямо на C
  • Более быстрые новые алгоритмы выборки по сравнению с нашей проверенной и проверенной выборкой с обратным преобразованием

Как мы генерируем нормально распределенные случайные выборки в SciPy?

Ниже приведен код для генерации 1,000,000 случайных чисел из стандартного нормального распределения.

43.5 ms ± 1.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Таким образом, функция rvs генерирует 1,000,000 отсчетов чуть более чем за 40ms. Для сравнения, мы смогли добиться этого в среднем за 2.3s, используя наш алгоритм, основанный на принципе выборки с обратным преобразованием. Чтобы понять разницу в скорости, нам нужно погрузиться в этот метод rvs.

Стоит отметить, что (вообще) в SciPy ядро ​​логики содержится в методах подчеркивания — поэтому, когда мы хотим заглянуть в rvs, на самом деле мы хотим увидеть код для _rvs. Методы без подчеркивания обычно реализуют некоторую проверку типа аргумента или установку по умолчанию перед передачей методам подчеркивания.

Прежде чем приступить к работе, давайте просто сделаем краткий обзор того, как SciPy организует функциональность распространения в библиотеке.

rv_generic и rv_continuous

Дистрибутивы SciPy создаются из аккуратной структуры наследования с:

  • rv_generic как класс верхнего уровня, предоставляющий такие методы, как get_support и mean
  • rv_continuous и rv_discrete наследуя от него более специфичные методы

Таким образом, в приведенном выше случае, когда мы инициировали наш класс нормального распределения snorm как stats.norm(), на самом деле он создает экземпляр rv_continuous, который наследует множество функций от rv_generic. Чтобы быть еще более конкретным, мы фактически создаем экземпляр rv_frozen, который является версией rv_continuous, но с фиксированными параметрами распределения (например, средним значением и дисперсией). Имея это в виду, давайте теперь заглянем внутрь метода rvs.

внедорожники

Когда мы запускаем магию ?? для snorm.dist._rvs, мы видим следующий фрагмент кода:

# ?? snorm.dist._rvs
def _rvs(self, size=None, random_state=None):
    return random_state.standard_normal(size)

Таким образом, кажется, что где-то в созданном нами классе распределения мы где-то присвоили объект random_state, и этот объект random_state содержит метод, который может возвращать числа, распределенные согласно стандартному нормальному распределению.

Оказалось, что объект random_state, который выдает эти случайные числа, на самом деле принадлежит NumPy. Мы видим это, глядя на исходный код rv_generic, который содержит в своем методе __init__ вызов утилиты SciPy. метод под названием check_random_state, который, если начальное число уже не передано, установит random_state как экземпляр np.random.RandomState. Ниже приведен этот фрагмент кода:

Переходим к NumPy

Таким образом, кажется, что волшебство, которое обеспечивает такую ​​молниеносную быструю выборку, на самом деле находится в NumPy, а не в SciPy. Это не должно быть таким уж шокирующим, поскольку SciPy намеренно построен на основе NumPy, чтобы предотвратить дублирование и несоответствия, когда две библиотеки могут предоставлять идентичные функции. Об этом прямо сказано в первой строке документации SciPy Intro здесь:

"SciPy — это набор математических алгоритмов и удобных функций, созданных на основе расширения Python NumPy".

Чтобы увидеть, что происходит, мы можем взглянуть на класс np.random.RandomState здесь. Мы можем видеть из использования:

  • cdef вместо def для объявления функции
  • расширение файла .pyx вместо .py

оба из которых указывают на то, что функция написана с использованием Cython — языка, очень похожего на Python, который позволяет писать функции почти с синтаксисом Python, но затем для эффективности компилируется в оптимизированный код C/C++. Как они сами выразились в документации:

«Исходный код преобразуется в оптимизированный код C/C++ и компилируется в виде модулей расширения Python. Это обеспечивает как очень быстрое выполнение программы, так и тесную интеграцию с внешними библиотеками C, сохраняя при этом высокую производительность программиста, которой хорошо известен язык Python».

В этом классе есть две вещи, на которые нам нужно обратить внимание, чтобы понять процесс выборки:

  • что он делает для генерации равномерно распределенных случайных чисел (PRNG)
  • какой алгоритм он использует для преобразования этих равномерно распределенных чисел в нормально распределенные числа

ГПСЧ

Как упоминалось в части I, создание случайной выборки требует некоторой формы случайности. Почти всегда это не настоящая случайность, а ряд чисел, сгенерированных генератором псевдослучайных чисел (PRNG). Как и в случае с алгоритмами выборки, существует множество доступных PRNG, и конкретная реализация, используемая здесь, подробно описана в методе __init__ из np.random.RandomState:

Как показано выше, когда класс инициируется, PRNG по умолчанию устанавливается как реализация алгоритма Вихрь Мерсенна — названный таким образом, что он имеет длину периода простое число Мерсенна (количество случайных чисел, которые он может генерироваться до того, как он начнет повторяться).

Процесс выборки

Где-то дальше по коду класса np.random.RandomState мы видим определение standard_normal, делающее вызов чего-то под названием legacy_gauss. Код C для функции legacy_gauss находится здесь, и для простоты просмотра мы покажем его здесь:

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

Резюме

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

  • функция SciPy под названием _rvs, написанная на python, инициирует
  • класс NumPy np.random.RandomState, написанный на Cython, который
  • генерирует равномерно распределенные числа, используя алгоритм Mersenne Twister, а затем
  • передает эти числа в функцию legacy_gauss, написанную на C, которая производит нормально распределенные выборки, используя полярный метод Marsaglia.

Вышеизложенное показывает, на что пошли умные люди, создающие SciPy и NumPy, для создания эффективного кода. У нас есть верхний уровень, вызываемый пользователями (такими как вы и я), который написан на python (для «продуктивности программиста» python), прежде чем более глубокие уровни инфраструктуры будут все чаще писаться как можно ближе к C (для скорости).

Почему SciPy вызывает функцию NumPy, которая считается «устаревшей»?

Потому что выборка — это отрасль математики/информатики, которая все еще движется вперед. В отличие от других областей, где определенные принципы были согласованы несколько столетий назад и с тех пор не претерпели изменений, эффективная выборка различных дистрибутивов все еще претерпевает новые изменения. По мере тестирования новых разработок мы хотели бы обновить наши процессы по умолчанию, чтобы включить эти улучшения.

Именно это произошло в июле 2019 года с NumPy 1.17.0, когда они представили 2 новые функции, влияющие на выборку:

Однако из-за стремления к обратной совместимости PRNG вместо внесения критического изменения они представили новый способ инициирования PRNG и переключили старый способ на ссылку на «устаревший» код.

Упомянутая здесь обратная совместимость — это желание функции PRNG генерировать одну и ту же строку случайных чисел с одним и тем же начальным числом. Два разных алгоритма не будут давать одинаковые случайные числа, даже если им задано одно и то же начальное число. Эта воспроизводимость особенно важна для тестирования.

Похоже, SciPy еще не был обновлен, чтобы использовать эти новые разработки.

Можем ли мы победить SciPy?

Учитывая то, что мы знаем сейчас о том, как в SciPy реализована выборка нормального распределения, можем ли мы превзойти ее?

Ответ положительный — используя последние разработки в области выборки, реализованные для нас в NumPy. Ниже приведена реализация выборки, в которой мы:

  • использовать последний PRNG
  • использовать новый алгоритм зиккурата для преобразования этих чисел в нормально распределенную выборку
# test scipy speed
%timeit snorm.rvs(size=n)
51 ms ± 5.08 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
# test numpy speed
%timeit nnorm.normal(size=n)
24.3 ms ± 1.84 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Таким образом, кажется, что сейчас мы примерно в 2x быстры, чем SciPy — что-то, что находится в ожидаемом диапазоне 2-10x, как подчеркивает NumPy в своем выпуске здесь.

Вывод: насколько это полезно?

Когда дело доходит до реализации пользовательской выборки распределения: очень полезно. Теперь мы полностью понимаем решение использовать скорость выборки, подобную SciPy, и можем соответствующим образом реализовать пользовательскую выборку распределения. Мы можем либо:

  • придерживайтесь чистой реализации обратного выборочного преобразования Python в части I (в конце концов, 2s неплохо для выборки 1,000,000 в большинстве контекстов)
  • написать нашу собственную процедуру выборки — и желательно написать эту процедуру выборки на C или Cython — что немаловажно

В следующей части мы рассмотрим именно это — реализовать эффективную пользовательскую функцию выборки распределения в инфраструктуре SciPy. Это дает нам лучшее из обоих миров — гибкость для реализации точного распределения по нашему выбору, а также использование эффективных и хорошо написанных методов, которые мы наследуем от классов rv_generic и rv_continuous SciPy.