Простое объяснение современных генераторов псевдослучайных чисел (PRNG) и их новой реализации NumPy.

Создавая образец нормально распределенных чисел, мне было любопытно, откуда они взялись — в частности, как компьютер может создавать числа, которые следуют выбранному распределению, будь то нормальное, экспоненциальное или что-то более необычное. Какой бы метод ни лежал в основе создания этих нормально распределенных значений (инверсионная выборка, Бокс-Мюллер или быстрый алгоритм Зиккурата), все они начинаются с одного фундаментального строительного блока: последовательности равномерно распределенных значений.

Тогда возникает вопрос: откуда они берутся? В большинстве случаев: от генератора псевдослучайных чисел — или PRNG. Заметив, что NumPy изменил свой PRNG по умолчанию еще в июле 2019 года в ответ на этот NEP (несмотря на то, что в Интернете все еще полно старых способов сделать это), я очень хотел понять:

  • почему они изменили это
  • как работают PRNG в первую очередь

а затем запишите мои мысли на простом английском языке с некоторыми полезными диаграммами и примерами кода. Это заметки, в значительной степени основанные на этой (довольно доступной) статье и этой связанной с ней лекции Мелиссы О’Нил — автора семейства ГПСЧ PCG еще в 2014 году, которые теперь формируют ГПСЧ по умолчанию в NumPy.

Почему «псевдо»?

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

Существуют способы убедиться, что эти входные данные являются «действительно» случайными, которые в основном включают аппаратное обеспечение для измерения действительно случайных явлений природы (например, атмосферного шума), но обычно мы «задаем» алгоритмы выбранным (неслучайным) начальным значением. Даже если мы используем машинное время (как целое число), когда генератор запускается как начальное значение, это не является 100% случайным. Однако мы, вероятно, должны задаться вопросом, стоит ли вообще стремиться к истинной случайности.

Хотим ли мы «настоящей случайности»?

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

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

  • не знаю начальное значение
  • знать «достаточно» последовательности невозможно на практике

Что тогда вызывает следующий вопрос:

Какие свойства PRNG желательны?

Скорость, размер и воспроизводимость

Здесь нет ничего особенного — нам нужно что-то, что могло бы быстро генерировать воспроизводимое большое количество случайных чисел и не занимать слишком много памяти, особенно если мы планируем запускать многие из них в разных потоках.

Единообразие

Если мы знаем, что сгенерированные случайные числа попадут в интервал [0, n], то мы хотели бы, чтобы каждое число в этом интервале имело равные шансы быть выбранным — в противном случае, хотя теоретически у нас много чисел для выбора на практике это будет намного меньше. Например, если у нас есть возможные числа в интервале [0, 255], но на практике наш PRNG выбирает только 3 или 250, то это совсем не случайно. С другой стороны, нам не обязательно нужен алгоритм, гарантирующий единообразие. Если это так, то если у нас есть диапазон n чисел для выбора, и мы уже выбрали n-1 чисел, то окончательное число вовсе не случайное — это будет просто число, которое еще не было выбрано. Работая в обратном направлении, мы можем увидеть, что любой алгоритм, гарантирующий единообразие, будет уменьшать случайность по мере того, как мы достигаем конца диапазона, также известного как «период».

Долгий период

Если мы согласимся с тем, что для того, чтобы заставить компьютер генерировать последовательность, мы должны дать ему функцию, которая преобразует предыдущее число в следующее число, то, если мы когда-либо получим число, которое мы видели раньше, мы начнем повторять последовательность. . Это гарантированно произойдет в какой-то момент, пока мы выбираем числа из ограниченного интервала, потому что, если мы выбираем числа в [0, n], то по определению n+1 число должно быть повторением предыдущего числа. Другими словами, в последовательности подходящего размера все ГПСЧ будут дублироваться и, таким образом, станут детерминированными. Чтобы предотвратить это, мы можем сделать количество чисел («период») до того, как это произойдет, намного больше, чем длина последовательности чисел, которую мы хотим выбрать.

Отсутствие шаблона

Кажется очевидным, но стоит добавить. Например, хороший алгоритм для удовлетворения всех вышеперечисленных свойств — просто «добавить 1»:

