О рандомизации квазислучайных последовательностей Холтона Текст научной статьи по специальности «Математика»

УДК 519.245 Вестник СПбГУ. Математика. Механика. Астрономия. 2017. Т. 4 (62). Вып. 4

MSC 65C05

О РАНДОМИЗАЦИИ КВАЗИСЛУЧАЙНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ ХОЛТОНА*

С. М. Ермаков

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

Рассматривается вопрос об оценке погрешности методов квази Монте-Карло с помощью рандомизации. Известное неравенство Коксмы—Хлавки позволяет судить об асимптотике погрешности, но совсем не пригодно для практического использования в процессе вычислений, так как вычисление входящих в него величин — вариации функции и дискрепанса последовательности чрезвычайно трудоемки и практически неосуществимы. По этой причине имеются многочисленные попытки использовать средства теории вероятностей для решения указанной задачи. Одним из распространенных подходов является случайный сдвиг точек псевдослучайной последовательности. Известны случаи практического использования этого подхода, но теоретически он мало исследован.

В данной работе показано, что полученные таким образом оценки являются оценками сверху, установлена связь с теорий кубатурных формул с одним случайным узлом. Подробно рассмотрен случай последовательностей Холтона. Анализируется преобразование Ван дер Корпута последовательности натуральных чисел, с помощью которого строятся точки Холтона. Показано, что кубатурная формула с одним свободным узлом, соответствующая последовательности Холтона, точна для некоторого класса ступенчатых функций. Класс явно описан. Полученные результаты позволят более эффективно использовать указанные последовательности при вычислении интегралов и поиске экстремума, а также могут служить отправной точкой для дальнейших теоретических исследований в области квазислучайных методов. Библиогр. 6 назв. Ил. 1. Табл. 1.

Ключевые слова: квазислучайные последовательности, метод Монте-Карло, последовательности Холтона, случайные квадратурные формулы.

Введение. Работа посвящена методам численного интегрирования «гибридного» типа, использующим методы Монте-Карло и квази Монте-Карло в совокупности. Известно, что одним из важных преимуществ метода Монте-Карло при вычислении интегралов по вероятностной мере является возможность оценивания погрешности в процессе этих вычислений. Метод квази Монте-Карло, являющийся детерминированным методом, этой возможности лишен. В связи с этим возникла необходимость создания «гибридных» методов, использующих рандомизованные последовательности с малым отклонением (discrepancy). Эти методы уже допускают оценку погрешности с помощью доверительных интервалов, но могут быть несколько менее точными, чем методы квази Монте-Карло в чистом виде. Оказалось, что, по крайней мере, некоторые из методов рандомизации имеют тесную связь со случайными кубатурны-ми формулами, теория которых достаточно хорошо развита (см., например, [1, 2]). В частности, алгоритмы, возникающие при упомянутой рандомизации, тесно связаны со случайными кубатурами с одним свободным узлом и являются их обобщениями. В данной работе эта связь подробно прослеживается на примере квазислучайных последовательностей Холтона [1], которые наиболее просты по своей конструкции из известных последовательностей. Указаны функции, для которых построенные формулы точны, и получено неравенство, которому удовлетворяет дисперсия кубатурной

* Работа выполнена при финансовой поддержке РФФИ (грант № 17-01-00267-а).

© Санкт-Петербургский государственный университет, 2017

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

1. Постановка задачи. Рассмотрим последовательность векторов X^X^..., Xn в единичном s-мерном гиперкубе Ds, Ds = {0 ^ xi ^ 1; i = 1,..., s}, Xi G Ds, l = 1,..., N. Отклонением (star discrepancy) последовательности

X,

.(i) i,

, x(i) j xs

называют величину

D*n=DN(Xi,..., Xn)

= sup N

X£DS

N

j=1i=1

x(i)-x(i)

8(Z )=

1, Z<0; 1/2, Z=0; 0, Z>0,

(1)

характеризующую отклонение распределения последовательности от равномерного.

Известно, что наилучший порядок убывания D*N есть O(lnsN) [1]. Последовательности, для которых такой порядок имеет место, называют квазислучайными или последовательностями с малым отклонением (low discrepancy). Величина D*N тесно связана с оценкой остатка приближенного интегрирования с помощью кубатурной суммы специального вида:

1 N

/ f(X)dX--YJf(Xl)

J А— 1

< Var/—^

(Xb...,XN)

N

(2)

(неравенство Коксмы—Хлавки [1]), где Var есть вариация функции f в смысле Харди—Краузе.

Таким образом, убывание остатка интегрирования с помощью квазислучайных последовательностей (метод квази Монте-Карло) есть O(lnsN/N), что выгодно отличается от случайных последовательностей, для которых остаток убывает как O(N-1/2). Имеется в виду случай, когда Xi —независимые равномерно распределенные в Ds векторы.

