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

УДК 519.245

Вестник СПбГУ. Сер. 1. Т. 1 (59). 2014. Вып. 4

МЕТОД МОНТЕ-КАРЛО И АСИНХРОННЫЕ ИТЕРАЦИИ*

А. В. Дмитриев, С. М. Ермаков

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9

Асинхронные методы, так же как и методы Монте-Карло, играют важную роль при многопроцессорных вычислениях, и исследование связей между ними, а в особенности расширение области их применимости, являются важными задачами. Рассматривается система вида х = Ах + /, где х — неизвестный вектор-столбец длины п, А — матрица п X п, / — вектор правых частей длины п. Случай, когда первое собственное число матрицы |А| по модулю меньше единицы, хорошо изучен. В этом случае установлена тесная связь между асинхронными итерациями и методом Монте-Карло. Если же первое собственное число матрицы |А| по модулю больше единицы, в то время как первое собственное число матрицы А по модулю меньше единицы, то схожие методы оказываются неприменимыми. Предложены модификации методов Монте-Карло и асинхронных итераций с частичной синхронизацией для решения задачи в указанном случае. Получены оценки погрешностей предложенных алгоритмов и условия сходимости предложенных методов. Проведены численные эксперименты. Библиогр. 6 назв. Ил. 4.

Ключевые слова: Метод Монте-Карло, асинхронные методы, асинхронные итерации, параллельные алгоритмы, статистическое моделирование.

Введение. Необходимость синхронизации операций, как известно, приводит к простою отдельных элементов оборудования многопроцессорных вычислительных комплексов. Это обстоятельство привлекло внимание исследователей к алгоритмам, не требующим синхронизации, так называемым асинхронным алгоритмам. В частности, были предложены асинхронные варианты методов итераций для решения систем уравнений (см. [1, 2]). Естественными свойствами асинхронности обладают также многие разновидности метода Монте-Карло для решения систем уравнений, и представляет интерес проследить связи и общность свойств этих разновидностей и детерминированных асинхронных методов, упомянутых выше.

Заметим также, что требование асинхронности налагает определенные ограничения на класс решаемых задач. Мы покажем ниже, как эти ограничения можно ослабить за счёт частичной синхронизации. Оказывается, что во всех случаях методика частичной синхронизации является общей, хотя метод Монте-Карло позволяет, вообще говоря, использовать больший арсенал приемов конструирования асинхронных алгоритмов.

Ограничимся рассмотрением методов итераций для решения систем линейных алгебраических уравнений вида

х = Ах + /, (1)

где х = (хх, Х2,..., хп)т — вектор-столбец неизвестных, / = (/1, /2,..., /п)Т — вектор правых частей и А = УП"^=1 —матрица системы.

Введем следующие обозначения: компоненты вектора х из Кп будем обозначать х^, г = 1,...,п, а последовательность векторов из Кп будем обозначать х(3), 3 =0, 1,. ...

Будем использовать следующие нормы:

* Работа выполнена при финансовой поддержке РФФИ в рамках научного проекта (грант № 14-01-00271-а).

• для вектора х = (х1,..., хп)

11х11 = тах ^

1<г<п

• для матрицы А = Ца^- 11^=1

п

(3)

1. Асинхронные итерации. В [2] Боде ввел определение асинхронных итераций для решения систем (в том числе и нелинейных) уравнений. Перепишем это определение для линейного случая.

Определение 1. Асинхронными итерациями для 'решения системы уравнений (1) с начальным вектором х(0) называется последовательность х(7), 7 = 0, 1, ..., векторов из Мп, определяемых соотношением

где 3 = {Jj | 7 = 1, 2,... } — последовательность непустых подмножеств {1, 2,... п} и Б = {(в 1 (), «2(7),..., «п(?))Ь' = 1, 2,... } — последовательность из Мп. Причем 3 и Б удовлетворяют следующим условиям:

• «¿(7) < 7 - 1 для 7 = 1, 2,...,

• «¿(7) является функцией от 7 и стремится к бесконечности, когда 7 стремится к бесконечности,

• г появляется бесконечное число раз в в множествах Jj ,7 = 1, 2,...

Асинхронные итерации для задачи (1), начинающиеся с х(0), для определенных 3, Б обозначим пятеркой (А, /, х(0), 3, Б).

Определению асинхронных итераций можно дать следующую интерпретацию. Пусть вычисления осуществляются на многопроцессорной системе и ¿1 < ¿2 < ¿з < • • • <tj < ... — моменты времени, в которые система производит вычисление очередной итерации. В момент времени tj один из процессоров вычисляет значение итерации х(7), которая отличается от х(7 — 1) набором компонент ^¿|г € Jj}. Эти компоненты вычисляются на основе компонент, полученных на предыдущих итерациях. Последовательность Б определяет номера предыдущих итераций и компоненты, используемые для вычисления х(7). Выбор компонент может быть обусловлен любым критерием, в том числе самым естественным — всегда использовать значения недавно вычисленных компонент.

