CC BY
69
2014 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 1. Том 1 (59). Вып. 1
МАТЕМАТИКА
УДК 519.644.2
ЭМПИРИЧЕСКАЯ ОЦЕНКА ПОГРЕШНОСТИ ИНТЕГРИРОВАНИЯ МЕТОДОМ КВАЗИ МОНТЕ-КАРЛО*
А. А. Антонов, С. М. Ермаков
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7/9
В работе рассматривается возможность применения теоретико-вероятностных методов к детерминированной процедуре оценки погрешности метода квази Монте-Карло вычисления многократных интегралов. Существующие методы оценки упомянутой погрешности являются неконструктивными. Известное неравенство Коксмы—Хлавки включает в качестве константы вариацию подынтегральной функции. Вычисление же этой вариации является более трудной задачей, чем исходная. Поскольку метод квази Монте-Карло использует в качестве оценки интеграла среднее арифметическое значение интегральной функции, то можно было бы ожидать, что распределение (в теоретико-числовом смысле) остатка приближенного интегрирования подчиняется нормальному закону. Дополнительная трудность состоит, однако, в том, что квазислучайные точки, рассматриваемые с теоретико-вероятностной точки зрения, являются зависимыми, и численная оценка 2-го момента остатка тем самым затруднена.
Авторы предлагают новый подход к проблеме оценки 2-го момента погрешности, основанный на полученных авторами результатах из теории случайных кубатурных формул. Приводимые в конце работы численные примеры свидетельствуют о перспективности разработанного метода оценки погрешности. Библиогр. 7 назв.
Ключевые слова: Метод Монте-Карло, квази Монте-Карло, доверительный интервал, случайные кубатурные формулы, функции Хаара, последовательности Соболя.
1. Введение. Известно, что важной особенностью метода Монте-Карло при вычислении интегралов является возможность эмпирической оценки погрешности (построение доверительного интервала). Метод квази Монте-Карло, строго говоря, такой особенностью не обладает. Известное неравенство Коксмы—Хлавки (см. [7]) требует оценки вариации интегрируемой функции и величины дискрепанса квазислучайной последовательности. Эти задачи существенно более трудные, чем исходная задача численного интегрирования. Последнее обстоятельство привело к использованию так называемых рандомизированных квазислучайных последовательностей (см. [5]). Их использование может, однако, существенно увеличить вычислительную работу.
Рассмотрим в качестве примера один из способов построения рандомизированных квазислучайных последовательностей (см. [6]). Пусть Х\, Х2,..., XN —ква-
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00769-а).
зислучайные ¿-мерные векторы, расположенные в ¿-мерном единичном кубе а Я.2,..., Дм —набор случайных векторов, равномерно распределенных в Если {У} обозначает вектор, компонентами которого являются дробные доли компонент вектора У, то сумма
1 м N
т= 1 п=1
служит оценкой интеграла / (X )йХ. Порядок погрешности будет равен 0(-^=Лщ), математическое ожидание величины £млг(/) есть искомый интеграл, и ее дисперсия может быть использована для построения доверительного интервала. Очевидно, что без использования рандомизации при том же объеме вычислительной работы мы имели бы погрешность порядка 0(), что в ра-3 лучше рандо-
мизированной оценки.
Возникает вопрос: нельзя ли обойтись без рандомизации при оценке погрешности? Поскольку квазислучайные последовательности имеют некоторые свойства, сходные со свойствами случайных чисел (аналог закона больших чисел, например, см. [4]), можно попробовать применить к ним формально статистические процедуры оценивания погрешности. Разумно предположить, что распределение в теоретико-числовом смысле средних арифметических /(-^п) близко к нормальному
со средним / (X )йХ. Что касается дисперсии, то ее оценка оказывается затруднительной в силу очевидной зависимости векторов {Хп}^^=1. Действительно, попытка формально оценить второй момент распределения средних приводит нас к необходимости оценивания ковариаций /(Хп) и /(Хп+&) для п = 1, 2,..., N и к = 1, 2,..., N — п. Очевидно, что корректная процедура оценивания требует привлечения более чем N значений подынтегральной функции, что опять же приводит к увеличению вычислительной сложности.
Заметим теперь, что сумма /(^0) может рассматриваться как реализа-
ция случайной кубатурной суммы + Квазислучайные узлы, од-
нако, обладают свойством расслоения. Если случайные точки попадают в гиперпараллелепипед объема ¿, содержащийся в единичном гиперкубе с вероятностью ¿, то из N квазислучайных точек такому параллелепипеду принадлежит [N/¿1 точек с вероятностью, близкой к единице. Для последовательностей Соболя, связанных с функциями Хаара (см. [1, 2]), указанное свойство выполняется при делении гиперкуба на параллелепипеды со сторонами, равными 2-т', I = 1, 2,..., в, при некоторых целых положительных ш8, и в этом случае зависимость между квазислучайными узлами интегрирования выражается в принадлежности с вероятностью 1 некоторых подмножеств узлов различным гиперпараллелепипедам. Указанные соображения делают оправданным привлечение к оценке погрешности квазислучайных методов интегрирования теории кубатурных формул со случайными узлами.
Далее известный результат (см. [3]) относительно дисперсии квадратурной суммы в случае точности соответствующей формулы для двух функций Хаара распространяется на общий случай и приводятся численные примеры, свидетельствующие о перспективности данного подхода.
2. Вывод формулы и анализ дисперсии. Пусть пространство X наделено мерой причем ^(Х) = 1. Разобьем X на N непересекающихся частей равной меры: X = Х\ и Х2 и ... и Хдг, Хг Г)Хз = 0Ч/г ф^ ц{Х\) = = ... = ц{Хм) = Пусть
функция / : X ^ К ^-интегрируема с квадратом (/ £ Ь2(^,)). Далее, пусть N = 2я. Построим N функций Хаара для Х1,Х2,..., Хм:
x G Xj ,j G {ist,..., ist + 2s-qi - 1}; x G Xfc, k G {ist + 2s-qi,..., ist + 2s-q*+1 - 1}; иначе;
где для i > 1 qi = [log2 i], ist = (i - 1) • 2s-qi+i - 2s + 1.
Утверждение. {Xi}^ — ортонормированная система в L2(d^,):
1) ¡ж х2 = 1 Vi;
2) fX XiXjd« = 0 Vi = j.
Используя аппарат случайных кубатурных формул, построим формулу, точную для N функций Хаара. В силу теоремы Ермакова—Грановского [3] она будет выглядеть, как SN[f] = в обозначениях Д(/, Q) = det Х2..., XN{xi)\f=i и A(Q) = A(xi,Q). Плотность совместного распределения узлов этой формулы будет равна ^(Д(д))2.
Определитель A(Q) равен сумме всевозможных произведений вида ±xi(xi1) X2(xi2)... xn(iw), где (ii, i2,..., iw) есть перестановка (1, 2,..., N), причем каждое такое произведение однозначно определяется своей перестановкой. Заметим, что если хотя бы две точки попадают в одно и тоже подмножество Xfc, то, поскольку для любого i значение Xi(x) постоянно Vx G Xfc, определитель A(Q) имеет две одинаковые строки и, следовательно, обращается в ноль. Это означает, что плотность совместного распределения узлов может быть ненулевой только тогда, когда все точки Xj попадают в разные подмножества Xfc. Дадим определение таким множествам.
Определение. Латинским множеством Lat(ii, i2,. .., iw), задаваемым перестановкой (ii, i2,. .., iw), мы будем называть множество, задаваемое следующим образом:
Xi = 1;
Xi = < -2—,
(х1, Х2,..., хм) £ £а1(«1, ¿2,..., ¿м) ^ V? £ (1,2,..., N) 3!к^ : ж^- £ Хд^
и при этом к^ = = ^2.
Определитель Д(/, ф) также является суммой произведений вида ±/(ж^) Х2(ж^2).. .хм(¿м). Если такое произведение не есть ноль, то оно равно ±/(ж^1) • См, где См > 0 — некая константа, не зависящая от перестановки (¿1, ¿2,..., ¿м). Если 1пу(«1, ¿2,..., ¿м) обозначает количество инверсий в перестановке (¿^¿2,...^м), то для точек (ж1, Ж2,..., жм) £ £^(¿1, ¿2,..., ¿м) будет выполнено
Д(/, Ф) =(-1)1п^•<2"-<*> • См • (/(Ж1) + /(Ж2) + ... + /(хм)) • С, Д(Ф) =(-1)1п^) • См • (1 + 1 + ... + 1) • С,
где константа С>0 отвечает за количество ненулевых произведений. Итак, кубатурная формула имеет вид
д(/,д) 1Л
а плотность совместного распределения узлов -
, \ I "Ж> если (жьж2,... ,ждг) <Е £^(гь«2, • • • ,«дг),
^>(х1 ,x2, • • • , ХЯ) Л
I 0, иначе.
Таким образом, нами доказана следующая
Теорема 1. Кубатурная формула Бя [/] с плотностью распределения узлов ^(хх, Х2, .. ., Хя) является случайной кубатурной формулой, точной для семейства функций Хаара Х1, Х2, .. ., Хя■
Хотя факт несмещенности оценки непосредственно следует из вышеприведенной теоремы, мы можем доказать его отдельно в качестве проверки:
Г
Е5дг[/] = / 5дг[/] • -д^" • • • • , хц)ц{<1х1 )ц{<1х2 ) ■ ■ .^((¿Ждг) =
Xя
1
N
X f (x1)^(dx1) + ... + J f (xN )^(dxN = J f (x)^(dx).
Ключевым вопросом здесь является величина дисперсии построенной квадратурной формулы. Заметим, что теорема Ермакова—Грановского дает точное значение для дисперсии в случае регулярных систем, но система Хаара таковой не является, что приводит нас к необходимости вычислять дисперсию непосредственно:
DSn[f] = J SN[f]dp - П Sn[f]dp I = A - B,
XN \XN J
где в силу несмещенности
B = I J f (x)y«(dx) X
Вычислим A: 1
A=j!p J (Яж1)+Яж2) + ■■■ +f(xN)fdcp =
J N • f 2(xi)l£ot(i1,i2,...,iN)(xi, x2, .. ., xn)^(dxi)^(dx2) . .. ^(dxN)+
XN
N N-2
N!
XN
+ 2 J f (xi)f (xj)1zat(ii,i2,...,iN)(xi, x2, .. ., xn)^(dxi)^(dx2) . .. ^(dxN) >
i<3 Xn '
~ I I * * ^xN)l^(dx1)fj/(dx2) .. ./л(dxN)+
N
N!
XN
+ N(N - 1) J f (xi)f (x2)l£0t(i1,i2,...,iN)(xi,x2,... ,xn)^(dxi)^(dx2). ..^(dxN)
XN
N N-2
m [N-C + N(N-1)-D],
Выражения С и Б посчитаем отдельно: С = J / 2(Ж1)1£04(41 )(Ж1, Ж2 , . . . , Жм )^(^Ж1 )^(^Ж2) . . . ^(¿Жм ) =
(ЛГ 1)1
N м -1
X
X"
_ (N - 2)!
Б = J / (Ж1)/(Ж2 )1£01(г1,г2,...,г" )(Ж1,Ж2, . . . , Жм )^(^Ж1 )^(^Ж2 ) . . . ^(¿Жм ) =
¡(х1)/(х2)ц((1х1)/л((1х2) =
Х1 =Х2
^ - 2)! I [ [ I
/(ж1)/(ж2)л*(сгж1)/х(сгж2) — / /(х1)/(х2)^((1х1)ц((1х2) > =
N м-2
X1 =Х2
2
^ - 2)! 1 ' ' 1 м
N м-2 _
/(ж)^(йж) | - J /(Ж1)^(ЙЖ1) J /(Ж2)м(ЙЖ2)
X /(ж)м(^ж)| - ^ /(ж)М(^ж)|
Таким образом, N м-2
N м
2
N!
= У (у f<yx)^l<ydx)^ >
что приводит к окончательному виду дисперсии:
2
1 м / \
е2(х)и(йх) — у | I т{х)и{ах)
Ю>ЗД] = А - В = 1 у /2(х)^х) " Е (7 ^ =
2 2 2 ^ е2(х)^х)--^ ( У /(х)/л((1х) \ +1 I у /(Х)[Л(0Х)1 -у ^ Т[х)щах)
= Юмс[Л + ^(а 1 + а2 + - . . + адг)2 - а2 - .. . - = ВМс[/] ~ У"^(а» - а3-)2,
гдеВмс[/] = 77 Jf2(x)jj1(dx) — | Jf{x)jj,{dx) J —дисперсия традиционного метода
X \x J
Монте-Карло, в то время как ai = jf (x)^(dx), i = 1, 2,..., N.
Xi
Таким образом, нами получен следующий результат.
Теорема 2. Кубатурная формула Sn [f] с плотностью распределения узлов y>(xi, x2, .. ., xn) имеет дисперсию, равную DSn [f ]. Эта дисперсия всегда не больше, чем дисперсия традиционного метода Монте-Карло.
Отметим также, что выполняется формальное требование DSn [f ] ^ 0, что немаловажно для описанного далее алгоритма.
3. Схема эксперимента и численные результаты. Полученная нами формула Sn [f ] имеет весьма простой вид. Сложность в ее практическом применении, однако, кроется в моделировании плотности y>(xi,x2,..., xn), что особенно актуально в случае высокой размерности. Алгоритм такого моделирования может иметь следующий вид:
1) компоненты вектора Y моделируются равномерно в каждом из множеств Xi,X2,..., Xn: y G U(Xi);
2) случайная перестановка (ii, i2,..., iN) выбирается равномерно из множества всевозможных перестановок чисел (1, 2,..., N);
3) финальный вектор X образуется под действием (ii ,i2,...,iN) на вектор Y.
Полученный вектор X будет распределен с искомой плотностью y>(xi, x2,..., xn).
Этому факту, безусловно, требуется некое формальное доказательство. Мы не будем его приводить здесь, поскольку не будем использовать этот алгоритм. Вместо этого мы остановимся на одном простом частном случае, который, тем не менее, исчерпывает собой подавляющее количество прикладных вычислений.
Пусть X есть d-мерный куб [0,1]d, а ^ — мера Лебега на нем. В качестве ХЬХ2, • • •, Хдг мы будем рассматривать параллелепипеды вида jj-) х [0,1] х ... х [0,1] (формально Хдг = 1] х [0,1] х ... х [0,1]). Это означает, что для N точек
с плотностью совместного распределения y>(xi, x2,..., xn) ровно одна должна иметь первую координату в отрезке [0, еще одна — в отрезке и так далее, на-
конец, ровно одна — в отрезке , 1]- Таким замечательным свойством обладают точки Соболя. Это означает, что мы можем использовать любые N подряд идущих точек из обобщенной d-мерной последовательности Соболя.
Зафиксируем общее количество вычислений функции, обозначив его за M. Для удобства примем M = k • N = k • 2s. Заметим, что для корректного сравнения результатов алгоритма с различными параметрами необходимо обеспечить постоянство M. Это означает, что при заданном M у нас есть два переменных параметра, k и s, но менять их мы можем только так, чтобы всегда выполнялось M = k • 2s. Далее, при выбранных k и s алгоритм оценивания будет выглядеть следующим образом:
1) возьмем M точек обобщенной последовательности Соболя xi с номерами i G
4,s = {1, 2,. .., k • 2s};
2) построим оценку искомого интеграла 5*м[/] = jjj f(xi)>
3) построим оценку дисперсии BS^^ [/] на основе величин &k,j = j f(xi)i где j = 1, 2,..., N, то есть мы оцениваем N коэффициентов a^j по k точкам каждый. На основе этой дисперсии строим доверительный интервал с уровнем доверия 95%.
Отметим двойственность выбора параметров к и е. С одной стороны, при в = 0 мы получаем не что иное, как традиционный метод Монте-Карло, а с ростом в мы ожидаем сокращения дисперсии. Следовательно, стоит выбирать в как можно большим. С другой стороны, при максимальном в имеем к = 1, то есть каждый из коэффициентов оценивается всего по одной точке, и встает вопрос о достоверности такой оценки. Компромиссным решением в этой ситуации, по-видимому, будет использование неких средних величин к и в.
Продемонстрируем рассмотренный алгоритм и его чувствительность к изменениям параметров к и в на численных примерах.
В качестве тестовой функции рассмотрим функцию вида /(хх, Х2,..., хо) = Пг=1 1+т • Константа выбрана таким образом, что значение интеграла от этой функции для любой размерности с! равно единице. Для наглядности мы будем рассматривать модуль истинного остатка в силу того, что оценки сверху симметричны относительно оси абсцисс. Мы будем сравнивать оценки, полученные при помощи формулы Бм [/] для различных к и в. Для всех последующих графиков приняты следующие обозначения:
• по оси ординат отложены значения модуля остатка интегрирования и полученных оценок;
• по оси абсцисс — количество вычислений подынтегральной функции;
• красным цветом выделен истинный остаток интегрирования;
• синим — эмпирические оценки, полученные предложенным методом с различными параметрами;
• зеленым — доверительный интервал Монте-Карло с уровнем доверия 95%.
20М) «ОЮ
Количество вычислений функции
Рис. 1. Эмпирические оценки метода квази Монте-Карло, одномерный случай.
2000 4000 6(300
Количество вычислений функции
Рис. 2. Эмпирические оценки метода квази Монте-Карло, пятимерный случай.
2000 4000 6000
Количество вычислений функции
Рис. 3. Эмпирические оценки метода квази Монте-Карло, десятимерный случай.
Как видно из иллюстрации 1, даже при s = 1 мы имеем значительно более точную оценку остатка по сравнению с Монте-Карло. Увеличение s ведет к дальнейшему уточнению оценки и лучше отражает поведение истинного остатка интегрирования.
С ростом размерности эффективность метода несколько падает, и выигрыш уже не выглядит настолько значительным. Так, на иллюстрациях 2 и 3 представлены случаи d = 5 и d =10. Однако становится заметным другое преимущество предложенного метода: оценка более устойчива при даже небольших s по сравнению с доверительным интервалом. Максимальное увеличение параметра s дает больший выигрыш, однако оценка приобретает ступенчатый вид, что не всегда желательно. Такой вид оценки может быть полезен, если нас устраивает оценивание всего в нескольких контрольных точках.
Литература
1. Соболь И. М. О распределении точек в кубе и приближенном вычислении интегралов // Журнал вычисл. матем. и матем. физ. 1967. 7(4). С. 784—802.
2. Соболь И.М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1969.
3. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975.
4. Halton J. H. Quasi-probability— why quasi-monte-carlo methods are statistically valid and how their errors can be estimated statistically // Monte Carlo Methods & Applications. 2005. 11. P. 203—350.
5. Kuipers L., Niederreiter H. Uniform Distribution Of Sequences. New York: Wiley, 1974.
6. L'Ecuyer P., Munger D., Tuffin B. On the distribution of integration error by randomly-shifted lattice rules // Electronic Journal of Statistics. 2010. 4. P. 950-993.
7. Niederreiter H. Random Number Generation and Quasi-Monte Carlo Methods. Dover: Dover Publications, 2006.
Статья поступила в редакцию 24 октября 2013 г.
Сведения об авторах
Антонов Антон Александрович — аспирант; anton.antonov@math.spbu.ru Ермаков Сергей Михайлович — доктор физико-математических наук, профессор; sergej.ermakov@gmail.com
EMPIRICAL CONVERGENCE BOUNDS FOR QUASI-MONTE CARLO INTEGRATION
Anton A. Antonov, Sergei M. Ermakov
St.Petersburg State University, Universitetskaya nab., 7/9, St.Petersburg, 199034, Russian Federation; anton.antonov@math.spbu.ru, sergej.ermakov@gmail.com
A possibility of a probabilistic approach to a deterministic error estimation procedure of quasi-Monte-Carlo method is presented. Existing estimation methods of the aforementioned estimation are non-constructive.
A well-known Koksma—Hlawka inequality includes a constant integrand variation, estimation of which is a more complex problem. Since quasi Monte-Carlo method estimates the integral with the mean integrand value, we can expect that the error distribution (in a probabilistic sense) converges to normal. An additional difficulty, however, is that quasi-random sequences, being treated as randomly generated, are statistically dependent, which complicates the numerical estimation of second moments.
Authors propose a new approach to the problem of estimating second order moments, based on new results, derived from a theory of random cubature formulas. Presented numerical examples indicate the superiority of the discussed error estimation method over classical ones. Refs 0. Figs 0. Tables 0.
Keywords: Monte Carlo and Quasi-Monte Carlo method, confidence interval, random cubature formulas, Haar functions, Sobol sequences.