Не погружаясь в детали сравнения случайных и квазислучайных методов (по этому поводу см., например, [3]), отметим только обстоятельство, что формула (2) мало пригодна для реальной оценки погрешности в процессе вычислений, тогда как средства математической статистики легко позволяют строить доверительный интервал для оценки интеграла. Достаточно предполагать интегрируемость с квадратом функции f.

Это обстоятельство привело к появлению различных приемов рандомизации квазислучайных последовательностей (scrambling) [4, 5]), один из которых (наиболее употребительный) состоит в следующем. К каждому вектору Xi добавляется равномерно распределенный в Ds случайный вектор. Если при этом некоторая компонента результата оказывается больше единицы, то берется ее дробная доля. Таким образом, для каждого l = 1, . . . , N будем иметь

Xi = {Xi + Si},

(3)

где X; — рандомизованный вектор, а 2; — случайный вектор, равномерно распределенный в Юя.

2. Оценка погрешности. Если векторы 2; независимы, то нетрудно показать, что / (X;) также независимы и одинаково распределены, а рандомизованная оценка ничем не отличается от оценки, полученной с помощью обычного метода Монте-

Карло. Содержательный результат может быть получен, если исходная последовательность векторов будет разбита на М групп (для простоты одинаковой численности, N = Мт) так, чтобы в каждой группе векторы были одними и теми же, но разными (независимыми) между группами. Имеем в к-й группе векторы Х(к-1)т+1, к = 1,..., М, I = 1,... ,т, и рандомизованные векторы Х(к-1)т+; = {Х(к-1)т+; + 2к }. В каждой группе можно построить несмещенную оценку интеграла

1 т

Кк,т(Л = ~ Е $ {Х{к-1 )т+г) • (4)

1 = 1

Практически для оценки дисперсии случайной величины Кк,т(/) достаточно М^10 испытаний. Легко видеть, что сумма (4) является случайной квадратурной суммой с одним свободным узлом. Теория кубатур с одним свободным узлом достаточно разработана [2]. Своеобразие в данном случае состоит в том, что величины кк,т (/) различаются при разном к не только реализациями случайного свободного узла, но и входящими в них квазислучайными векторами X;. Имеем

/ 1 м \ 1 м

^ м Е = Ж Е ПхьМ = **(/)> Х=Мт, (5)

V к=1 ) к= 1

что и позволяет построить доверительный интервал для среднего при заданном уровне доверия р:

у

2 ' е"2^2 ¿и.

1 N

у дх) ¿х--][>„,т(/)

Юа к=1

< щШ. (в)

л/Ж

Легко понять, что можно также оценивать дисперсию отдельной случайной величины Кк,т при фиксированном к, но привлекая N векторов Х^, г = 1,..., М, можно ожидать большей точности оценки. Что касается оценки детерминированной погрешности, то неравенство (6) служит, очевидно, лишь ориентировочной ее оценкой — оценкой сверху.

3. Последовательности Холтона. Далее подробно рассмотрим частный случай кубатурной формулы

У"/(X) ¿X « Кк,т(/) (7)

для последовательности квазислучайных векторов Холтона [6] — одной из простейших последовательностей такого рода.

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

акгк 1, то рг (п) = ^^ акг к (га&а1 шуегэе Гипсйоп). (8)

"к 1, то (п) = ^2 а г к к=1 к=1

Обсудим некоторые свойства этого преобразования.

р

С.1. Если п принимает значения 0,1, 2,.. . ,гй-1, то ^у (п) принимает все значения из множества п/гй. Доказательство тривиально, если представить п в гй-ичной системе счисления. Имеем один разряд, который и принимает значения 0,... , гй-1.

С.2. Если Я2,..., ат в (8) фиксированы, а Я1=/ — любое целое, 0 ^ I ^ г — 1, то (п) € [(I — 1)/г, 1/г) = Д;. Действительно, при преобразовании ^ последний разряд в представлении становится первым.

Рассмотрим теперь поведение ^у (п) при случайном сдвиге. Легко проверить следующее утверждение.

С.3. Если у1,..., уг —это г равноотстоящих точек (Ж = г), у; = у1 + (I — 1)/г на отрезке [0,1], у1 € [0,1/г), то точки у; = {у; + а} также равноотстоящие, и каждая из них принадлежит одному и только одному из отрезков Д;, I = 1, 2,..., г. При г = г^, целом к > 1 получаем, что если точка у; € Д(й), Д(й) = ((I — 1)/гх, I/г ^), то она также

А 1)

очевидно принадлежит и Д; .

При г^ < г < 1 часть отрезков Д(й 1 ) не будет содержать ни одной из г точек. При этом расстояние между точками не может быть меньше 1/г^.

