Рекуррентные алгоритмы Монте-Карло

Содержание

• История вопроса

Развитие методов Монте-Карло связано с применением стохастических алгоритмов в области теории переноса излучения ([8], [13]). В частности, схема Неймана-Улама ([10]) возникла в связи с проблемой точности моделирования решения задачи переноса излучения. При этом схема Неймай^-УлаЖ' «состоит в моделировании цепей Маркова и, с другой стороны, тесно связана с итерационным процессом

Хп+1 = АХп +

Развитые в связи с этим алгоритмы расходовали небольшой объем машинной памяти, однако для их применения требовалось выполнение условия сходимости мажорантного процесса, что резко ограничивало область применимости схемы Неймана-Улама.

Это обстоятельство потребовало развития новых методов, применимых к более широкому классу уравнений. Существенный шаг в этом направлении был сделан в работах [5], [3]. Во второй из них исследование было проведено в применении к разностному аналогу волнового уравнения, сформулированы условия стохастической устойчивости, приведен пример неустойчивой схемы.

Понятие устойчивости играет большую роль в исследовании свойств стохастических алгоримов. Во многих случаях неустойчивость, проявляющаяся в экспоненциальном росте дисперсии оценок, делает применение схемы Неймана-Улама невозможным.

Альтернативой схеме Неймана-Улама могут служить, в частности, методы, использующие аппарат стохастических дифференциальных уравнений. Для решения эллиптических уравнений построены эффективные алгоритмы такого рода ([18]). Строго говоря, такие алгоритмы не укладываются в схему Неймана-Улама. Безусловно, представляет интерес разработка общего языка, позволяющего с единой точки зрения взглянуть на эти два типа стохастических алгоритмов. Определенный прогресс в установлении такой общности был достигнут в работе [4].

Однако, задачи этого типа достаточно сложны. Их анализ осуществляется сравнительно просто после введения дискретного времени. Заметим, что при решении практических задач дискретизация по времени является стандартным приемом и служит достаточно общим инструментом исследования. В особенности она важна, когда применяются метод расщепления или метод дробных шагов ([21]), позволяющие заменить исходную задачу последовательностью более простых задач.

Если исходное уравнение имеет вид

Ltu + v (t), (0.1) где ^ - линейный оператор, а и — искомая (вектор-) функция, то после дискретизации с шагом Д£ по времени получаем приближенные соотношения: ип+1 ип

Д£

Ьпип + Vй (явная схема) или (0−2) ип+1 ип ип+1 ип

--- = Ьп--1ип+1 + Vй (неявная схема). (0−3)

-А V

Если же представить в виде суммы = + ь[2 то можно построить явно-неявную схему:

0.4)

Полученное рекуррентное соотношение определяет рекуррентный характер связанных с ним стохастических алгоритмов.

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

Возникает также вопрос о предельном поведении серий последовательностей оценок рекуррентного алгоритма и связанных с ними марковских процессов при Д£ —> 0. В предлагаемой работе получены условия сходимости конечномерных распределений таких процессов к конечномерным распределениям некоторого диффузионного процесса. Полученные результаты позволяют проследить нюансы взаимосвязи предельных переходов по пространственным и временной переменным с основными классами оценок Монте-Карло.

В диссертационной работе получены некоторые обобщения результатов, опубликованных ранее в работах [3], [4]. Эти обобщения потребовали, в частности, введения некоторого нового класса цепей Маркова, с помощью которого строится и обосновывается рекуррентный алгоритм.

Структура диссертации

В главе 1 описан рекуррентный алгоритм построения несмещенных оценок? п, п > 0 величин ип, п > 0, удовлетворяющих рекуррентному соотношению ип+1 = 41)"п+1 + Ь&ип + vn+1 (0.5) с заданными значениями и°, уп, п > 0 и известными операторами т (1) г (2) ¦ип? -ип

