Журнал вычислительной математики и математической физики, 2020, T. 60, № 11, стр. 1807-1814

1. ПРОБЛЕМА

Методы Монте-Карло широко используются в различных областях прикладной математики: вычисление многомерных интегралов, разыгрывание столкновений в методах молекулярной динамики и задачах реакторной защиты, оценка отказов в сложных конструкциях и многие другие. Все эти методы используют некоторые последовательности случайных чисел с заданной функцией распределения. Чаще всего такие задачи сводятся к использованию последовательности случайных чисел ${{\gamma }_{n}}$, равномерно распределенных на отрезке [0, 1]. Такие последовательности имитируют наборами псевдослучайных чисел, получаемых с помощью некоторых математических алгоритмов. Все эти алгоритмы являются по сути эвристическими.

Любой предлагаемый алгоритм обычно тестируют – проверяют, насколько выполняются те или иные свойства, характерные для чисто случайных последовательностей. Например, близко ли выборочное среднее от ${{\gamma }_{n}}$ к 1/2; существенна ли корреляции между ${{\gamma }_{n}}$ и ${{\gamma }_{{n + 1}}}$, равномерность распределения точек в кубах разной размерности и т.д. В [1] тестировалось 13 популярных генераторов псевдослучайных чисел. При этом использовалась стандартная библиотека diehard-тестов, включавшая проверку равномерности заполнения единичного куба с числом измерений до 15. Автор работы констатирует, что последовательности удовлетворяют тестам, однако, он считает, что это не гарантирует хорошего качества последовательностей. Поэтому систему формальных математических тестов нельзя считать достаточно надежной.

Математики дают путевку в жизнь только тем генераторам, которые проходят подобные тесты. Затем физики, которые не всегда хорошо знают математику, но очень часто ее обожествляют, доверчиво используют эти последовательности в своих расчетах. Рассмотрим критически эту проблему. Любой математический алгоритм является детерминированной процедурой, поэтому любая псевдослучайная последовательность не может быть истинно случайной и должна содержать какие-то закономерности. Вопрос лишь в одном – насколько они существенны. Мы протестировали ряд широко употребительных алгоритмов и убедились в том, что их вряд ли можно считать надежными. Данная работа посвящена иллюстрации этих выводов.

2. ПРОСТЕЙШИЙ ВИЗУАЛЬНЫЙ ТЕСТ

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

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

Заметим, что в этом сложном способе мы наблюдаем не всю трехмерную картину, а лишь ее проекцию на плоскость. Поэтому гораздо проще с самого начала ограничиться плоским случаем – исследованием равномерности распределения точек в единичном квадрате. Получаемую картину легко оценить визуально. При этом отметим одну существенную техническую деталь. Размер маркеров на экране или при выводе на печать надо выбирать так, чтобы в местах максимального сгущения точек они почти перекрывались бы. Для этого суммарная площадь всех маркеров должна составлять 15–25% от площади единичного квадрата (таким образом, размер маркера зависит от числа выведенных точек).

Нередко визуальный контроль объявляют субъективным. Однако известно, что глаз зачастую лучше выявляет закономерности, чем формальные математические методы. Например, глаз легко оценивает выход последовательности точек на асимптоту. Формальные математические критерии требуют гораздо большего числа точек для надежного установления такого выхода. Это хорошо известно тем, кто разрабатывает методики автоматического контроля точности расчетов.

Рассмотренные генераторы. Мы исследовали 13 генераторов псевдослучайных чисел, широко распространенных в отечественной и зарубежной литературе, а также в коммерческих пакетах:

• вихрь Мерсенна (Mersenne twister) [2];

• быстрый вихрь Мерсенна (SIMD-oriented fast Mersenne twister) [3];

• мультипликативный конгруэнтный генератор (Multiplicative congruential generator) [4];

• мультипликативный генератор Фибоначчи с запаздыванием (64-bit multiplicative lagged Fibonacci generator) [5],

• комбинированный множественный рекурсивный генератор (combined multiple recursive generator) [6];

• генератор Philo 4 × 32 [7];

• генератор Threefry 4 × 64 [7];

• генератор Марсалья (Marsaglia’s SHR3 shift-register generator) [8];

• модифицированный генератор Subtract-with-Borrow (modified Subtract-with-Borrow generator) [9];

