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

2004_ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА_Сер. 1_Вып. 2

МАТЕМАТИКА

УДК 519.676

А. В. Адамов, С. М. Ермаков

О СТОХАСТИЧЕСКОЙ УСТОЙЧИВОСТИ МЕТОДА МОНТЕ-КАРЛО

(случай операторов)

Введение

Данная статья является продолжением работ [1]—[4]. В этих работах рассматривались, главным образом, вопросы построения и стохастической устойчивости рекуррентных алгоритмов для получения несмещенных оценок векторов Уп, удовлетворяющих рекуррентной процедуре

Уп = АоУп + АУп-1 + /п, п > 1 (1)

с заданными векторами У0, /п € К и йх ¿-матрицами Ао,Ах.

Одной из важных характеристик рекуррентного стохастического алгоритма является поведение еоу(£п) с ростом п. Если существует постоянная С такая, что выполняется условие

||еоу(£п)|| < С для всех п > 1, то алгоритм называется стохастически устойчивым.

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

В данной статье некоторые из результатов работы [1] обобщены на рекуррентную процедуру

= Ао^п + А^п-1 + Iп, п > 1, где Ао, А1 —линейные интегральные операторы, — функции, заданные на мно-

жестве О С КЛ Предложен рекуррентный алгоритм для получения последовательностей оценок фп функций (рп, использующий для аппроксимации функций на каждом шаге процедуру непараметрического оценивания. Получены достаточные условия стохастической устойчивости предложенного алгоритма. Заметим, что поскольку оценки фп являются смещенными, то для устойчивости алгоритма, кроме условия

||соу(<уЗп)|| < С для всех п > 1

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

© А. В. Адамов, С. М. Ермаков, 2004

Схема Неймана—Улама

Классическая схема Неймана—Улама позволяет получить оценки линейного функционала (у, к) от решения уравнения

у = ^ + I, (2)

где А — линейный интегральный оператор, функции у, I являются элементами некоторого банахова пространства Б, к € Б *. Мы будем рассматривать случай

Ау(х) = J а(х,у)у(у)л(<у), х,у € Б С Кл,

где л — а-конечная мера, Б С Ь'2(<л), (•, •) —скалярное произведение. Как известно, схема Неймана—Улама применима в случае

II А || < 1, где Ау(х) = У \а(х,у)\у(у)л(<у). (3)

Если условие (3) выполнено, мы можем по траектории некоторой марковской цепи построить оценки £ = £(у, к) такие, что

Е £(у,к) = (у, к) и Е £2(у,к) <

Процедура аппроксимации функции у в Б, которую мы будем рассматривать в данной работе, состоит в следующем. Реализуются алгоритмы метода Монте-Карло в схеме Неймана—Улама для оценки значений величин {(у, к) г = 1,..., М}, а затем по полученным приближенным значениям {£(у, Ы)} с помощью некоторого преобразования О происходит восполнение функции у. В итоге оценка функции у выглядит следующим образом:

у(х)= О(£(у,кг),...,£(у,кК ),х), х € Б. (4)

Замечание 1. Частным случаем оценки (4) является оценка ядерного типа

1 И

Ф(х) = — ^2К(х,Хг)у(Хг)/д(Хг), (5)

г=1

где К(х, г) — функция, заданная на множестве Б х Б и близкая к ¿-функции 6(х — г); х\,...,х^ —независимые реализации случайного вектора х в Б, имеющего плотность распределения д(х), у(х^) —оценка по поглощению значения функции у в точке х^.

Поскольку у является решением уравнения (2), случайный оператор (обозначим его [(I — Ао)-1]), сопоставляющий функции I оценку у, можно рассматривать в качестве оценки оператора (I — Ао)-1.

Условие (3) ограничивает во многих случаях применение схемы Неймана—Улама. Кроме того, существует важный класс задач, где схема Неймана—Улама также не может быть применена на практике, несмотря на то, что все условия ее применения формально выполнены. Например, рассмотрим итерационный процесс

уп = Аоуп + А1уп-1 + Iп, п > 1, (6)

где Ао, А1 —линейные интегральные операторы.

Введем матрицу-оператор (размерности п х п)

Л(п)

и вектор-функции (размерности п)

/ /п \

/ п-1

р(п)

Ао А1 0. .. 0 0

0 Ао А1 . .. 0 0

0 0 Ао . .. 0 0

0 0 0 . .. Ао А1

0 0 0 . .. 0 Ао

/ 2

\Аофо + / V

ф(п)

( фп \

фп-1

ф

Н(п)

0 0

V0/

\ ф1 /

Тогда итерационный процесс (6) записывается в виде

ф(п) = л(п)ф(п) + р(п), п > 1. (7)

Применяя схему Неймана—Улама к уравнению (7), можно получить оценки

£ = £(ф(п), Н(п)) такие, что Е £ = фп, к) и Е £2 < +ж,

и затем построить оценку фп функции фп.

Заметим, что условие применимости схемы Неймана—Улама (3) для уравнения (7) сводится к условию У Ао|| < 1, и мы формально можем применять схему Неймана— Улама, даже если || А1Ц > 1. Однако, в последнем случае возможна ситуация, когда ковариация оценки фп будучи конечной при фиксированном п с ростом п будет расти экспоненциально. Поэтому описанный выше метод становится неприемлимым при достаточно больших п.

Следовательно, если мы хотим строить оценки фп для больших п, мы должны интересоваться поведением еоу(фп) при п ^ ж, при этом важно, чтобы метод был стохастически устойчивым, то есть чтобы выполнялось условие

||еоу(фп)|| < +ж при п ^ж.

Рекуррентный алгоритм

Рассмотрим рекуррентную процедуру

Фп = Аофп + А1фп-1 + /п, п > 1, (8)

где ф0,/п € Б — заданные функции.

Будем строить рекуррентную последовательность оценок фп для функций фп, удовлетворяющих соотношению (8). Положим

ф0 (х) = ф0(х).

Пусть построена оценка фп-1 функции фп-1. На п-м шаге строим операторные оценки [(I — Ао)-1А1]п и [(I — Ао)-1]п с помощью процедуры, описанной в пункте 1, и полагаем

фп(х) = (I — Ао)-1А1

ф

п1

(х) + (I — Ао)-1 /п(х).

(9) 5

п

= ЕО(((1 - Ао)-1А1уп-1 + (1 - Ао)-1/п, Ы), 1 < г < М,х) +

Теорема 1. Пусть выполняется условие

\\(1 - Ао)-1А1|| < 1. (10)

Тогда существует число N0 такое, что при любом N > N0 рекуррентная процедура (9) стохастически устойчива.

Доказательство. Введем обозначение для ошибки

£п(х) = фп(х) - рп(х)

и исследуем поведение первого и второго момента ошибки еп с ростом п. Имеем

Ееп(х) = Е фп(х) - рп(х) = 5п(х), (11)

где §П — систематическая погрешность на п-м шаге.

Обозначим Тп — а-алгебра, порожденная фп, Тъ — а-алгебра, порожденная Н^, 1 < г < N, Еп и Еъ — условные математические ожидания относительно Тп и Тъ соответственно, ¿1 —оператор ошибки восполнения

ёрх) = О((р, Н1),..., (р, Нм), х) - р(х).