Для этого сначала вводится специальный класс цепей Маркова, множество состояний которых — целочисленная решетка, ограниченная слева, справа и снизу. Траектории цепей устроены так: сначала движение по горизонтали, затем переход на верхний слой, после чего цепь обрывается. Движение вниз запрещено.

При этом параметры марковской цепи связаны с операторами г (1) г (2)

Ьп', Ьп условиями согласования.

По аналогии с классической схемой Неймана-Улама выделены четыре основные оценки: «сверху вниз», «снизу вверх», по поглощению и по столкновениям. Доказана их несмещенность, вычислены вторые моменты.

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

В главе 1 также подробно обсуждаются примеры применения рекуррентного алгоритма: для решения уравнений в частных производных, интегральных уравнений, а также систем алгебраических уравнений. В последнем примере применение рекуррентного алгоритма позволяет расширить по сравнению со схемой Неймана-Улама границы применимости методов Монте-Карло.

Глава 2. посвящена вопросам устойчивости рекуррентного алгоритма. При этом мы ограничиваемся исследованием спектральной устойчивости. Так, слабая стохастическая устойчивость определяется как свойство алгоритма, заключающееся в линейном росте нормы матрицы ковариации при росте п. Сильная стохастическая устойчивость определяется как ограниченность последовательности норм матриц ковариации оценок при всех п.

Таким образом, задача сводится к отысканию линейного рекуррентного соотношения, связывающего элементы матриц ковариа-ций оценок? п, п > 0: со^п+1 = 7гпсоу? п + рп. (0.6)

Для этого используется следующий прием: матрица ковариаций «растягивается» в вектор, составленный из ее элементов.

Устойчивость рекуррентного алгоритма, для оценок которого выполняется соотношение типа (0.6), определяется значением нормы операторов 71п и ограниченностью последовательности рп. Получены достаточные, а также в некоторых случаях необходимые и достаточные условия слабой и сильной стохастической устойчивости.

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

Операторы, входящие в соотношения (0.6), получены в явном виде для общего случая и для некоторых частных, более удобных в практических задачах.

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

В главе 3 рекуррентный алгоритм рассматривается в применении к системам линейных дифференциальых уравнений. После дискретизации по времени с шагом At они сводятся к (явно-неявной) рекуррентной процедуре.

Рассмотрим последовательность серий случайных величин, где каждая серия есть последовательность оценок £п величин, удовлетворяющих этой рекуррентной процедуре, полученная с помощью рекуррентного алгоритма при фиксированном Д£. Рассмотрим кусочно-линейный процесс с вершинами в точках п, пДг), п > 0.

В главе 3 изучаются условия сходимости конечномерных распределений процессов и распределений функционалов от £л к соответствующим распределениям марковского процесса являющегося решением некоторого стохастического дифференциального уравнения.

При некоторых условиях на параметры алгоритма получены достаточные условия сходимости и явный вид предельного стохастического дифференциального уравнения.

Полученные результаты служат теоретической основой для построения алгоритмов приближенного (в смысле сходимости конечномерных распределений) решения стохастических дифференциальных уравнений произвольной размерности. Пример такого типа завершает главу 3.

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

В этом смысле результаты главы 3 являются обобщением результатов, полученных в работах [11], [16], [12].

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

Показано, что оценка метода частиц является асимптотически несмещенной при N оо, где N — число частиц. Получена оценка порядка трудоемкости метода. Это позволяет сравнивать алгоритм метода частиц с другими стохастическими алгоритмами.

Об обозначениях

Перечислим обозначения, принятые в диссертационной работе.

1) Под точкой жЕ®-^ подразумеваем вектор-столбец. Для всякого вектора ж = (жх5., ж^) =

1. хт — транспонированный вектор (строка) —

2. |ж| - модуль вектора:

М = ух1 + --- + х1

3. Обозначим символом Кк к-й базисный вектор: < при % = 1,. с?, А- = 1,., с?. г [ 0 иначе

2) Для всякой матрицы, А порядка с1х. с1:

1. Ат — транспонированная матрица-

