Оптимизация функционального двустороннего геометрического метода Монте-Карло Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, Специальный выпуск, 2006

ОПТИМИЗАЦИЯ ФУНКЦИОНАЛЬНОГО ДВУСТОРОННЕГО ГЕОМЕТРИЧЕСКОГО МЕТОДА МОНТЕ-КАРЛО*

Т. Е. Булгакова Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия

e-mail: ceta@gorodok.net

A discrete-stochastic algorithm for approximation of an integral with laboriously calculated integrand which also depends on a parameter is investigated. To reduce the calculation cost it is proposed to use the double-sided geometrical method. Recommendations for the choice of conditionally optimal parameters of the algorithm are formulated. The results of test calculations are presented. These results confirm the correctness of the choice of conditionally optimal parameters.

1. Дискретно-стохастический алгоритм

В данной работе рассматривается задача приближения интеграла, зависящего от параметра

I(х) = I д(х, у) у, X е X е Ет. (1)

Изучается дискретно-стохастический алгоритм [1], в котором по параметру х вводится сетка {х15..., хм}; в каждом у зле х^ функция (1) вычисляется методом зависимых испытаний

/(х.) = [ /(У) ¿у « /П1(х0 = 1 £с}М), (2)

>У

j=l

где

Г(м) - аЫ- £) - 9

с а- т ,

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

м

I(х) И Ьм1(х) = ^2 1п 1 (х»)^г(х). (3)

* Работа выполнена при финансовой поддержке Президентской программы "Ведущие научные школы" (грант № НШ-4774.2006.1).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

По аналогии с [1] в качестве базисных функций {^¿(х)} возьмем конечные элементы Стренга — Фикса [2, 3] с кусочно-линейной образующей функцией (или "функцией-крышкой").

2. Двусторонний геометрический метод

Зафиксируем номер г и рассмотрим стандартный алгоритм (2) для случая, когда функция д(х^ у) как функция от /-мерной перемеиной у = (у1;..., у) является трудновычислимой. В работе [4] для этого случая в качестве альтернативы стандартному методу (2) представлен следующий алгоритм, названный .модифицированным, геометрическим методом Монте-Карло.

Принимаем, что 0 < д(х^ у) < ф(х^, у). Предлагается ввести дополнительную координату уг+1 и в (/ + 1)-мерном пространстве (у1,..., у1, у 1+1) = (у, у 1+1) рассмотреть область = У х (0,ф(хг, у)), а в ней — случайную точку (<^,...,6,6+0 = (£,6+0, распределенную согласно плотности

Муъ •••,У1, У1+1) = /г (у, У1+1) = ТГ^Т' у е¥-

ф,(хг, у)

Из этого соотношения следует, что вектор £ по-прежнему распределен согласно плотности f(у) в У, а случайная величина равномерно (условно) распределена па интервале (0,Ф(хг, £)).

Алгоритм 1. Реализуем п2 зпачепий /-мерного случайного вектора £ € У согласно плотности f (у).

2. Для ^'-й реализации £вектор а £ моделируем случайную в еличину £(+1, равномерно распределенную на отрезке [0,ф(х^,£где ] = 1,... ,и2.

3. Строим оценку интеграла:

1 П2

7(х0 « 4 (xi) = - £ (f'i], (f'i] = ФЫ, £w) ХУ(0 (£0), ей), (4)

n2 Л q 3 = 1

гДе Xv(i) (y,yi+i) — индикатор множества

Yq

Y(i) = {(y, У1+1); y G Y, yi+i G [0,q(Xi, y)]}

— подграфика функции q(xj, y).

Алгоритм 1 из [4] является модифицированным по сравнению с геометрическим методом из [5], который сформулирован для случая ^(xj, y) = const. Несложно показать, что

Eq(Xj, £) = Е^(хг, £) хусо (£,£i+i) = I(xj), т.е. оценки Z(1'j),Z(2'j) из (2) и (4) являются несмещенными и, кроме того,

Dq(Xj, £) < D^(Xj, £) Ху({) (£,£i+i), (5)

т. е. алгоритм 1 не дает уменьшения дисперсии по сравнению со стандартным методом Монте-Карло.

В работах [4, 5] указано, что применение этого алгоритма может оказаться полезным в следующих (достаточно редких) случаях: вычисление значений у) является трудоемким, а проверку неравенства

ей <д(х* Й (6)

(при выполнении которого Ху« Снг) = 1) удается упростить, заменив его на равносильное.

В работах [6, 7] представлена еще одна эффективная модификация геометрического метода — двусторонний геометрический .метод, который предполагает сокращение количества вычислений подынтегральной функции д(х^ у) при подсчете интеграла I(х^). Кроме мажоранты строим миноранту подынтегральной функции: ф(х^, у) < д(х^ у) < — (х^, у). Алгоритм двустороннего метода аналогичен приведенному выше алгоритму, только при проверке, попадает ли реализованная случайная величина (£(й), е+1) в подграфик Уд(г) подынтегральной функции д(х^ у), сначала проверяется, попадает ли эта точка в подгра-фик миноранты: < ф(£й)). Положительный результат избавляет нас от трудоемкой проверки неравенства (6). При этом в соотношении (4)