• rand на языке C/C++;

• модифицированная последовательность Лемера [10];

• последовательность чисел из сборника [11];

• последовательности Соболя [12], [13].

Этот список частично совпадает с последовательностями, тестированными в [1]. Дадим краткие характеристики этих алгоритмов.

Вихрь Мерсенна – это весьма сложный алгоритм. В нем генерируются две различные независимые последовательности 32-битовых целых чисел, из которых составляется одно 64-битовое число. Поэтому окончательный результат соответствует точности double precision. Такая комбинация двух последовательностей обеспечивает период 2¹⁹⁹³⁷-1, далеко превосходящий все ожидания современных вычислителей. Ясно, что при такой длине невозможно использовать целую замкнутую цепочку чисел. Этот алгоритм считается одним из наиболее совершенных.

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

Мультипликативный конгруэнтный генератор – это алгоритм Лемера с некоторыми конкретными константами:

Он использует одну 32-битовую последовательность целых чисел, и его период не может превышать $M$. Это несравненно меньше, чем для вихря Мерсенна.

64-битовый мультипликативный генератор Фибоначчи с запаздыванием использует одну 64-битовую последовательность целых чисел. Его период составляет приблизительно 2¹²⁴. Это много больше, чем нужно для современных приложений. Константы запаздывания равны $l = 63$, $k = 31$.

Комбинированный множественный рекурсивный генератор использует две 32-битовые последовательности целых чисел. Рекурсивная процедура обеспечивает длину периода 2¹⁹¹.

Philox 4 × 32 – генератор с использованием четырех независимых последовательностей 32-битовых чисел. По-видимому, этот генератор ориентирован на задачи криптографии. В каждой последовательности нахождение ${{m}_{{n + 1}}}$ через ${{m}_{n}}$ осуществляется с помощью сети Фейстеля с некоторыми наборами ключей. Период составной последовательности равен 2¹⁹³.

Threefry 4 × 64 – это адаптация блочного криптографического алгоритма Threefish из хэш-функции Skein. Она использует четыре 64-битовых последовательности и обеспечивает период 2⁵¹⁴.

Генератор Марсалья использует линейный конгруэнтный алгоритм (1) с параметрами $A = 69069$, $M = {{2}^{{32}}}$, $C = 1234567$. Вычисления производятся с 32-битовыми числами. Каждое полученное по формуле (1) число подвергается следующей процедуре. Сначала оно сдвигается влево на 5 регистров и суммируется с исходным числом. Затем полученное число сдвигается на 17 регистров вправо и оба этих числа суммируются. Новое число сдвигается на 13 регистров влево и также выполняется суммирование. Все указанные суммирования выполнялись с 64 разрядами. Период полученной последовательности оценивался как 2⁶⁴. Использование в качестве $M = {{2}^{{32}}}$ степени числа 2 удешевляет вычисления, но такой выбор может ухудшить качество полученных чисел.

Модифицированный генератор Subtract-with-Borrow аналогичен генератору Фибоначчи с задержкой (константы задержки 27 и 12). Модификация, описанная в [10], позволила увеличить период до 2¹⁴⁹².

Генератор rand языка C/C++ использует алгоритм Лемера и 16-битовые целые числа. Период такой последовательности не может превышать 2¹⁶, т.е. он очень невелик и годится лишь для учебных расчетов. Несмотря на это, генератор очень популярен (быть может, потому, что он является стандартным генератором в C/C++).

Модифицированный алгоритм Лемера использует три независимые последовательности 16-битовых целых чисел. Он был распространен около 1990 г., когда широко использовались 16-разрядные PC XT. Период такого генератора не может превышать 2⁴⁸.

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

В 1955 г. американской корпорацией RAND была опубликована книга [11], содержащая миллион случайных цифр и сто тысяч стандартных нормальных отклонений. Эти случайные цифры были получены с помощью электронной рулетки, т.е. физического прибора, а не математического алгоритма.