Вспомним теперь алгоритм получения векторов последовательности Холтона [1, 6]. Выбираются в взаимно простых чисел г 1,г2,... ,гя, и для й = 0,1,... последовательность вычисляется по формулам

Хй = (^ И,...,^И), й = 0,1,..., (9)

а рандомизованные векторы — по формулам

Х] = ({^1 (й) + а(])},..., {^ И + а^}) , ¿=1,..., к. (10)

Рассмотрим сначала одномерный случай. Как и ранее будем последовательно выбирать группы по т точек Хй и строить ).

Кубатурная формула

У"/(X) йХ « (/), (11)

как уже отмечалось, является обобщением достаточно полно исследованных [2] формул с одним свободным узлом. Для этих формул один из узлов (первый узел X 1) является случайным вектором с заданным распределением, а узел с номером й определяется как Хй = Тй(Х1), где Т — некоторое отображение, заданное на множестве векторов размерности в. Иначе, Хй = Т(Х^_1). Отличие этих формул от введенных нами в рассмотрение формул вида (10) состоит в том, что отображение, с помощью которого вычисляется Хй по Х^_1, зависит от номера й: Хй = Тй(Х^_1) или Хй = П^=2 Т((Х1). Множество отображений Т зависит также, вообще говоря, от номера группы к.

Заметим, что при с! = 0,..., г1 — 1 точки Хй будут принадлежать интервалам Дй = ((й — 1)г, й/г] последовательно (свойство С.1). Далее, при й = 0,..., сг1 — 1, где с — заданное натуральное число, в каждый отрезок Дй попадут с точек последовательности.

При этом, если с = С1 г], с, ] —натуральные числа, то С1 точек будет в каждом

(^+1) = (а_ 1)/г]+1 л/г^1 й

из отрезков Д]+1) = ^(й — 1)/г]+1, й/г когда с кратно некоторой степени г1 .

. Далее мы рассматриваем лишь случай,

В этом случае имеем

1 сг 1-1

Хк,т = — > , сг1 1 ^=0

и соответствующая формула точна для всех ступенчатых функций, постоянных на

Обобщение на многомерный случай осуществляется достаточно просто, когда число точек т кратно каждому из т^, у = 1,..., в. Рассмотрим простейший пример: в = 2, Т1 = 2, Т2 = 5, т = 10. Координаты 10 точек приведены в следующей таблице.

п 0 1 2 3 4 5 6 7 8 9

411 0 1/2 1/4 3/4 1/8 5/8 3/8 7/8 1/16 9/16

0 1/5 2/5 3/5 4/5 1/25 6/25 11/25 16/25 21/25

Расположение точек схематически отражено на рисунке.

Каждая точка принадлежит параллелепипеду со сторонами длины 1/5 х 1/2 и случайный сдвиг с циклическим переносом сохраняет это свойство (свойство С.3).

-2 • 4 •10

.1 • 7 '3 -9 "5

1/5 1/5 3/5 4/5 1

Из приведенных рассуждений очевидным образом следует теорема. Теорема 1. Кубатурная формула с одним свободным узлом (11), полученная в результате рандомизации отрезка длины т последовательности Холтона (9) при т = П^=1 Щ, точна по крайней мере для ступенчатых функций, постоянных на

гиперпараллелепипедах вида Д^...^ = < х^ = 1,...,в}, где Щ = г1-'

г = 1,..., в, рг — натуральное число.

Следствие. Дисперсия кубатурной суммы (4) при т = П неравенству

Щ удовлетворяет

к,т

<

1

|/2(х) ¿х - ¿Ц /(X)

1 = 1 \Аг

¿X

где I = (я,...,я) — мультииндекс я = 1,..., Д».

Последнее утверждение вытекает из неравенства Минковского и свойств расслоенной выборки.

Полученные результаты позволяют качественно сравнить поведение рандомизо-ванных последовательностей Холтона и Соболя. Последние, как правило, связаны с квадратурными формулами, точными для функций Хаара. Можно надеяться также, что эти результаты помогут связать для последовательностей Холтона выбор чисел т со свойствами подынтегральной функции.

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

2

т

последовательностей. Поскольку кубатурные суммы ) при различных к зави-

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

Если Кп[/] =5^; А; (У1,..., УП)/(У) является кубатурной суммой со случайными узлами, то для оценки интеграла можно рассматривать сумму

где величина кт,;(/) использует в качестве случайного параметра сдвига вектор Y;. Простейшим примером такого подхода может быть применение метода противоположной переменной (antithetic variate), в котором используется случайная кубатур-ная формула

где I — вектор, все компоненты которого равны единице.

Другие удобные для применения случайные кубатурные формулы можно найти в монографии [2].

Литература

1. Кейперс Л., Нидеррейтер Г. Равномерное распределение последовательностей. М.: Наука, 1985. 408 с.