сГ = (ху«+(йей)), (7)

V ф /

гДе XV({) (у,У1+1) — индикатор множества уф

= {(у,У1+1); у е г, е [0, ф(х4,у)]}

— подграфика миноранты подынтегральной функции, а Хум (у, Уг+1) — индикатор множе-

Ф,ч

ства

У2 = {(у,У1+1); у е г, е [ф(хг, у),д(х*, у)]}

— области между графиком миноранты и графиком подынтегральной функции д(х^ у).

Если мажоранта — (х^, у) и миноранта ф(х^, у) достаточно близки, то значения подынтегральной функции д(х^ у) вычисляются редко. Следуя работе [7], будем предполагать, что мажоранта — (х^, у) и миноранта ф(х^ у) представляют собой кусочно-постоянные функции, определяемые значениями минимума ф™ и максимума -0™ функции д(х^ у) в ячейках сетки Вт, па которые разбита область иптегрировапия У.

3. Сравнение точности оценок и оптимизация по числу точек разбиения. Случай I = 2

Пусть для простоты У = [0,1] х [0,1] и / (у1, у2) = 1 при (у1;у2) е У (в этом случае д(хг,у1,у2) = д(х^ у2)). Строим равномерную сетку с шагом Н = 1/К Полагаем фт < д(х^ у2) < при (у1; у2) е Вт, где ячейка имеет вид

Вт = {(У1, уг) : (4™ - 1)Н < У1 < Н, (^2т) - 1)Н < У2 < Н}.

Здесь ¿1т), ^2т) — целые числа в пределах от 1 до К^, при этом т = 1,..., К2.

В этом случае оценка (4), (7) интеграла I(х^) выглядит следующим образом. Разыгрывается случайная точка (£ь£2,£3) согласно плотности /(у15 у2, у3), равной ^—г при

(у!,у2) ^ Вт и у3 € (0, ф^:*) и нулю иначе; соответствующие моделирующие формулы таковы: £1 = £2 = а2, и если (£ь£2) € Вт, то £3 = фт а3, где а — реализации стандартного случайного числа а (т. е. случайной величины, равномерно распределенной в интервале (0,1)).

Введем случайную величину £(2'г), равную ф^г при (£ь£2) € Вт и выполнении одного из условий: £3 < фт или фтт < £3 < д(хг, , £2); иначе С(2'г) = 0. Математическое ожидание случайной величины £ (2'г) равно

К2

ес (2,г) = Фтт х (р{£з < фй, £2) € вт >+р{фт} < £з <

к2

< д(Хг, £1, £2), (£1, £2) € Вт}) = фт) х Р{£з < д(Хг,£1,£2), (£1^2) € Вт}-