2. (Над, А — диагональ матрицы А, (Над, А Е

3. А (А) — максимальное по модулю собственное число матрицы А.

Единичную матрицу будем обозначать буквой I.

3) Обозначения, связанные с интегральными уравнениями.

1. Для всякой ограниченной функции /(ж) и заряда /х (с?ж) = ! /(х)р ((1х).

2. Дельта-меру, сосредоточенную в точке Х: будем обозначать

3. Символом [х будем обозначать полную вариацию заряда ц.

4. Для двух зарядов (мер) V символом /IX V будем обозначать их прямое произведение. Будем также использовать обозначение ц2 = /IX/I.

5. Символом X будем обозначать тождественный оператор.

6. Для всякого интегрального оператора, А символом Х (Л) будем обозначать наибольшее по модулю собственное число оператора Л.

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

Для любых двух векторов х, у G

1. [х, у]е — вектор, составленный из произведений соответствующих компонент векторов х и у: х, у]е = (xiyi,., xdyd)]

Для векторов, составленных из квадратов компонент вектора, будем применять обозначение

Х^ ^ = [ж, хе — (XI,. .. , X.

И вообще, для любой степени п будем обозначать х^- ^ = (Xi,., xd) .

2. | - вектор, составленный из частных соответствующих компонент векторов ж и у:

X /Хх хЛ у V2/i г/d/ '

3. а (х, у) — матрица dxd, составленная из попарных произведений компонент векторов х и у: а (х, у) =

Вместо а (х, х) будем для краткости писать а (х).

5) Введем обозначения для операций поэлементного умножения и деления матриц. Для любых dxd-матриц, А = = H&zjll

1. [А, В]е — матрица, составленная из произведений соответствующих элементов матриц, А ж В:

А, В]е = \а^Ь{ й

2. — матрица, составленная из частных соответствующих элементов матриц, А и В

А В агЗ

Ьц г, 3=1

Для матриц, составленных из квадратов элементов матрицы, будем применять обозначение

А&trade- = [А, А]е = а гэ 1 г, 3=1

И вообще, для любой степени п будем обозначать

АЫ = а^ к*

3. Для всяких матрицы, А и вектора х будем обозначать сI, А х агз х3

6) Пусть, А = Ца^-Ц — матрица с1хс1. Тогда будем обозначать 1. Вектор длины с12

А = (а>ц) Е

2. Матрица в (А) = \в (А)гЛз, Г1 1 Ог. чйт

При этом порядок, в котором элементы матрицы выписываются У в строку, не имеет значения, важно только, чтобы, А и в (А) были согласованы с точки зрения этого порядка. Для исследования устойчивости (глава 2) нам понадобятся специфические обозначения.

7) Пусть А, В- матрицы с1х<1, х €

1. Обозначим за 6(х) диагональную матрицу <1х (1, для которой сИад5(х) = х.

2. Обозначим за (3(х, А) матрицу с? хс?, для которой

3(х, А) ц — + ХуАц при всех г,] = 1,., с?.

8) Для матриц А, В порядка обозначим:

1. к (А) — матрица (12хс12 такая, что 4 |? ^ - 5 ^ | о иначе-

2. ш (А, В) — матрица с?2Хб?2 такая, что

АхдВуъ ~Ь AjsBij, 5 — cv{A, B) ij. s>t =

О иначе.

9) Введем еще два обозначения.

1. Для любой симметричной матрицы, А введем вектор длины 1)/2 string (А) =

А22,, Addi

-^-31?^42, , Ad, d-2,, Д*,] в строку последовательно выписаны элементы главной диагонали, умноженные на ½, затем нижней поддиагонали и т. д.) Вектор string (А) взаимно-однозначно задает симметричную МАТРИЦ, у А.

2. Для любого вектора х = (xi,., ж^) обозначим 1опд (х) = (а?!,., Xd, 0,., 0) — вектор длины d (d+l)

Т. е. вектор 1опд{х) получается из х дополнением нулей до d (d+1) ДЛИНЫ —

Последние два обозначения будут использоваться только графе 2.4.