Такая схема подразумевает отсутствие синхронизации между процессорами, однако для сходимости асинхронных итераций к решению задачи (1) при любом выборе х(0), 3, Б необходимо, чтобы матрица А удовлетворяла определенным условиям (см.

Теорема 1. Произвольные асинхронные итерации (А,/, х(0), 3, Б) сходятся к решению задачи (1) тогда и только тогда, когда

х(7) = |

x¿(7 — 1), если г € Jj, ЕП=1 хк(в к(7)) + /¿, если г € Jj,

(4)

[1, 2]).

А1(|А|) < 1,

(5)

где |А| — матрица, составленная из модулей элементов матрицы А, а Л1(|А|) — первое собственное число матрицы |А|.

Таким образом, если обычный (синхронный) процесс сходится— |Л1 (А) | < 1, но Л1 (|А|) > 1, то по крайней мере асинхронные итерации некоторого вида обязательно расходятся. Можно попытаться исправить положение, осуществляя после некоторой группы асинхронных итераций определенное количество синхронных, которые уменьшают ошибку (частичная синхронизация).

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

Оценка возможной получаемой выгоды существенно зависит от вида матрицы А ив общем случае может быть достаточно грубой. По-видимому, наиболее эффективным здесь может быть численный эксперимент. Тем не менее мы докажем лемму, которая указывает границы роста ошибки в асинхронном случае при Л1 (|А|) > 1.

Пусть х(7), 7 = 0,1, 2,... —последовательность асинхронных итераций, определенная набором (А, /, х(0), 3, Б), а X — решение системы (1). Тогда х(7) можно представить в виде х(7) = х + Дх(7), где Дх(7) —последовательность асинхронных итераций, заданная пятеркой (А, 0, Дх(0), 3, Б).

Лемма 1. Если для матрицы А выполнено |Л1(А)| < 1 и Л1(|А|) > 1, то

||Дх(к)|| < ||А||к||Дх(0)|| (6)

для к = 0,1, 2,...

Доказательство. Проведем доказательство по индукции. Для к = 0 имеем

|| Дх(0)|| = ||А||0|| Дх(0)||,

и, следовательно, база индукции доказана. Пусть теперь (6) выполнено для всех к < т. Покажем, что (6) выполняется при к = т. Согласно определению асинхронных итераций если 7 € Jm, то

п

Дxj(т) ^^ ajr Дхг(вг(т)),

Г=1

а если 7 € Jm, то

Дх^- (т) = Дх^- (т — 1).

Пусть Дх — вектор длины 2п, такой что Дх^- = Дх^- (т), Дхп+ = 0 при 7 € Jm и Дх^- =0, Дхп+ = Дх^- (т — 1) при 7 € Jm. Обозначим вектор (Дх1 (в1(т)),..., Дхп(вп(т)))Т за Дх. Для Дх, Дх и Дх(т — 1) справедливо равенство __