т=1

Вычисляя вероятность

Р{£3 < д(хг,£1,£2), (£1,£2) € Вт} =J J ДУъ^Ы ^3

вт о

д(ъ,У1,У2)(1у1(1у2,

фт ^

г В

получаем

к2

Б((2,г) = £ фЩ д(Хг, У1, У2) - 12(Хг) > Б((М), (8)

т=1 В

Вт

где Б£(1'г) = Бд(хг,£15£2) —дисперсия оценки для стандартного метода Монте-Карло (2):

к2

БС(1,г) = ^ I д2(Хг, У1, У2) ^2 - 12(Хг) = / д2(хг, ^1,^2) ^2 - 12(хг). (9)

т=1

Вт У

Таким образом, подтверждено соотношение (5) для рассматриваемого случая. Найдем условие на количество точек разбиения по каждой из координат области интегрирования так, чтобы при его выполнении оценка £(2'г) из (7) оказалась бы тем не менее предпочтительней оценки £(1'г) из (2). Будем предполагать, что подынтегральная функция д(хг, у2) удовлетворяет условию Липшица с константой Ь по координатам ^ и у2 (одной и той же для всех г)- Пусть а — время, которое затрачивается па сравнение двух чисел, Ь — время вычисления функции д(хг,у15у2) в некоторой точке (уьу2) (полагаем, что эта величина не зависит от г) и с — время получения стандартного случайного числа с помощью соответствующего датчика случайных чисел. Полагаем также, что а ^ Ь, с ^ Ь (т. е. функция д(хг, у2) является трудновычислимой — "равномерно" по г).

Затраты стандартного метода Монте-Карло (2) равны 51 = 2п1с + п1Ь. Здесь первое

£1 = а1 £2 = а2

затраты на вычисление значений д(хг,£1,£2).

Затраты двустороннего геометрического метода (4), (7) равны

4г) = Ь(К + 1)2 + 3аК2 + 3сп2 + п? (ар + (1 - р)(а + Ь)).

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

К2 К2 ,(»)

р = Р((£ь£2, £з) € у«) = Р{£з < Ф«}Р{(£ь £2) е вт} = Л2 ^.

т=1 ш=1 фт

Отсюда

"фгг

< + 1)2 + 3aK2 + Зеп? + n? ( a +

m=1 rm

2Lbh \

у'

где ф(г) = ш1п1<то<к2 -0ТО. Полагаем

+ I)2 + 3a.fi2 + Зсп? + п? (^а + < + щЬ, (10)

что соответствует случаю, когда затраты двустороннего геометрического метода (4), (7) меньше затрат стандартного метода Монте-Карло (2).

Согласно "правилу трех сигм", при приближении математического ожидания Е£ сред-

га

ним арифметическим (1/п) ^ £ для достаточно большого п с большой вероятностью

3 = 1

(около 0.997) выполнено соотношение 72 = Зл/0(/п, где 72 — опустимый уровень погрешности (индекс "2" нужен для дальнейших рассуждений — см. разд. 4). Используя последнее равенство и соотношения (8), (9), получаем

п1] = ~2 I I 92(*г, 2/1, у2) <1у1 ¿у2 - /2(хг) ] , \у

Л I g(xi,yi,y2)dy1dy2

Y2 * m=1 B

Отсюда

9 K2

. Y1 I ~ gУъ У2)) gdyi dy2'

Y2 m=lB

Bm

Используя условие Липшица для функции g(xj, yl5 y2), оцениваем разность ^г —g(xj, yl5 y2) < 2Lh и, следовательно,

(i) {г) ^ 18Lhljxj) 18L/(x,) щ — щ < -5-= Aih, Ai =-5-= const. (11)

Y2 Y2

Поэтому п^ < п1г) + А^/К^. Отсюда, предполагая, что К + 1 « из неравенства (10) получаем выражение