которые будут:

  • молниеносно быстрый, маленький и воспроизводимый (нам просто нужно знать X_0 значение, с которого мы начали)
  • униформа: каждое число в заданном интервале будет выбрано один раз
  • длительный период: он никогда не будет повторяться до тех пор, пока мы никогда не «обернем его вокруг» на каком-то значении, чтобы убедиться, что размер числа не станет большим, т.е. сопоставим большое число, y, обратно к нулю и начнем снова

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

Этот бит является ключевым — все PRNG будут детерминированными, если вы знаете «состояние» алгоритма (например, последнее случайное число и функцию для преобразования его в следующее случайное число). Суть в том, что без этой информации числа будут казаться случайными — точно так же, как индейка думает, что Рождество — это «Черный лебедь», а фермер — нет (случайность зависит от вашего набора информации).

LCG

Линейные конгруэнтные генераторы (LCG) — один из старейших PRNG, и, к счастью, довольно простой для понимания.

Другими словами — чтобы получить следующее число в последовательности, мы берем предыдущее число и:

  • умножьте на некоторую константу a
  • добавьте к нему другую константу c
  • взять остаток при делении на какую-то другую константу m

Пока ничего слишком «информатического» со страшными словами или фразами, такими как «матричная линейная рекуррентность» или «простое число Мерсенна». Давайте выберем некоторые значения для {a, c, m} и посмотрим, как выглядит результат.

В частности, давайте сформируем генератор в стиле генератора Лемера 1951 года, который предельно прост — a=3, c=0, и чтобы гарантировать, что мы генерируем 8-битные числа, установим m=2^8=256. Этот последний бит просто означает, что числа останутся между [0, 255] и, таким образом, будут соответствовать 8 битам.

3, 9, 27, 81, 243, 217, 139, 161, 227, 169, 251, 241, 211, 121, 107, 65, 195, 73, 219, 145, 179, 25, 75, 225, 163, 233, 187, 49, 147, 185, 43, 129, 131, 137, 155, 209, 115, 89, 11, 33, 99, 41, 123, 113, 83, 249, 235, 193, 67, 201, 91, 17, 51, 153, 203, 97, 35, 105, 59, 177, 19, 57, 171, 1, 3, 9, 27, 81, 243, 217, 139, 161, 227, 169, 251, 241, 211, 121, 107, 65, 195, 73, 219, 145, 179, 25, 75, 225, 163, 233, 187, 49, 147, 185, 43, 129, 131, 137, 155, 209

Итак, какие проблемы мы уже видим?

  • все числа нечетны, поэтому мы вообще не касаемся всех чисел одинаково (любая наблюдаемая закономерность является антитезой случайности)
  • период короткий — вы можете увидеть примерно половину пути, мы возвращаемся к началу, а затем начинаем повторяться

Если LCG настолько ужасны, то актуально ли это для PRNG сегодня?

Потому что LCG не ужасны. Приведенная выше LCG ужасна, но это связано не столько с LCG как семейством самих генераторов чисел, сколько с тем, как мы параметризовали этот генератор. Фактически, следующая LCG называется drand48 и лежит в основе java.util.Random:

но имеет принципиальное отличие от приведенной выше спецификации LCG.

Выводит не просто «состояние», а функцию от него

В первом примере мы просто сгенерировали следующее число в последовательности и вывели его. Если это было 255, то следующим числом в нашей последовательности было 255. Не смешное дело. В приведенной выше реализации LCG это не так — мы имеем следующий yield seed >> 16. В python это побитовый оператор, который сдвигает все биты на 16 позиций вправо, в результате чего самые правые 16 бит удаляются.

Здесь мы можем взять пример — если у нас есть число 1017, мы можем представить его в двоичном виде как 11 1111 1001 (промежуток просто для удобства чтения) — если мы делаем 1017 >> 3, то мы получаем 111 1111 (то есть 127), т. е. мы сдвигаем все к справа на 3 места и отбросить первые 3 бита справа. Это всего лишь один из таких типов функций, иллюстрирующий способ улучшения результатов. Теперь у нас есть следующая настройка для нашего алгоритма LCG:

  • генерировать следующее число в последовательности — это «состояние» LCG (просто терминология)
  • используйте это «состояние» для генерации «вывода» — это число, которое затем используется для формирования части последовательности