д. = (А1 0\/ Дх N Дж ^ 0 А^ уДх(т — 1У ,

где А1, А2 — матрицы п х п. При 7 € Jm 7-я строка матрицы А1 равняется 7-й строке матрицы А, а при 7 € Jm 7-я строка матрицы А1 состоит из нулей. А2 —диагональная матрица, для которой 7-й элемент диагонали равен 0, если 7 € Jm, и равен 1 в противном случае.

Заметим, что ||Дх|| = ||Дх(ш)||, ||А1| < ||А||, и в силу того, что ||А|| > |Лх(|А|)| > 1, выполнено неравенство || А21 < ||А||. Тогда справедливо неравенство

||Дх(т)|| =

А1 0

0

А2

Дх

Дх(т — 1)

< | А|

Дх

Дх(т — 1)

Пусть ||(Дх , Дх(т — 1)Т )|| = |Дх^(в')| для некоторого натурального в' < т и 1 < 3' < п. Так как в' < т, в силу предположения индукции будет выполнено

Дх

Дх(т — 1)

а следовательно

< ||Дх(в')|| < ||АГ ||Дх(0)| < ||А||т || Дх(0)||

|| Дх(т)|| < ||А||т|| Дх(0)||,

и лемма доказана.

Теперь легко доказать следующую теорему.

Теорема 2. Если итерационный процесс состоит из последовательности групп т асинхронных итераций и затем I синхронных, то при достаточно большом I он

сходится. При этом сходится не медленнее, чем Ае ^^^ + , где Ае удовлетворяет

неравенству |Лх(А)| < Ле < 1.

Доказательство. Из леммы 1 следует, что за т асинхронных итераций ошибка может возрасти не более чем в ||А||т раз. Поскольку |Лх(А)| < 1, для Уе > 0, такого, что е < 1 — |Лх(А)|, найдется такое 1о, что для VI > 1о будет выполняться неравенство ||Аг|| < (|Лх(А)| + е)1 < 1. Обозначим (|Лх(А)| + е) за Ле. Нетрудно видеть, что |Лх(А)| < Ле < 1. Норму ошибки после т + I итераций (т асинхронных и I синхронных) можно оценить следующим образом:

||Дх(т + 1)|| < ||Аг||||А|Г||Дх(0)|| < Л^||А|П|Дх(0)||.

При достаточно большом I имеем Л^ ||А||т < 1, что и доказывает первую часть теоремы.

Мы видим, что за т + I итераций ошибка уменьшается в Л^||А||т раз. Такой результат мы имели бы при геометрической сходимости с параметром г, если бы гш+1 = л^||А||т. То есть

г = (Ле|А|т)

1

№

что доказывает вторую часть теоремы.

2. Алгоритмы метода Монте-Карло. Известны различные алгоритмы метода Монте-Карло для решения системы вида (1). Мы будем рассматривать простейшие алгоритмы, состоящие в оценивании сходящегося ряда Неймана

(Н,*) = Н,^Ак /

V к=0 /

с помощью выбранной соответствующим образом оценки — оценки на траекториях моделируемой цепи Маркова. Здесь Н — заданный вектор, а х — решение системы (1). Как известно [3], при

Лх(|А|) < 1 (7)

Л

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

Как было показано в [3], нарушение условия (7) и использование асинхронного алгоритма приводит к стохастической неустойчивости — экспоненциальному росту дисперсии. Эта трудность, вообще говоря, может быть преодолена за счёт увеличения вычислительной работы, но её экспоненциальный рост делает алгоритм нереализуемым. Альтернативой является запоминание промежуточных результатов (синхронизация).

Предлагается в случае Лх(|А|) > 1 («дешевых» асинхронных алгоритмов и относительно «дорогих» синхронных), как и в случае обычных асинхронных итераций, использовать смешанный алгоритм с частичной синхронизацией. Как и ранее, будем рассматривать подход, в котором чередуются асинхронные и синхронные алгоритмы. Для формального исследования такого комбинированного алгоритма мы подробно рассмотрим схемы оценивания частичных сумм ряда Неймана в асинхронном и синхронном случаях.

2.1. Оценки частичных сумм ряда Неймана. Рассмотрим аналог схемы Неймана—Улама для нахождения частичной суммы ряда Неймана

т

к ,

Sm "У ] Ак /

fc=0

где m (Е N задано. Случай m = то подробно изучен [5]. Случай конечного m требует некоторой модификации стандартных рассуждений.

Пусть h — заданный вектор, h = (hi,..., hn)T. Алгоритм предполагает задание стохастической матрицы P = ||"j=i, удовлетворяющей условиям согласования:

(а) pi,j > 0, если ai,j = 0, и распределение q = (q1,..., qn) такое, что qi > 0, если hi = 0,

или

(б) pi j > 0, если aj i = 0, и распределение p = (p1,... ,pn) такое, что pi > 0, если / =0.

Также зададим вектор g = (gi,..., gm)T, такой что gj > 0, i = 1,..., m, и

m

Egj = 1.

¿=1

Введем случайную величину т, распределенную по закону

1 ... m

gi . . . gm

и рассмотрим случайные величины

^ _ hi(iai0,il ■•■Ч-1)'г/'т j-g-j

qi0 pi0 , ¿1 . . .piT-1 , ¿T gT

_ fjoajl,jo ■ ■ ■ ajT,jr-lh~jT ^

pjoPjo,jl • • -Рзт~1,3т 9т

где ¿0 ^ ¿i ^ • • • ^ iT соответствует цепи Маркова (q, P), а jo ^ ji ^ • • • ^ jT — цепи Маркова (p, P).

Легко проверить, что

m

EZ = (h,^Ak f ) = (h,Sm), k=i

m

EZ * = (£ fT (AT )k ,h) = (sm ,h). k=i

Столь же просто вычисляются дисперсии оценок (8), (9), если учесть, что

m n h2 Я2 ■ a2 • f2 /h2 m / A2 \ k f2 \

EZ2 = ^ ^ niouiQ,n ■ ■ ■ = / n_ (—) — j (10)

k=i io ,ii ,...ik=i qioPio'ii •••pik-i 3k \q Й=Д p7

hj0ljl jk = l PóoP3o,ñ---Pók-uñ9k \¿[ P \ Г i ' и )' где операции возведение в квадрат векторов и матриц, как и деление матрицы на матрицу и вектора на вектор, выполняются поэлементно. В дальнейшем будет отдельно оговариваться использование подобных операций над матрицами и векторами. Оценкой всех компонент Sm может быть вектор z = (Zi, • • •, Zn)T, где Z¿, ¿ = 1, • • •, n, получаются из (8), если полагать h равным вектору с i-й компонентой, равной единице, и остальными, равными нулю.

Как видно из (10) и (11), поведение дисперсии с ростом m определяется первыми собственными числами матриц A2/P и (AT)2/Р соответственно. Так, если выполнено Ai(A2/P) < 1 и Ai((At)2/Р) < 1, то соответствующие дисперсии будут ограничены при стремлении m к бесконечности. В противном случае будет наблюдаться экспоненциальный рост дисперсии.

Теперь для заданного вектора y = (yi, • • • , yn)T введем оценки п = n(y) и п* = П*(у) такие, что

En = (h,Am+iy),

En* = (yT (AT )m+i ,h)

Оценки n и n* строятся на траекториях цепей Маркова (p, Р) и (q, Р) длины m +1 и выражения для них имеют вид

hÍ0 ai0,il • • • ®im,im+1 Vi-m + 1 /10\

'П=-1-—, (12)

qi o pio

_ Vj0ajl,j0 ' ' ' ajm+l,jmh'jm+1 ^^

pjo pjo ,j1 • • •pjm,jm+1

Как и ранее, дисперсии оценок (12), (13) вычисляются просто, с учетом того, что выражения для вторых моментов имеют вид

m+i N

Evz = I —, ( ^ ) У2

где операции возведение в квадрат векторов и матриц, как и деление матрицы на матрицу и вектора на вектор, выполняются поэлементно. Поведение дисперсии, также как и для оценок (8), (9), с ростом т определяется первыми собственными числами матриц А2/Р и (Ат)2/Р.

Таким образом, при помощи пар оценок п, С и п*, С * для заданного вектора у и натурального т может быть получена оценка вектора Ат+1у + ^т=о /. Поручая моделирование цепей Маркова и вычисление оценок по ним разным процессорам, мы вновь получим асинхронный алгоритм.

2.2. Оценки суммы ряда Неймана. Оценки (8), (12) и (9), (13) могут служить основой алгоритмов с частичной синхронизацией для получения решения системы (1) без ограничений вида (7). Алгоритмы состоят в последовательном вычислении оценок векторов у (7), где у (7) связаны соотношением

т— 1

у(0) = /, у(7) = Ату(7 - 1) + £ А/, (14)

к=0

у(7 + 1) = Ау(7) + /,..., у(7 + 1) = Ау(7 + 1 - 1) + /, 7 = 1,1 + 2, 21 + 2,... (15)

или соотношением

т

у(0)= /, у(7)= Ат+1 у(7 - 1) + £ А/, 7 = 1, 2,... (16)

к=0

для некоторого фиксированного натурального т. Первый случай соответствует частичной синхронизации, когда аналогично описанному выше подходу для асинхронных итераций после асинхронного алгоритма следуют синхронные. Если 7 устремить к бесконечности, то у (7) в (14)—(16) будет стремиться к искомой сумме ^^=о /.

При использовании метода Монте-Карло оценками векторов у(7) будут случайные векторы 2(7) = (£1(7),..., £п(7))Т. Скалярное произведение (Л., у(7)) может оцениваться с помощью оценок вида (8), (12) —в этом случае оценкой для (Н, у(7)) будет

П(2(7 - 1))+ С + (Н,/),

или оценок вида (9), (13) —в этом случае (Н, у(7)) будет оцениваться как

П*(2(7 - 1))+ С * + (Н,/).

В каждом случае берется среднее из некоторого количества N реализаций оценок. Оценка всех компонент у (7) может быть получена, если полагать Н равным векторам с одной из компонент, равной 1, и остальными, равными 0.

При этом во многих случаях интерес представляет оценка у(7) для Т < 7 < (разностные схемы).

Заметим, что один и тот же вектор (£т в (16) и £т— 1 в (14)—(15)) оценивается на каждой итерации, поэтому можно организовать алгоритм таким образом, чтобы оценка этого вектора уточнялась с каждой новой итерацией, используя оценки с предыдущих шагов итераций.

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

Рассмотрим теперь итерационный процесс у (7) с начальным вектором у(0) и

у(7) = - 1)+ !, 7 = 1, 2,..., (17)

где В - матрица п х п, а ! — вектор длины п. Пусть при каждом 7 вычисление слагаемых в правой части (17) происходит с помощью рандомизированной процедуры. Таким образом, вместо последовательности у (7) возникает последовательность случайных векторов 2(7) = (£1 (7),..., £„(7))Т, 7 = 0,1, 2,..., с начальным вектором 2(0), связанных соотношением

2(7)= В,-2(7 - 1) + (18)

где В, = ||„&=1 —случайные матричные операторы, В(7) — случайные векторы. При этом для любого натурального 7 выполнено Е2(7) = у(7) = (у1(7),... , у„(7))Т, ЕВ, = В, ЕЮ(7) = !

Следующая лемма определяет характер поведения ковариации векторов ошибок Е(7)=2(7) - у(7).

Лемма 2. Пусть случайные операторы В,, векторы В(7) и 2(к) независимы в совокупности при любых 7 =0,1, 2,... и к < 7 в том смысле, что случайные величины а1, а2 «з, где а1 — произвольный элемент оператора В,, а2 — произвольный элемент В(7), а3 — произвольный элемент 2(к), независимы в совокупности. Тогда для матрицы ковариации вектора Е(7) справедливо соотношение

уессоуЕ(7) =В ® ВТуессоуЕ(7 - 1) + Е(Д, ® Д^)уессоуЕ(7 - 1)+

+ Е(Д, ® ДТ)уес(у(7 - 1)у(7 - 1)Т) + уеееоу^(7), (19)

где Д, = В, - В, ¿(7) = В(7) - !, уес — операция векторизации матрицы, а <8> — операция кронекеровского произведения матриц.

Доказательство. Подставим в (18) выражения для В,, В(7) и 2(7):

у(7) + Е(7) = (В + Д,)(у(7 - 1) + Е(7 - 1)) + ! + ¿(7).

Поскольку у(7) подчинен (17), получим выражения для Е(7) и Е(7)Т:

Е(7) = ВЕ(7 - 1) + Д,у(? - 1) + Д,Е(7 - 1) + ¿(7) (20)

и

Е(7)Т = Е(7 - 1)ТВТ + у(7 - 1)ТДТ + Е(7 - 1)ТДТ + ¿(¿)Т. (21)

Далее необходимо перемножить правые и левые части равенств (20), (21) и вычислить математическое ожидание от всех членов получившегося равенства. При этом стоит отметить, что математические ожидания многих слагаемых правой части будут равны нулю. Так, ЕД,Е(7 - 1)у(7 - 1)ТДТ = 0 в силу того, что Е(7 - 1) не зависит от Д, и ЕЕ(7) = 0. Учитывая эти соображения, имеем

соуЕ(7) = ВсоуЕ(7)ВТ + Е(Д,соуЕ(7 - 1)ДТ) + Е(Д,уС? - 1)уС? - 1)ТДТ) + coу¿(j).

Теперь, применяя операцию векторизации к матрицам соу£(_?'), соу£(_?' — 1), — 1)у(^ — 1)Т и ), получим (19), что и доказывает лемму.

Из леммы 2 можно вывести

Следствие 1. Для стохастической устойчивости алгоритма (18) необходимо и достаточно, чтобы модуль первого собственного числа матрицы В ( ВТ + Е(Дд ( ДТ) был меньше единицы.

Известно [6], что собственными числами В ( ВТ являются Л^(В)Л^(В), г, к = 1,..., п. Заметим также, что ||В ( ВТ|| = ||В|| ||ВТ||.

Для удобства введем следующие обозначения: Мд = Е(Дд ( ДТ), Вд = В ( ВТ + Мд и р = Е(Дд ( ДТ)уее(у(^ — 1)у(^' — 1)Т) + уессоу^'). Тогда выражение для ковариации из леммы 2 примет вид

уессоуЕ) = (В ( ВТ + Мд )уессоу£— 1) + Г>д.

Элементы матрицы Мд имеют порядок 1/Мд, где N —количество моделируемых траекторий на ^'-й итерации, то есть ЛЛ^ = где Л4' — матрица с элементами,

не зависящими от Жд.

Лемма 2 и следствие из нее позволяют сделать определенные выводы относительно стохастической устойчивости рандомизированных процедур для итерационных процессов (14)-(15) и (16).

Теорема 3. Если |Лх(А)| < 1 и на каждой итерации (Ц)-(15) вычисления осуществляются при помощи метода Монте-Карло, то для любого натурального т существует такое N', что для УЖд > N', ] = 1, 2,..., где Жд — количество моделируемых траекторий цепи Маркова для оценки у(^) в (Ц)-(15), рандомизированный итерационный процесс (Ц)-(15) будет стохастически устойчивым.

Доказательство. Рассмотрим (14) и заметим, что матрица В из леммы 2 равна Ат, при этом в явном виде матрица Ат не известна и оценивается при помощи метода Монте-Карло. Учитывая, что ||Вд|| < ||В|| ||ВТ|| + -^-ЦЛ^'Ц, норму ковариации ошибки после вычисления (14) методом Монте-Карло можно оценить как

||уессоу£(.Л|| < (\\Ат\\ \\(Ат)Т\\ + ^ессоу^ - 1)|| + ||2?,||.

Согласно (15) последующие I итераций ковариация ошибки будет изменяться в соответствии с

уес соу£(^' + 1) = ( А <8> Ат + -гр— М') уессоу£(^') + Рд+ь V Жд+1 )

Так как |Лх(А)| < 1, верно |Л1 (А ( АТ)| = Л2 (А) < 1 и существует такое N', что для УЖд+х > Ж' первое собственное число матрицы А <Е> Ат + ^у1 Л4' по модулю будет меньше единицы. Тогда, согласно следствию 1, рандомизированный итерационный процесс будет стохастически устойчивым.

Оптимальный выбор параметров т, I, N' из теоремы 3 зависит от матрицы А, вида оценок и свойств многопроцессорной системы. Как видно, при достаточно большом N' предложенная процедура соответствует случаю частичной синхронизации в детерминированном случае, за исключением того, что на каждом шаге итерации возникает случайная ошибка. Помимо очевидного увеличения числа моделируемых

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

Похожую теорему можно сформулировать для (16).

Теорема 4. Если |Ai(A)| < 1 и на каждой итерации (16) вычисления осуществляются при помощи метода Монте-Карло, то для любого натурального m существует такое N', что для VN,+1 > N', где N, —количество моделируемых траекторий цепи Маркова для оценки y(j) в (16), рандомизированный итерационный процесс (16) будет стохастически устойчивым.

Доказательство. При применении метода Монте-Карло для вычисления очередного y(j) в выражении (16) ковариация ошибки, согласно лемме 2, будет изменяться следующим образом:

veccov£(j) = (лт <g> (Ат)Т + jj-M'^ veccov£(j - 1) + Vj.

Сделанное предположение о матрице A, а именно |Ai(A)| < 1, означает, что выполнено неравенство |A1 (Am)| = |A1((Am)T)| = |A1 (A)|m < 1. Модуль первого собственного числа матрицы Am ® (Am)T равен |A1(A)|2m. Этот факт позволяет заключить, что для любого натурального m существует такое N', что для VNj+1 > N' модуль первого собственного числа матрицы Am Cg> (Ат)Т + ^у1 ЛЛ' будет меньше единицы. А это в свою очередь согласно следствию 1 означает, что процедура вычисления (16) методом Монте-Карло будет стохастически устойчивой.

При использовании (16) и метода Монте-Карло может быть построен стохастически устойчивый алгоритм. Он будет обладать большей асинхронностью, чем алгоритм, построенный на основе (14)-(15), но у такого алгоритма при равном числе моделируемых траекторий ковариация оценок будет больше.

3. Численные эксперименты. Для асинхронных итераций численные эксперименты проводились для случайным образом смоделированной матрицы A размерности 4 х 4, такой что A1(A) = 0.76, A1(|A|) = 1.64, ||A|| = 2.04.

На рисунке 1 изображена норма вектора ошибок для асинхронных итераций с множеством J = {Jjтаким, что k Е Jj для k = 1,2, 3,4, если j mod k = 0, а множество S такое, что sj(j) = j — 1, г = 1,..., n, j = 1, 2,... В данном случае наблюдается экспоненциальный рост ошибки.

На рисунке 2 изображена норма вектора ошибок для асинхронных итераций c частичной синхронизацией. Теорема 2 гарантирует сходимость нормы ошибки к нулю при m = 2, l = 6, однако сходимость будет наблюдаться уже при значениях m = 2, l = 3.

Для метода Монте-Карло численные эксперименты проводились для случайным образом смоделированной матрицы A размерности 10 х 10, такой что A1 (A) = 0.79, A1 (|A|) = 3.21, ||A|| = 8.45.

На рисунке 3 сплошной линией изображена норма вектора ошибок при решении (16) методом Монте-Карло с использованием оценок (8) и (12) для параметров

Итерации Итерации

Рис. 1. Рис. 2.

т = 2, N = 50000, 7 = 1,..., 30. Пунктирной линией изображены доверительные интервалы.

На рисунке 4 сплошной линией изображена норма вектора ошибок при решении (14)-(15) методом Монте-Карло использованием оценок (8) и (12) для параметров т = 2,1 = 4. При этом N = 50000 при расчете (14) и N = 5000 при расчете каждой итерации в (15). Пунктирной линией изображены доверительные интервалы.

Рис. 3. Рис. 4.

Заключение. Одновременное выполнение условий |Л1 (А)| < 1 и Л1(|А|) > 1 может создать определенные трудности при использовании как детерминированных, так и стохастических асинхронных методов для решения системы (1).

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

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

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

Литература

1. Chazan D., Miranker W. Chaotic relaxation // Linear Algebra and its Applications. 1969. Vol.2, N2. P. 199-222.

2. Baudet G. M. Asynchronous Iterative Methods for Multiprocessors // J. ACM. 1978. Vol.25, N2. P. 226-244.

3. Ермаков С. М. Метод Монте-Карло в вычислительной математике (Вводный курс). Невский диалект, Бином. Лаборатория знаний, 2009. С. 192.

4. Ермаков С. М. Метод Монте-Карло и асинхронные вычисления // Тезисы 1-й международной конференции общества Бернулли. Том 1. 1987. С. 462.

5. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. С. 472.

6. Ланкастер П. Теория матриц. Наука, 1973. С. 282.

Статья поступила в редакцию 26 июня 2014 г.

Сведения об авторах

Дмитриев Алексей Валерьевич — аспирант; alx.dmitriev@gmail.com Ермаков Сергей Михайлович — доктор физико-математических наук, профессор; sergej.ermakov@gmail.com

MONTE CARLO METHOD AND ASYNCHRONOUS ITERATIONS

Alexey V. Dmytriev, Sergey M. Ermakov

St.Petersburg State University, Universitetskaya nab., 7-9, St.Petersburg, 199034, Russian Federation; nkk@math.spbu.ru, nevolga@gmail.com

Asynchronous methods as well as Monte Carlo methods play an important role in using multiprocessor (multicore) computations and investigation of the connections between them and, especially, expansion the field of their applicability are important problems. System of the form x = Ax + f is considered, where x is unknown column-vector of length n, F is a column-vector of the right-hand sides of length n, A is n X n-matrix of the system. Case when the first eigenvalue of |A| less than 1 is well studied and in that case there are close connections between asynchronous iterations and Monte Carlo method. If the first eigenvalue of |A| is greater than 1, whereas the first eigenvalue of A less than 1, then similar methods appear to be inapplicable. New modifications of Monte Carlo method and asynchronous method with partial synchronization are proposed. Estimates of the errors and convergence conditions for proposed methods are obtained depending on the ratio are obtained. Numerical experiments are performed. Refs 6. Figs 4.

Keywords: Monte Carlo methods, asynchronous methods, asynchronous iterations, parallel algorithms, statistical modeling.