Автоматика и телемеханика, Л- 5, 2007
РАС Б 02.60.СЬ, 02.70.-c. 02.70.Tt, 02.70.Uu
© 2007 г. „М.А. МАРЧЕНКО, канд. физ.-мат. наук, Г.А. МИХАИЛОВ, д-р физ.-мат. наук, чл.-корр. РАН (Институт вычислительной математики и математической геофизики СО РАН, Новосибирск)
РАСПРЕДЕЛЕННЫЕ ВЫЧИСЛЕНИЯ ПО МЕТОДУ МОНТЕ-КАРЛО1
Обсуждаются вопросы эффективной параллельной реализации некоторых алгоритмов метода Монте-Карло. Описывается параллельная модификация генератора базовых псевдослучайных чисел, равномерно распределенных в единичном интервале. Описывается технология распределенных вычислений в сети персональных компьютеров с использованием разработанного авторами комплекса программ МОКС.
При решении задач математической физики по методу Монте-Карло для искомого функционала у подбирается вероятностное представление в виде математического ожидания от некоторой случайной величины
где ш - траектория моделируемого случайного процесса. Здесь имеется в виду, что вероятностное представление может быть приближенным, т.е. соответствующая случайная оценка имеет ненулевую детерминированную погрешность. Величина Е £ приближенно оценивается выборочным средним для достаточного числа независимых реализаций Сп случайной величины
Однако случайные оценки часто имеют большую трудоемкость. Здесь под трудоемкостью С(£) понимается средний объем вычислительной работы, необходимый для достижения заданной вероятностной погрешности, который пропорционален величине
где t(() - среднее время ЭВМ, необходимое для моделирования одного выборочного значения, Б £ - дисперсия случайной оценки [1].
1 Работа поддержана грантами Российского фонда фундаментальных исследований 06-0100586 и 06-01-00046, грантом НШ-4774.'2006.1 программы "Ведущие научные школы", грантом Президента РФ МК-1777.'2005.1, грантом ШТЛ8 05-109-5267, грантом Лаврептьевского конкурса молодежных проектов СО РАН.
1. Введение
у « ЕС = ЕС(ш),
п=1
С (С )= КС) б с,
С целыо понижения трудоемкости моделирование независимых выборочных значений можно распределить по процессорам параллельной вычислительной системы. Очевидно, что главным критерием осуществимости такой параллельной реализации является возможность "поместить" вычислительную программу моделирования выборочного значения в память каждого процессора.
При распределении вычислений по процессорам следует допускать возможность реализации различных объемов выборки на различных процессорах с использованием статистически оптимального способа осреднения результатов по формулам вида
M м
С nmZm / Y^ nm■■
m=1 m=1
где M - число процессоров, nm - объем выборки для то-го процессора, Zm - соответствующее среднее значение. Таким образом, можно добиться обратно пропорциональной зависимости величины трудоемкости случайной оценки от числа процессоров (при условии, что используемые процессоры имеют одинаковую производительность) [2].
Описанную методику параллельного моделирования можно легко модифицировать с целью эффективной оценки функционалов, зависящих от параметра x G X:
v(x) к E Z(x; и).
При этом моделируемое распределение случайной величины и может зависеть, а
x
зовать для оценки функционалов, возникающих в результате осреднения базовых функционалов по случайному параметру а задачи, например по случайной плотности среды. Вероятностное представление при этом имеет вид двойного математического ожидания
V к EEZ(и,а).
Примеры подобных задач приведены в разделах 2 и 3.
Как правило, при параллельной реализации необходимый объем выборки базовых случайных чисел очень велик, поэтому целесообразно использовать "длиннопо-риодные" псевдослучайные последовательности. Удовлетворяющий этому требованию ''параллельный" генератор описан в разделе 4. Этот генератор был проверен с помощью содержательных статистических тестов.
Описанная методология параллельной реализации соответствует концепции использования инфраструктуры GRID, которая быстро развивается в последнее время [3]. Отметим, что существует достаточное число разработок GRID-систем, предназначенных для распределенных вычислений по методу Монте-Карло (см., например, [4 6]). В настоящей работе предлагается система MONC, обладающая большей простотой и удобством в использовании, чем имеющиеся аналоги (см. раздел 5).
2. Глобальная оценка решения
Для глобальной оценки функции
V(x) = J g(x,y)P (dy)
Y
в ограниченной области D можно оценить ее значения в узлах сетки с шагом h методом Монте-Карло и затем осуществить какое-либо восполнение, например линейное
непрерывное. Такую оценку обозначим через у(х). Вообще говоря, функция у(х) -это случайное поле, распределение которого связано с объемом выборки (т.е., с числом реализаций в методе Монте-Карло), с типом базовой оценки £ и с шагом Н. Задача сходимости оценки у к у в некоторой метрике может быть решена на
основе рассмотрения величины
у) = Е \\у(х) - у{х)\\ь{в),
где Ь(Б) - соответствующее Банахово пространство, так как Р(|£| > е) < Е \£\/е Уе. Эта задача связана с трудоемкостью алгоритма, т.е. со средним числом операций, которые обеспечивают выполнение неравенства В (у, у) < ¿.Оценить В (у, у) наиболее просто в пространстве Ь2(Б):
в2 (<,<)
1 f ^
E M [<(x) — <(x)]2dx
2
<
\ KD
2
^ j E[<(x) — <p(x)\l dx = j D <p(x)dx + j [<(x) — E <p(x)\l dx.
D D
у(х)
в Б. Тогда можно построить такую аппроксимацию сеточной функции, что выполнено соотношение
J[<(x) — E<(x)]2dx < C0h4
D
При этом
(1) в2 (<,<) < d/n + Coh4.
Цель настоящего раздела уменьшение трудоемкости функционального алгоритма метода Монте-Карло для заданного способа аппроксимации сеточной функции. Здесь используется лишь информация о порядке соответствующей аппроксимацион-ной погрешности в зависимости от шага сетки. Если этот порядок можно повысить, то следующие далее оценки изменяются очевидным образом.
В рассмотренных условиях задача приближенной минимизации трудоемкости может быть сформулирована в виде:
(2) S0 = nt0h-k ^ min, d/n + C0h4 = ö2,
n,h
где к - размерность пространства, n - обьем выборки, h - шаг сетки, смысл d виден из (1) и t0 = t mes (D), где t - затраты на одну реализацию mes (D) - мера Лебега множества D. Оптимальный порядок величин h, n и So таков:
h*0 х ö1/2, n*0 x ö-2, SO x ö-(2+k/2).
При этом B(<,<p) < ö и, следовательно, <р сходится к < в метрике Ь2.
C(D)
теоремы вложения пространства W2 (D) в пространство C(D) при 21 > к. Если используются производные первого порядка, то по необходимости к = 1, т.е. рассматривается решение одномерного уравнения (или многомерного уравнения на данной
линии). Здесь
\W - <K
j [p(x) — <p(x)]2dx + j [i'(x) — ip' (x)]2dx . D D
Поскольку дисперсия суммы независимых случайных величин равна сумме их дисперсий. то при независимой оценке значений решения в узлах решетки необходимо использовать неравенство
D ip'(x) <
d(x) nh2
Кроме того если ср' кусочно-постоянна, то J [сp'(x) - E ¿p'(x)]2dx < Ch.
Таким образом, асимптотически при h ^ 0
"с - ^ <% + Ch.
Эта оценка приводит к следующей задаче минимизации трудоемкости:
(3)
S*! = ntoh
-к
di
min, —2 + C1h2 = S2. n,h nh2
—>
Оптимальные порядки здесь таковы:
hl x S, n! x S-4, SI x S-5.
Итак, в одномерном случае трудоемкость глобальной оценки решения в С (В) квадратична сравнительно с оценкой в Ъ2(0). Трудоемкость оце нки в С (В) может быть существенно уменьшена путем зависимой оценки величины у, обеспечивающей соотношение
\i(x) — lp(x + h)\ <Ch, h ^ 0, с вероятностью 1. При этом вместо (3) возникает задача S2 = nt0h-1 ^ min, d2/n + C2h2 = S2,
n,h
решение которой дает
h*2 х ^ n2 х S-2, S2 х S-3.
Имеются различные способы коррелирования оценок по методу Монте-Карло [1]. Например, целесообразно использовать одни и те же псевдослучайные числа для различных точек, дополнительно распределяя эти числа специальным детерминированным способом между отдельными траекториями (см. раздел 4). При этом серии траекторий, исходящих из различных точек фазового пространства, можно моделировать на различных процессорах.
3. Метод "двойной рандомизации"
1. Далее рассматриваются рандомизированные оценки вероятностных моментов решений задач со случайными параметрами. Пусть уравнение Ьф = / решается методом Монте-Карло на основе моделирования некоторого случайного процесса с траекториями и. Это означает, что случайные величины £к(и) строятся так, что
Е 6 (и) = Лк, к = 1, 2,...,т,
где Лк — оцениваемые функционалы от ф.
Пусть оператор Ь и функция / зависят от случайного поля а (например, случайная среда в теории переноса, случайная сила в теории упругости и т.д.). Тогда
£к = £к(и,а), Лк = Лк(а)
Е[& (и,а)\а] = Лк (а),
где траектории и и а, вообще говоря, зависимы. Рассмотрим задачу оценки величин
Лк =ЕЛк (а), Якз =Е[Лк (а)Л^ (а)], к,] = 1,...,т,
Еа
Для решения этой задачи целесообразно использовать метод 'Двойной рандомизации", который строится на основе следующих представлений:
Е Лк(а) = ЕаЕш£к(и, а) = Е(ш,а)£к (и, а)
(4)
Е[ЛкХа)Лз (а)] = Е(Ш1,Ш2,а)[€к(и1,а)&(и2,а)}:
где и! и и2 условно-независимые траектории, моделируемые для одной реализа-а
ние. Ясно, что надо предположить существование полного математического ожидания, представленного в (4), т.е.
Е(Ш,<т)\£к (и,а)\ < Е(ш1,Ы2,а) [\£к (и1,а)^о (и2,а)\\ < +<»•
Соотношения (4) показывают, что для оценки Лк достаточно строить только одну траекторию для фиксированного а, а оценка величины Н^ требует построения двух условно-независимых траекторий. Для оптимизации такого алгоритма естественно использовать метод расщепления (см., например, [1]).
2. В качестве частного случая использования формулы (4), можно рассмотреть рандомизацию оценки по столкновениям (см. [1])
N
с = ^2 ЯпЦхп),
= 0
где
^ _/ (х0) ^ к(хи-1,хи)
Чо = —,—г, Чи = Чи-1—,-т.
п(хо) р(Хи-!,Хи)
6 Автоматика и телемеханика, о
161
Пусть ](хп— 1,хп), /о, Нп — независимые несмещенные оценки соответствующих величин к(хп—1,хп),/(х0),И,(хп) (например, случайные оценки интегралов, выражающих эти величины), <кп — несмещенные оценки весов Qn и К\ - интегральный оператор с ядром Е \](х' ,х)\. Если р(К1) < 1, Е \Н\ € Е \/\ € Ь1, то [1]
N
Е] = Е£, где ] = ^ ^
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.