Диссертация (1150625), страница 16
Текст из файла (страница 16)
Спрактической точки зрения, параметр можно выбирать из диапазона = 2, . . . , 10.Обобщение алгоритма на другие численные типы с плавающей точкой (например, binary32)может быть получено тривиальным образом: для этого достаточно только скорректировать представление (3.11). Естественно, диапазон допустимых значений параметра изменится в зависимости от того, сколько бит отводится под мантиссу в рассматриваемом типе данных.8810.750.50.250110211111111012011101102111121021120111102101111210211102011021102111121111111101120111111211102101120111102111102111211101120202111011111110212201111111120101111211102111211111101021111111111110220202101112011111101021111121120111211112001011211201111121100.250.50.751Рисунок 3.7: Распределение точек Холтона3.4Гибридная битовая рандомизация для метода „блужданийпо сфере“В предшествующих разделах мы ограничивались рассмотрением исключительно задач численного интегрирования.
В этих условиях рандомизированный квази-Монте-Карло почти всегдаоказывается существенно (с точки зрения асимптотики сходимости) лучше, чем наивный илирасслоенный Монте-Карло, что, в общем, только подтверждает распространённую точку зренияи известные результаты (например, Шюрер Р. [23]). Однако в настоящий момент существуетбольшой пласт задач, в которых замена Монте-Карло на квази-Монте-Карло не только не эффективна, но и вообще неосуществима.
Это, в основном, задачи моделирования специальныхраспределений или случайных процессов. Это связано в первую очередь с тем, что последовательные точки квазислучайных последовательностей сильно коррелированы, что приводит кпадению скорости сходимости дисперсии или же, что более потенциально опасно, к сильному смещению оценок (см. Морокофф У. и Кафлиш Р. [70]). Универсальный способ адаптацииквазислучайных конструкций в этих случаях в настоящий момент неизвестен, поэтому в такихзадачах по-прежнему актуальны традиционные для Монте-Карло методы понижения диспер-8910.750.50.250211131120121031101010010211010200121011221112220040121000110001101011201120131021221012201130010010112201111102102212011102021010210001012021011121230211011021231113000001220020002102231111111012001102001012012200311210221031002021212101011202111200112310100.250.50.751Рисунок 3.8: Результат расслоения (точки Холтона) для = 3, мелкие ячейкисии, в частности расслоение.
Однако полученный нами результат об эффективности гибриднойбитовой рандомизации может служить ключом к успешному применению low-discrepancy последовательностей и в таких условиях. В частности, мы покажем, каким образом с его помощьюможно получить уменьшение дисперсии при использовании метода „блуждания по сферам“ длярешения задачи Дирихле.Рассмотрим постановку внутренней задачи Дирихле для оператора Лапласа. Пусть требуетсянайти решение уравнения∆ = 0(3.14)в ограниченной области ⊂ R , удовлетворяющее на границе Γ = условию⃒⃒⃒ = (3.15)Γдля наперёд заданной непрерывной функции (). Искомая функция должна быть дважды дифференцируема в области и непрерывна вплоть до границы: ∈ 2 () ∩ ( ∪ Γ).При наличии некоторых ограничений на область и границу Γ решение задачи всегда существует и единственно (Владимиров В.С., [71]).90Оценка дисперсии метода2-62-12МетодMCMC_stratHybrid_1Hybrid_2Hybrid_3Hybrid_4Hybrid_5Hybrid_10Hybrid_26RQMC2-182-242526272829210211212Количество вычислений подынтегральной функции213214215Рисунок 3.9: Гибридная битовая рандомизация, размерность = 5Одним из способов численного решения рассматриваемой задачи является так называемыйсферический процесс (процесс блуждания по сферам).
Введём рад обозначений:– Ω̄ – замыкание области Ω;– (, Γ) – расстояние от точки ∈ Ω̄ до границы Γ;– (, ) – расстояние между двумя точками и ;– Γ = { ∈ Ω̄ : (, Γ) < } – -окрестность Γ;– ( ) = { ∈ Ω̄ : (, ) = (, Γ)} – сфера максимального радиуса, целиком лежащаявнутри Ω̄.Сферический процесс определяется как однородная марковская цепь { }, выходящая иззаданной точки 0 . Переход между последовательными состояниями цепи задаётся следующимобразом: каждая следующая точка имеет равномерное распределение на сфере (−1 ). Цепьобрывается после шага , если ∈ Γ .Определённая таким образом цепь имеет несколько важных свойств (Мюллер М., [72], Ермаков С.М.
и Михайлов Г.А. [5], Ермаков, С.М., Некруткин, В.В., Сипин, А.С. [73]):1) траектория блуждания по сферам сходится к границе области с вероятностью 1;2) для каждого шага вероятность попадания в Γ ограничена снизу величиной2;42 (,Γ)3) среднее количество переходов не превосходит величины |ln()|, где – константа, независящая от .91Оценка дисперсии метода2-10МетодMCMC_stratHybrid_1Hybrid_2Hybrid_3Hybrid_4Hybrid_5Hybrid_10Hybrid_26RQMC2-152-202-252526272829210211212213Количество вычислений подынтегральной функции214215Рисунок 3.10: Гибридная битовая рандомизация, размерность = 20Пусть теперь необходимо получить оценку решения задачи во внутренней точке 0 . Это решение оценивается с помощью усреднения известных граничных значений в точках выхода понекоторому количеству смоделированных траекторий процесса.
Основанием для такой процедуры служит результат о несмещённости такой оценки, впервые полученный Мюллером М. в [72]на основе теоремы о среднем и формул Грина. Классический алгоритм Монте-Карло предполагает, однако, обрыв траекторий в окрестности Γ замену точки обрыва на ближайшую точку Γ,что вносит смещение в оценку в силу того, что > 0.Основная сложность для метода Монте-Карло по сравнению с численным интегрированиемзаключается в том, что в этом случае речь идёт о траекториях неограниченной длины, то есть,фактически, об интегралах, имеющих сколь угодно большую размерность. Это затруднение приводит к тому, что в настоящий момент не существует надёжной адаптации квази-Монте-Карло,которая не приводила бы к непредсказуемому росту смещения схемы.
Среди предпринимаемыхпопыток построения гибридных схем моделирования марковских цепей достаточно известен метод Array-RQMC (Лекуер П., Леко К. и Туффин Б., [74]), но этот метод предполагает наличиенесмещённой схемы и поэтому не подходит для рассматриваемой задачи.Ключевым шагом для генерации траектории сферического процесса является подзадача моделирования положения каждой следующей точки, равномерно распределённой на -мерной сфере.Для этого можно использовать следующий алгоритм, предложенный Марсальей Г. [75]. Возьмёмвектор = (1 , 2 , . . . , ), состоящий из нормально распределённых независимых компонент,√︀т.е.
∀ ∼ (0,1). Положим = 12 + 22 + . . . + 2 , тогда вектор равномерно распределён92на единичной -мерной сфере. Распределение на сфере произвольного радиуса моделируется израспределения на единичной сфере при помощи тривиального масштабирования.Описанный алгоритм будет служить основой для традиционного Монте-Карло подхода, скоторой мы будем сравнивать предлагаемую адаптацию.
Эта адаптация заключается в том, чторазные траектории будут моделироваться нами не независимо друг от друга, как это делается в классическом алгоритме, а исходя из идеи расслоения траекторий на нескольких первыхшагах, при этом остаток траектории моделируется обычным способом. Нами уже предложенспособ генерации такого рода расслоений на основе гибридной битовой рандомизации, поэтомуцелесообразно использовать его в данной задаче. В самом деле, пусть нам требуется смоделировать траекторий сферического процесса.
В случае Монте-Карло мы запустим независимыхпроцессов, где каждый переход происходит в соответствии с процедурой, описанной выше. Впредлагаемой гибридной вариации мы поступим иначе:1) Возьмём точек из последовательности Соболя размерности , где – размерность задачи, а – параметр алгоритма, отвечающий за число гибридных шагов:⎡⎢⎢⎣⎡⎢⎢⎣⎡⎢⎢⎣1,11,2...1,............,1,2...,+1,1+1,2...+1,............2,12,2...2,............(3.16)(−1)+1,1 (−1)+1,2 . . . (−1)+1,............,1,2...,2) Проведём гибридизацию массива по следующему правилу: первый блок из рядов получает случайных бит, во втором заменяется −1 бит и так далее, последний блок получает1 случайный бит.
Иными словами, последовательные блоки размером по рядов проходят93процедуру гибридной битовой рандомизации с параметрами = , = − 1, . . . , = 1:⎡⎢() ⎢⎣⎡⎢( − 1) ⎢⎣⎡⎢(1) ⎢⎣˜1,1˜1,2...˜1,............˜,1˜,2...˜,˜+1,1˜+1,2...˜+1,............˜2,1˜2,2...˜2,............(3.17)˜(−1)+1,1 ˜(−1)+1,2 . . . ˜(−1)+1,............˜,1˜,2...˜,3) Каждый из столбцов высоты преобразуем при помощи функции квантилей (обратнойфункции распределения) стандартного нормального распределения.
Отмасштабируем этиточки на сферу желаемого радиуса. Получим точки , ( – номер шага, – номер траектории), которые задают переходы для первых шагов по всем траекториям: 1,1 1,2 . . . 1,.........(3.18)... ,1 ,2 . . . ,4) Для тех траекторий, которые не оборваны после шага , последующие переходы моделируются методом Монте-Карло.Перечислим параметры алгоритма Монте-Карло (MC) и предложенного гибридного квазиМонте-Карло (HQMC). В рамках последующего численного эксперимента мы зафиксируем всеимеющиеся параметры, кроме размерности, по которой будем проводить перебор значений. Значения параметров алгоритма указаны в таблице 3.1.Таблица 3.1: Параметры алгоритмов моделирования сферического процесса и диапазонырассматриваемых значенийПараметрРазмерность задачи Точность выхода на границу Количество траекторий Количество внешних повторов Количество шагов гибридизации MC++++−HQMC+++++Перебираемые значения10, 200.0001100012010Как и в экспериментах в задаче численного интегрирования, внешние повторы необходимыдля того, чтобы нивелировать эффект случайности для одного запуска алгоритма.
Плюс илиминус в ячейках таблицы означает, присутствует ли данный параметр в этом алгоритме или нет.94Для демонстрации работы двух рассматриваемых алгоритмов мы рассмотрим конкретныйслучай задачи Дирихле для оператора Лапласа, подразумевающий специфическую область Γ испециального вида функцию (), гармоническую в этой области. А именно, положим() = ‖ − ‖2−(3.19)для некоторой фиксированной точки ∈ и > 3. Убедимся, что функция является гармонической всюду, кроме окрестности сингулярной точки . Первая частная производная имеетвид2−() =2(︃∑︁( − )2)︃ 2−−122( − ) ==1(︃= (2 − )( − )∑︁)︃− 2( − )2.=1Аналогично, вторая частная производная равна(︃ )︃− 2)︃− +22∑︁∑︁ 22() = (2 − )( − )− (2 − )( − )(− )( − )2( − ) =22=1=1(︃ )︃− 2 (︃)︃(︁ ∑︁)︁−1∑︁222= (2 − )( − )1 − ( − )( − ),(︃2=1=1и, суммируя по ,⎛∑︁ 2(︃∑︁() = (2 − )( − )22=1=1)︃− 2∑︀⎞2( − ) ⎟⎜=1⎜⎟⎜ − ∑︀⎟ = 0.⎠⎝( − )2=1Следовательно, достаточно вынести точку сингулярности за границу области, чтобы задачабыла поставлена корректно.
Поэтому определимΩ = ∖ (, ),где – единичный гиперкуб, (, ) – -мерный шар с центром в точке и радиусом . Вдальнейших экспериментах во всех рассматриваемых размерностях зафиксируем точку сингулярности = (0.7, 0.7, . . . , 0.7), а радиус выберем таким, чтобы объём этого шара был равен18(он будет различным в разных размерностях).