Используя тот факт, что О является линейным преобразованием, имеем

Е фп(х) = ЕО(£(р, Н) 1 < г < М,х ) =

= ЕО(ЕкЕп-1^(р, К), 1 < г < N х ) =

= ЕО(Ек((1 - Ао)-1А1 фп-1 + (1 - Ао)-1/п, Ы), 1 < г < М,х ) =

\о)-1А1 рп-1 + (1 - Ао) + ЕО(((1 - Ао)-1А1^п-1, Ы), 1 < г < N х ) = рп(х) + ё1 рп(х) + + (1 - Ао)-1А1ёп-1(х) + ё1(1 - Ао)-1А1ёп-1(х). Отсюда и из равенства (11) получаем рекуррентное соотношение для погрешностей

ёп(х) = ((1 - Ао)-1А1 + ё1(1 - Ао)-1А1)ёп-1(х)+ё1 рп(х). (12)

Замечание 2. Для оценки ядерного типа (5) оператор ¿1 имеет вид

1 1

5цр(х) = К(х,хн)(р(хн)/д(хн) - у(х) =

г=1

= Jр(г)К(х,г)3,г - р(х). (13)

Из (13) следует, что если функция близка к ё-функции ё(х - г) и функция р удовлетворяет условию Липшица на множестве Б, то величина ё1р(х) равномерно близка к 0 при всех х € Б.

Относительно восполнения будем предполагать, что 11ё 1 \ стремится к 0 при N ^ ж. Тогда учитывая условие (10) можно утверждать, что найдется такое число N1, что выполняется неравенство

(1 + \\ё1\\) \\(/ - Ао)-1А1\\ < 1. (14)

Из соотношения (12) имеем

п 1 - (1+ \ё1\) ||(1 - А0) 1А1|| п

Таким образом, малость нормы ошибки восполнения обеспечивает малость нормы систематической ошибки.

Считая ошибку восполнения малой, будем анализировать поведение вторых моментов случайной ошибки. Обозначим

Еп(х,у) = соу(еп(х),еп(у)).

Покажем, что величина Нп подчиняется некоторому рекуррентному соотношению

Яп(х,у)= и Еп-1(х,у)+рп(х,у), (15)

причем для оператора и = иN выполняется условие ЦиЦ < 1 при N > Ыо. Очевидно, что

Кп(х,у) = соу^^Ы),