2. Ермаков С. М. Метод Монте-Карло и смежные вопросы. 2-е изд. М.: Наука, 1975. 472 с.

3. Ермаков С. М. Метод Монте-Карло в вычислительной математике. Вводный курс. СПб.: Невский диалект. М.: Бином. Лаборатория знаний, 2009. 192 с.

4. L'Ecuyer P., Lecot C., Tuffin B. A Randomized Quasi-Monte Carlo Simulation Method for Markov Chains // Operations Research. 2008. Vol. 56, N4. P. 958-975. https://doi.org/10.1287/opre.1080.0556

5. Chi H. Scrambled Quasi-Random Sequences and their Application. A dissertation for the degree of Doctor of Philosophy. The Florida state university college of arts and science, 2004. URL: http://diginole.lib.fsu.edu/islandora/object/fsu:182066/datastream/PDF/view (дата обращения: 05.07.2017).

6. Niederreiter H. Quasi-Monte Carlo methods and pseudo-random numbers // Bull. Amer. Math. Soc. 1978. Vol. 84, N6. P. 957-1041. https://doi.org/10.1090/S0002-9904-1978-14532-7

Статья поступила в редакцию 16 мая 2017 г.; рекомендована в печать 22 июня 2017 г. Сведения об авторе

Ермаков Сергей Михайлович —доктор физико-математических наук, профессор; [email protected]

ON THE HALTON QUASI-RANDOM SEQUENCES RANDOMIZATION

Sergey M. Ermakov

St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation; [email protected]

The problem of estimating the error of quasi-Monte Carlo methods with the use of randomization is considered. The well-known Koksma—Hlawka inequality allows us to estimate the asymptotics for the error, but it is not suitable for practical use in the process of computation, since calculation of the

n

]ГАг(У1,...,У„)кт,г (/),

г=1

Ds

quantities occurring in it is a variation of the function and discrepancy of the sequence, is an extremely time-consuming and impractical process. For this reason, there were numerous attempts to solve the aforementioned problem with the probability theory methods. One common approach is to randomly shift the points of the pseudo-random sequence. Cases of practical use of this approach are known, but theoretically it has been scantily studied.

In this paper it is shown that estimates obtained this way are the upper estimates. A connection with the theory of cubature formulas with one random node is established. The case of the Halton sequences is considered in details. The Van der Corput transformation for sequence of natural numbers is analyzed, with its help the Halton points are constructed. It is shown that the cubature formula with one free node corresponding to the Halton sequence is exact for some class of step functions. The class is explicitly described. The obtained results allow us to use the mentioned sequences more effectively for integral calculation and extremum finding, and can also serve as a starting point for further theoretical studies in the field of the quasi-random methods. Refs 6. Fig. 1. Table 1.

Keywords: quasi-random sequences, Monte Carlo methods, Halton sequences, random quadrature formulas.

References

1. Kuipers L., Niederreiter H., Uniform Distribution of ¡Sequences (John Wiley & Sons, New York, London, Sydney, Toronto, 1974, 340 p.).

2. Ermakov S.M., Monte Carlo method and related matters (2th ed., Nauka Publ., Moscow, 1975, 472 p.) [in Russian].

3. Ermakov S. M., Monte Carlo method in computational mathematics. Introductory course (Nevskii dialekt Publ., St. Petersburg; Binom publ., Moscow, 2009, 192 p.) [in Russian].

4. L'Ecuyer P., Lecot C., Tuffin B., "A Randomized Quasi-Monte Carlo Simulation Method for Markov Chains", Operations Research 56(4), 958-975 (2008). https://doi.org/10.1287/opre.1080.0556

5. Chi H., ¡Scrambled Quasi-Random Sequences and their Application (A dissertation for the degree of Doctor of Philosophy. The Florida state university college of arts and science, 2004). Available at: http://diginole.lib.fsu.edu/islandora/object/fsu:182066/datastream/PDF/view (accessed July 5, 2017).

6. Niederreiter H., "Quasi-Monte Carlo methods and pseudo-random numbers", Bull. Amer. Math. Soc. 84(6), 957-1041 (1978). https://doi.org/10.1090/S0002-9904-1978-14532-7

Для цитирования: Ермаков С. М. О рандомизации квазислучайных последовательностей Холтона // Вестник СПбГУ. Математика. Механика. Астрономия. 2017. Т.4(62). Вып. 4. С. 570576. https://doi.org/10.21638/11701/spbu01.2017.404

For citation: Ermakov S. M. On the Halton quasi-random sequences randomization. Vestnik SPbSU. Mathematics. Mechanics. Astronomy, 2017, vol. 4 (62), issue 4, pp. 570-576. https://doi.org/10.21638/11701/spbu01.2017.404