Особняком стоят последовательности Соболя [12], [13]. Это название широко распространено в зарубежной литературе, а сам автор назвал их ЛП_τ-последовательностями. Эти числа следует называть не псевдослучайными, а квазислучайными. Они строятся так, что при магических числах точек $N = {{2}^{k}}$ проекции всех точек на любую ось единичного куба различаются и образуют равномерную сетку с шагом ${{2}^{{ - k}}} = {1 \mathord{\left/ {\vphantom {1 N}} \right. \kern-0em} N}$. При промежуточных значениях $N$ часть узлов этой сетки остается незанятой. Период такой последовательности неограничен. В [13] приведены таблицы для построения этих последовательностей при размерности пространства $p \leqslant 13$; но в настоящее время такие таблицы разработаны для $p\sim 4000$.

3. РЕЗУЛЬТАТЫ ТЕСТИРОВАНИЯ

Зависимость от N. Естественно ожидать, что при увеличении N распределение точек будет все более равномерным. Подробные расчеты были проведены для генератора Мерсенна для $N = 100,\;300,\;1000,\;3000,\;10\,000$; при дальнейшем увеличении N точки сливались, и визуальное исследование становилось проблематичным. Начало последовательности было одинаковым во всех этих расчетах. Типичные примеры для $N = 1000$ и $N = 3000$ приведены на фиг. 1 и 2 соответственно. На фиг. 1 диаметр точек в $\sqrt 3 $ раз больше, чем на фиг. 2, чтобы суммарная площадь точек была одинаковой.

На обоих рисунках хорошо видны области сильного сгущения и сильного разрежения точек. Важно то, что фиг. 2 строится по той же последовательности с тем же началом, т.е. он получен наложением на фиг. 1 следующего отрезка последовательности длиной 2000. Тем не менее области сгущения и разрежения точек оказываются в других участках квадрата. При дальнейшем увеличении числа точек места сгущения и разрежения также смещаются. Но даже при очень больших $N = 10\,000$ неравномерность заполнения остается значительной. Это свидетельствует о существенной неравномерности заполнения квадрата и о не слишком высоком качестве генератора.

Начало последовательности. На фиг. 3 приведен расчет для генератора Мерсенна с $N = 3000$, но со случайно выбранным началом последовательности. Сравнение с фиг. 2 показывает, что неравномерность заполнения, по-прежнему, достаточно велика, а области наибольшей и наименьшей плотности занимают другие места. Аналогичная картина была при других начальных точках последовательности. Таким образом, все участки последовательности характеризуются неравномерностью распределения точек.

Другие генераторы. На фиг. 4–14 приведены результаты тестирования других 11 генераторов псевдослучайных чисел. Для всех выбрано $N = 3000$. Видно, что картина оказывается такой же, как для генератора Мерсенна. При этом отметим, что одна из этих последовательностей RAND (см. фиг. 14) использует не математические псевдослучайные числа, а электронную рулетку – электронный источник шумов, от которого следовало бы ожидать истинной случайности. Тем не менее конструктивные особенности датчика внесли, по-видимому, неслучайность в результаты.

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

Последовательность Соболя. Результаты для нее при магическом числе $N = 2048$ представлены на фиг. 15. Наблюдаемая картина качественно отличается от фиг. 1–14. Видно, что заполнение имеет характер регулярного кружева. Точки в среднем распределены в квадрате гораздо равномернее, чем для псевдослучайных последовательностей, хотя на малых участках квадрата размерами $\sim {{N}^{{ - 1/2}}}$ их распределение неравномерно.

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

Сетка точек при $N \ne {{2}^{k}}$ получается из магической сетки выбрасыванием части точек. Интуитивно следует ожидать, что результаты вычисления многомерного интеграла при немагическом числе точек окажутся хуже, чем при магическом.

4. ЗАКЛЮЧЕНИЕ

В настоящее время вряд ли существует какой-то алгоритм построения псевдослучайных чисел, который можно с полной уверенностью использовать для расчета любых случайных процессов. В литературе встречались высказывания о том, что хорошие случайные свойства проявляютcя только на отрезках очень большой длины $N\sim {{10}^{9}}$. Вероятно, при меньших $N$ известные алгоритмы можно применять для тех или иных конкретных задач, но такие ситуации нужно дополнительно обосновывать. Например, для многомерных кубатур можно уверенно рекомендовать последовательности Соболя с магическими числами точек.

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

Авторы искренне благодарны И.М. Соболю, И.А. Козлитину, Д.Д. Соколову за полезные обсуждения.