где

^(х) = фп (х) — Е фп(х).

Для случайной ошибки еп выполняется

еп(х)= (I — Ао)-1А1 еп-1(х) + (I — Ао )-1А1 Е фп-1(х) +

+

(I — Ао)-1 /п(х) — Е фп(х) = ДпЕ фп-1(х) + ДШп(х) +

+ (I + 5^ — Ао) А^п (х) + Дп^п-1(х),

где случайные операторы Дп и Дп имеют вид

Дп

Дп Д2

(I — Ао)-1А1

Е

п

1

(I — Ао)-1А1

п

(I — Ао)-1 — Еп-1 (I — Ао)-1

Учитывая, что Еп-1Дп = 0, г = 1, 2, и Ееп(х) = 0, из (16) получаем

Кп(х, у) = со^(I + 61)Ц — Ао)-1А1 еп-1) (х, у) + + соу(Дпеп-1)(х,у) + рп(х,у),

где

рп(х, у) = соу(ДпЕ фп-1 + ДЩп)(х, у)

(16)

(17)

(18)

Далее для интегрального оператора А определим оператор М(А)

М(А)г(х,у) = J! а(х,и)г(и,ю)а(у,ю)л(йи)л(йю), где функция г € ^(¿л®л), тогда

соу( (I + 5^ — Ао)-1А1еп) (х,у) =

= М ((I + 5l)(I — Ао)-1А^Кп-1(х, у).

Замечание 3. Используя неравенство Гельдера, нетрудно показать, что имеет место равенство ||М(А)|| = ||А||2.

Далее мы будем рассматривать такие восполнения, для которых выполняются следующие условия: существует линейный оператор С = См такой, что для любой функции д € Б имеет место равенство

Еп-1Дп1д(х)Дп1д(у) = Сд(х)д(у), (19)

П

71

п

п

-1 71

и

\\L\\^ 0 при N ^ж. (20)

Замечание 4. Можно показать, что оценка ядерного типа (5) обладает свойствами (19)-(20).

Из (17)-(19) получаем

Rn(x, y) = (m((! + W - Ao)-1Ai) + с) Rn-1(x, y) + pn(x, y).

Таким образом, оператор U из рекуррентного соотношения (15) имеет вид

U = M ((I + Si)(I - Ao)-1Ai) + с.

Принимая во внимание условия (14) и (20), а также замечание 3, имеем

\\U\\ ^ \\M((I + Si)(I - Ao)-1Ai)\\ = \\(I + Si)(I - Ao)-1Ai\\2 < 1. □

Summary

A. V. Adamov, S. M. Ermakov. On the stochastic stability of Monte Carlo method (operator case).

Problems of construction and stochastical stability of recurrent Monte Carlo algorithms for approximation of the functions satisfying a recurrent procedure

n л n I л n — 1 I rn ^ -I

ip = Aoip + Anp + f , n > 1,

where Ao,A:l are linear integral operators are considered.

A recurrent algorithm is introduced, where a nonparametric estimation procedure for approximation of functions on each iteration levels is used. The sufficient conditions of stochastic stability in the sense of boundedness of a covariation function of the estimators are obtained.

Литература

1. Ermakov S. M, Wagner W. Monte Carlo difference schemes for the wave equations // Monte Carlo Methods and Appl. 2002. Vol. 8. N 1. p. 1-29.

2. Гладкова Л. А., Ермаков С. М. Рекуррентные алгоритмы Монте-Карло для решения кинетических уравнений // Статистические модели с приложениями в эконометрике и смежных областях. СПб.: изд-во НИИХ СПбГУ, 1999. С. 50-75.

3. Ермаков С. М., Гладкова Л. А. О стохастической устойчивости рекуррентных алгоритмов методов Монте-Карло // Математические модели. Теория и приложения. СПб.: изд-во НИИХ СПбГУ, 2002. Вып. 2. С. 132-149.

4. Вагнер В., Ермаков С. М. Стохастическая устойчивость и параллелизм метода Монте-Карло // Докл. Акад. наук. 2001. 179 №4. С. 439-441.

5. Ermakov S. M. Stochastical stability, Neumann-Ulam scheme and particle methods // IVth IMACS Seminar on Monte Carlo Methods: Abstracts (September 15-19, 2003, Berlin). Berlin: WIAS, 2003. p. 55.

Статья поступила в редакцию 14 октября 2003 г.