К?(3а + Ь) + К^(а + с - Ь)п}г) + --+ —^ < 0. (12)

ф (г) ф (г)

Для каждого конкретного случая соотношение (12) дает выражение для параметра Так как для метода (4), (7), как правило, выполнено соотношение К ^ п^, затраты на построение мажоранты и миноранты относительно невелики, и можно записать следующее асимптотическое выражение для параметра К

2 (а + Зс)Агф^ + 2ЬЬп(1) 2ЬЪАг Кг(а + с-Ъ)П1 +Кг-^-+ 7^Г<0-

Анализ корней левой части позволяет сделать вывод о том, что последнее неравенство имеет место при

Кг > ^У ('2ЪЬп?/ф« + (а + 3 с)А^2 + 8 ЫА^? (Ь -а- с)/ф®+

+ (2ЬЬп5г)/ф (г) + (а + 3с)Аг) ^ (2п1г)(& - а - с)) . (13)

Таким образом, в рассматриваемом случае при выполнении неравенства (13) двусторонний геометрический метод (4), (7) эффективнее стандартного метода (2).

4. Согласованный выбор параметров

В теории дискретно-стохастических методов важна проблема оптимального (согласованного) выбора соответствующих параметров [1, 8]. Для алгоритма (3) такими параметрами являются числа узлов М сетк и (х^} и числ о ^выборочных знач ений Ставится сле-

дующая задача условной оптимизации [1,8]: найти минимум трудоемкости ттМга1Б(М,п1) при Т(в) (М, п1) = 7, где 7 — фиксированное положительное число, а Т(в) — верхняя граница погрешности в норме функционального пространства В. Общая схема решения этой задачи такова: из соотношения для Т(в) (М, п1) один го параметров (например, п1)

М

для Б, при этом получается функция одного переменного (М), которая и исследуется на минимум.

Рассмотрим такую задачу для простейшего метода зависимых испытаний для С-подхода (т.е. для случая В = С(X)) [1, 8]. Полагаем, что трудоемкость имеет вид произведения Б(М, п1) = Н1Мп^, детерминированная компонента погрешности равна ) = Н2/М2/т, а стохастическая — 8^ = Имеем уравнение [1, 8]

^<«.«.> = 5^ + ^ = 7. <">

Соображения о выборе констант Н1, Н2 и Н3 сформулированы в [8]. Имеем

п1

Нз2

(7 - Н2/М2/т)2

и требуется найти минимум функции

§т = я,я|м

(Y - H2/M2/m)2'

Приравнивая нулю производную

dS (M) HiH| / Я2 /ш + 4

dM (7-Я2/М2/™)3 I 7 M2/'

m

для условно-оптимальных параметров получаем выражения

Mopt = V2, nopt = (15)

- _ ffiff2m/2ff2(m + 4)2+m/2 _2_m/2 16m™/2 7

Заметим, что если нас интересует только порядок по 7 оптимальных параметров Mopt и nopt, т.е. соотношения вида

Л/Г ^ ^ ^ —2 о ^ -2—m/2

Mopt ~ 7 , nopt ~ 7 , S opt ~ 7 ,

и трудоемкость S пропорциональна произведению Mn1; то достаточно приравнять детерминированную и стохастическую компоненты погрешности и получить требуемый порядок из соотношения (15).

Последнее соображение позволяет рекомендовать величину y2 = 7/2 (здесь 7 — по-прежнему допустимый уровень погрешности) при подсчете константы A из (11) для выбора параметров алгоритма типа (3):

м

I(x) И LmI(x) = 1«2 (xi)^i(x) (16)

i=1

с оценками (7) в узлах сетки.

В целом выбор параметров алгоритма (16) предлагается осуществлять следующим образом. Для заданного уровня погрешности 7 из равенства

#2 Нз 7

72

М2/™ ^/щ 2 выбираем