Вот как мы можем иметь «16-битный генератор» с «8-битным выходом» — потому что «состояние» генератора — это 16-битное число, а выход — 8-битное число, где этот 8-битный число создается путем применения некоторой функции к 16-битному состоянию. Как мы увидим, создание PRNG с другим выходом в состояние может значительно улучшить их статистические свойства.

Рандограммы: визуализация случайности

Чтобы интуитивно понять, насколько случайны различные алгоритмы, нам нужен способ визуализировать эту случайность. Для этого снова позаимствуем идею из статьи Мелиссы О’Нил. На практике генераторы случайных чисел намного больше, чем мы будем делать (64-битное состояние с 32-битным выходом), но идея заключается в следующем:

  • создать генератор с 16-битным состоянием, т. е. начальное значение/состояние находится в диапазоне[0, 65535] (где верхняя граница 2**16-1)
  • вывести 8-битное число, полученное из этого 16-битного состояния, т.е. каждое число в выходной последовательности будет в [0, 255]
  • возьмите последовательность и сгруппируйте соседние точки в пары, т.е. [x_0, x_1], [x_2, x_3], etc
  • эти пары образуют {x, y} координаты на графике 256 x 256
  • используйте PRNG для генерации координат 2^16 и нанесите их на график, где, если у нас нет пар для координаты, тогда эта координата белая, и чем больше пар у нас есть для данной координаты, тем более черной будет эта координата. быть

В некотором смысле это даст нам красивую картину случайности. Если у нас есть хороший алгоритм, у нас будет график с множеством точек, который в целом выглядит случайным. Это будет намного легче увидеть на примере. Давайте возьмем хорошо параметризованное «16-битное состояние, 8-битный выход LCG» и нарисуем несколько «рандограмм».

Верхняя левая диаграмма показывает, что мы получаем, когда берем крайние правые 8 битов 16-битного номера состояния. Например, если нашим «состоянием» для итерации является число 65413, представленное в двоичном виде как 1111 1111 1000 0101, то мы выведем крайние правые 8 бит — 1000 0101 (или 133). Мы делаем это для всех чисел, группируем их попарно и наносим на график.

Мы видим, что числа совсем не кажутся случайными — они образуют аккуратные прямые линии. Это теорема Марсалья, которая показывает проблему с LCG, когда у нас слишком маленький период (и поэтому получить повторяющиеся значения, подобные этому). Однако по мере того, как мы продвигаемся выше, 16-битное состояние начинает выглядеть немного лучше. На нижней правой диаграмме все еще есть четкая структура, но мы справляемся намного лучше с точки зрения охвата пространства.

Поэтому, глядя на это, мы можем сделать следующее наблюдение: чем выше группа из 8 битов в 16-битном номере состояния, сгенерированном LCG, тем более случайными они кажутся.

Введите PCG

Несмотря на то, что LCG все еще находят широкое практическое применение, они не были PRNG по умолчанию для NumPy до 2019 года. Вместо этого до NumPy 1.17 использовался алгоритм под названием Вихрь Мерсенна — в частности, MT19937 — названный из-за его (абсолютно колоссальной) длины периода, представляющей собой простое число Мерсенна (2**19937 - 1 — степень 2 минус 1). Однако с выпуском NumPy 1.17 он переключился на PRNG по умолчанию, являющийся PCG - Permuted (Linear) Congruential Generator.

PCG — это семейство генераторов, созданных Мелиссой О’Нил, которые умело используют наблюдения, сделанные выше, особенно диаграммы. Идея такова:

  • вывод функции состояния, а не непосредственно состояния, по-видимому, увеличивает случайность
  • В LCG явно отсутствует случайность в младших битах (верхний левый график), но старшие биты имеют тенденцию быть «более случайными» (нижний правый график).
  • если у нас есть, например. 16-битное состояние выводит 8-битное число, тогда нам нужно выбрать только 8 бит для вывода
  • почему бы нам не использовать несколько старших битов 16-битного состояния, которые являются наиболее случайными, чтобы выбрать, какую функцию мы применим к оставшейся части 16-битного состояния для генерации 8-битного вывода?
  • другими словами, давайте используем наиболее случайную часть нашего состояния, чтобы случайным образом выбрать функцию преобразования для применения к остальной части состояния — своего рода рандомизированный алгоритм.

Давайте посмотрим на простой пример.

ПКГ РС