M =(2H2)m/2 7—m/2, ni = (2Нз)27—2 И n11} = ... = n1M) = щ. Согласно формуле (И) вычисляем

Аг = 18Ы^г\ 1ф(ъ)= i^y)dy.

Y

Последний интеграл несложно вычислить в случае, когда ф(х^, у) — кусочно-постоянная функция (при этом сетка по у берется произвольной и относительно "грубой"). Далее вычисляем п2 = п1 + Л^/К^, где

Щ = (2ЪЬт/ф® + (3 с + а)А^2 + 8 ЪЬА^Ъ -а- с)/ф®+

+ (2Ып!/ф® + (а + 3е)Л^ ^(2щ(Ь - а — с)) (17)

(см. формулу (13)); затраты Ь и константа Липшица Ь полагаются равными для всех узлов {хД. Выбираем

К = шах(К1,..., Км), п2 = тах(п21),..., п2,М^) и используем далее полученные значения в алгоритме (16).

5. Результаты численных экспериментов

Численные эксперименты по сравнению оценок интеграла, зависящего от параметра, в узлах сетки проводились на примере функции

д(х, у1, у2) = ехр(—х + у1 + у2), при этом I(х) = ехр(—х)(е — 1)2.

Эта функция не является трудновычислимой, поэтому в экспериментах при определении каждого значения функции использовалась задержка времени 7.81 ■ 10-5 с. Верхний уровень погрешности 7 задавался равным 0.03. При этом а = 2 ■ 10-8, Ь = 8.038 ■ 10-5, с = 4.01 ■ 10-8с. Использовались следующие оценки констант, участвующих в выражениях для оптимальных значений параметров М и щ". Н2 = 1.00321 Н3 = 1.7217. В таблице приводятся погрешности (^ и е2) и времена счета (^ и ¿2) для методов (3) и (16) соответственно при М = 8 п1 = 13174. Число КорЛ вычислялось по формуле (17). Исследовано изменение ситуации при отклонении от предложенного КорЛ.

Параметр ^opt = 65 К = 55 К = 75 К = 25 К = 95

«2 46737 44906 36444 82985 31545

ei 4.48 • 10"3 13.15 • 10"3 9.70 • 10"3 1.34 • 10"3 3.22 • 10"3

h 8.382 8.542 8.562 8.503 8.482

е2 2.48 • 10"3 4.69 • 10"3 12.03 • 10"3 7.62 • 10"3 3.43 • 10"3

¿2 4.356 3.896 4.937 5.825 7.050

Тестовые эксперименты подтвердили преимущество алгоритма (16) по сравнению с алгоритмом (3) и показали, что выбор параметров, описанный в разд. 4, дает близкую к оптимальной трудоемкость алгоритма (16).

Список литературы

[1] Михайлов Г.А., Войтишек A.B. Методы Монте-Карло. М.: Академия, 2006.

[2] Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.

[3] Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. М.: Наука, 1981.

[4] Войтишек A.B., Дятлова Е.Г., Мезенцева Т.Е. Геометрический метод Монте-Карло и его модификации // Матер. V Междунар. семинара-совещания "Кубатурные формулы и их приложения". Красноярск: Изд-во КГТУ, 2000. С. 46-54.

[5] Соболь U.M. Численные методы Монте-Карло. М.: Наука, 1973.

[6] Войтишек A.B., Ухинов С.А. Использование существенной выборки в методе Монте-Карло // Сиб. журн. вычисл. математики. 2001. Т. 4, № 2. С. 111-122.

[7] Каблукова Е.Г. Двусторонний геометрический метод Монте-Карло // Тр. конф. молодых ученых. Новосибирск, 2002. С. 76-81.

[8] Войтишек A.B. Дискретно-стохастические численные методы: Дис. ... д-ра физ.-мат. наук. Новосибирск, 2001. 264 с.

Поступила в редакцию 15 сентября 2006 г.