Приведенные выше 9 диаграмм делают следующее: вычисляют 8-битный вывод из 16-битного состояния, где эти 8-битные выходные данные генерируются заранее определенным битовым сдвигом (от 0 для верхнего левого до 8 для нижнего правого ). Но как насчет фиксированного смещения вместо случайного смещения?

Откуда у нас эта случайность? Из нескольких верхних битов нашего 16-битного состояния. Другими словами, если у нас есть состояние 57277 (1101 1111 1011 1101), мы могли бы:

  • используйте верхние 2 бита, 11, для определения сдвига - в данном случае 3
  • применить этот сдвиг к другим битам, 01 1111 1011 1101 s.t. вместо выбора крайних левых 8 бит, 0111 1110, мы сдвигаем 3 вправо и выбираем биты 1111 0111

В некотором смысле то, что мы делаем, глядя на приведенные выше 9 диаграмм, использует случайность графиков 8–15, 9–16, чтобы случайным образом выбрать, выбираем ли мы число из графиков {4-11} - {7-14}. Давайте посмотрим, как это выглядит:

Очевидно, что там все еще есть какая-то структура, но улучшение огромно, если просто взять 2 старших бита 16-битного состояния и использовать их для выбора перестановки для остальной части состояния — отсюда и название Permuted (Linear) Congruential Generators. Но это всего лишь одно из таких преобразований — битовый сдвиг.

А как насчет других трансформаций? Существует целый ряд доступных преобразований (исключающее ИЛИ, вращение и т. д.), которые создают семейство PCG, где старшие биты случайным образом выбирают, какую перестановку применить к линейно сгенерированному состоянию. Давайте посмотрим на 2 других таких перестановки, которые можно использовать.

Вращение

Одно из таких преобразований, которое мы можем (произвольно) выбрать, — это «вращение» битов. Так же, как раньше с помощью оператора >> мы сдвигали биты вправо и отбрасывали переполнение, при вращении мы сдвигаем одно направление, но тогда вместо отбрасывания переполнения мы оборачиваем его в другую сторону.

Если у нас есть число 113, представленное в двоичном виде как 111 0001, мы можем выполнить «правое вращение» 2, чтобы создать 011 1100 или 60. Здесь мы взяли 2 самых правых бита (01) и повернули их к самому началу и сместили все вниз.

xorshift

Еще одно преобразование, которое мы могли бы применить, — это «xorshift». Опять же, давайте проиллюстрируем тот же пример. Снова взяв 113 (111 0001), мы могли бы:

  • «сдвинуть» его вниз на величину, например. 2, чтобы получить 001 1100
  • применить функцию побитовое исключающее ИЛИ к исходному номеру и смещенному номеру

В этом случае мы должны вычислить xor (или ^ в Python) на 111 0001 и 001 1100, чтобы получить 110 1101 (1 с, когда только 1 бит равен 1, 0, если оба равны 1 или 0).

PCG XSH-RS и PCG XSH-RR

Теперь давайте посмотрим на рандограммы для двух распространенных ГКП и посмотрим, как они соотносятся с приведенными выше. Они есть:

  • PCG XSH-RS: сначала вычислите операцию xorshift, затем случайным образом сдвиньте полученные биты
  • PCG XSH-RR: сначала вычислить операцию xorshift, а затем случайным образом повернуть полученные биты

Опять же, структура все еще есть, но они заметно улучшились. Эта структура существует, потому что мы используем маленькие генераторы. Точно так же, как старшие 8 бит более случайны, чем нижние 8 бит в 16-битном состоянии, верхние 16 бит более случайны, чем нижние в 32-битном состоянии. Кроме того, поскольку мы используем все большее и большее состояние, мы естественным образом увеличиваем период. Сочетание этих двух факторов позволяет даже очень большим (96-битное состояние, 32-битный вывод) LCG пройти тест Big Crush — обширный набор статистических тестов, составленный Пьером Л'Экюйером и Ричард Симард, чтобы проверить PRNG на наличие желаемых свойств, упомянутых ранее.

Учитывая дополнительную случайность случайно выбранных перестановок, PCG работают намного лучше, чем LCG, и в результате им не нужны такие большие состояния для прохождения набора тестов. Именно по этой причине они были приняты в качестве PRNG по умолчанию в NumPy — с точным алгоритмом PCG XSL RR 128/64.