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

УДК 519.245

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

0 МЕТОДЕ МОНТЕ-КАРЛО

В СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ*

С. М. Ермаков1, А. В. Тросиненко2

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

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

2 АО «Технологии для авиации»,

Российская Федерация, 199178, Санкт-Петербург, Малый пр. В. О., 54, корп. 5, лит. П

Работа посвящена решению систем линейных алгебраических уравнений (СЛАУ) на компьютерах с распределенной памятью. Предполагается наличие М вычислительных узлов, каждый из которых имеет ограниченную быструю память, а обмен данными между узлами занимает значительное время.

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

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

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

В работе построены и исследованы оценки решения СЛАУ для рассматриваемого случая. При некоторых дополнительных условиях на матрицу системы доказаны достаточные условия сходимости. Библиогр. 8 назв.

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

Введение. Рассматривается задача численного решения больших систем линейных алгебраических уравнений (СЛАУ) на компьютерах с распределенной памятью. При использовании «наивного» подхода к организации параллельных вычислений каждый вычислительный узел должен иметь свою полную копию матрицы системы.

Предположим, что имеется т вычислительных узлов, каждый из которых имеет ограниченную быструю память, и обмен данными между узлами занимает значительное время. Если требуется решить «большую» систему линейных алгебраических уравнений методом итераций при условии, что элементы матрицы и вектора правой части невозможно разместить в памяти одного узла, то возникает проблема эффективного использования оборудования в промежутке времени между обменами. Одно из решений этой проблемы состоит в привлечении фоновых вычислений. Если этого не делать, возникает вопрос: может ли каждый из узлов использовать доступные ему данные в промежутке между обменами так, чтобы уменьшить общую невязку. В общем случае легко получить отрицательный ответ на этот вопрос, достаточно для т = 2 рассмотреть систему из двух уравнений (см. пример в приложении). Тем более

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

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

удивительно, что существует достаточно широкий класс систем, для которых это все же оказывается возможным.

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

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

X = АХ + ^ X=(хъ...,хп)т, Г=(Ь,...,и)Т, А=\К;, (1)

при условии р(|А|) < 1 применительно к вычислительной системе с распределенной памятью. Предполагается, что речь идет о достаточно большой системе уравнений, для решения которой может быть выгодным использование метода Монте-Карло. Рассматривается случай разделения матрицы между вычислительными узлами на непересекающиеся блоки строк при одинаковом разбиении на блоки индексов строк и столбцов. Далее рассматривается модификация метода простой итерации, основанная на этом разбиении, состоящая из двух вложенных итерационных процессов, такая, что только внешние итерации сопряжены с обменами сообщениями между узлами. Предлагается пример матрицы, для которой данная модификация приводит к расходящемуся процессу. Показывается, что для случая неотрицательной матрицы А сходимость не нарушается и, более того, скорость сходимости в терминах количества обменов информацией между узлами не ухудшается.

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

1. Преобразование системы. Как известно, метод простой итерации

Хк+1 = АХк + Г (2)

сходится к решению X системы (1) при любом выборе начального приближения Хо при условии р(А) < 1. На основе этого метода существует ряд оценок Монте-Карло для скалярного произведения (Н,Х) при условии р(|«4|) < 1 [5]. В частности, можно, беря в качестве Н векторы, один элемент каждого из которых равен единице, а все остальные — нулю, оценить все компоненты решения системы (1). Пусть индексы строк и столбцов матрицы А, а также векторов X и Г разбиты на блоки Б^,

i = l,...,m, что задает следующее разбиение матрицы и векторов:

Л -

(All . . . Aim

. , X =

\Am1 . . . Am

( X(1) \ /F(1)'

■ ) ' F =1

Vx (m)J \F(

m)

Если выполнять метод простой итерации для матрицы, разбитой между m вычислительными узлами на блоки строк, то после каждого умножения матрицы на вектор следует обмен частями посчитанного вектора. Поскольку хочется уменьшить количество обменов, была рассмотрена следующая модификация метода простой итерации в предположении, что выполняется условие p(Aii) < 1 при всех i.

1. Взять начальное приближение Yo = col (Yg^,..., Yo(m)) (вектор-столбец, составленный из Y)(i)).

2. Параллельно выполнять внешние итерации по к по схеме: на i-м из m узле

— собрать посчитанные на предыдущем шаге части вектора Yk;

— выбрать начальное приближение Yk+1 0;

— выполнять внутренние итерации по l вида

Y(i) = A--Y(i) + V A-Y (j) + F(i)-1k+1,l+1 = Aii Yk+1,l + Aij Yk + F ;

j=1,...,m j=i

1 -\s(i) Л/'ii)

— для некоторого lkii разослать значение Yk_1 ¡k . в качестве Yk+1.

3. Для некоторого к0 взять Yko в качестве ответа на задачу.

Заметим, что если во всех внутренних итерациях выбирать lk,i = 1 и Yk+1 o =

(i)

Yk ), получится обычный метод простой итерации. Если же все lk,i выбрать равными некоторому p, получится частный случай метода two-stage multisplitting. Если внутренние итерации (которые являются методом простой итерации для матрицы Aii) продолжать до сходимости, то при p(Aii) < 1 в их результате получатся (приближенные) решения уравнений

Yk+1 = AiYk+1 + Е Aij Yj + F(i), i = 1,...,m, (3)

j=1,...,m

j=i

(i)

относительно Их решения можно записать как

Y

k(__)1 = (I - Aii)-^Ai1lk(1) + ••• + Aiii-1Yfc(i-1)+

+ Ai,i+1 Yki+1) + ••• + AimY(m) + F(i)),

i = 1,...,m.

Обозначим через А^^ матрицу, диагональные блоки которой совпадают с соответствующими блоками А, а все остальные блоки — нулевые. В таком случае данный

итерационный процесс можно переписать в терминах метода multi-splitting [6], положив A = I - A, K = m и матрицы Bi одинаковыми и равными I — Adiag • При этом в качестве Di нужно взять матрицы, в l-м диагональном блоке которых стоит единичная матрица, а все остальные блоки — нулевые. Таким образом, в этом случае внешние итерации можно переписать в матричной форме:

Yk+i = BYk + F,

где

В

( _0

A22A21

A11A12 0

Аг\_А\т\

А22А2г,

F

\AmmAmi Amm Am2

0

/ ATiF^ \ A^pi 2)

1 "4 p(rn) ,

\AmmF /

(4)

(5)

положив Ац = (I — Ац)^1. В терминах [6] В — это матрица Н, а вектор Т равен Будем говорить, что В = Вз(А), где Б — разбиение множества индексов на т блоков (Б1,Б2,..., Бт). Рассмотрение точных решений уравнений (3) представляет интерес, поскольку, в отличие от метода простой итерации с конечным числом шагов, многие способы решения систем линейных алгебраических уравнений с использованием метода Монте-Карло дают несмещенную оценку решения, то есть, не смотря на наличие случайной ошибки, матожидание равно точному решению системы.

Утверждение 1. Если справедливы условия р(А) < 1, р(Ац) < 1 и р(В) < 1, итерационные процессы (2) и (4) сходятся к единственным решениям уравнений X = АХ + Г и У = ВУ + Т соответственно, и эти решения совпадают.

Доказательство. Очевидно, что решения будут единственными. Если бы это было не так, вектор разности двух решений одного и того же уравнения был бы собственным вектором (матрицы А или В соответственно), отвечающим собственному числу 1. Сходимость к решению процессов (2) и (4) следует из условия сходимости для метода простой итерации. Остается доказать, что решения совпадают. Это следует из единственности решений и того факта, что решение исходного уравнения X является стационарной точкой (4), поскольку каждый его компонент X^ является стационарной точкой внутренней итерации для г-го блока строк. ■

Данное утверждение говорит о том, что применение рассмотренной модификации итерационного процесса при выполнении определенных условий не изменяет решение. Следующая теорема проясняет условие р(В) < 1.

Теорема 1. Если все элементы матрицы А неотрицательны, выполняется следующее условие: если р(А) < 1, то р(В) ^ р(А) < 1.

Доказательство. Для доказательства этой теоремы можно сослаться на лемму из раздела 3 статьи [7]. Действительно, (I — Аа^)-1 = I + Ашав + А«^ + ••• > 0 (здесь и далее подразумеваются поэлементные неравенства), поскольку А^^ ^ 0 и матричный ряд сходится, поскольку р(Амад) < 1, так как 0 ^ А^ад ^ А. Также поскольку А — А<цаг ^ 0, то В = (I — Аа^)~1(А — А<тв) ^ 0. Значит, пара матриц (I — Аатв, А — А<лав) является слабым регулярным разбиением матрицы I — А. Очевидно, что (I, А) также является слабым регулярным разбиением I — А. Поэтому, так как ^—АГ1 = !+А+А2 + ••• > 0 и А—А<^ < А, то р((!—А^)-1(А—А<-ше)) < р^-1 А), то есть р(В) ^ р(А), что и требовалось доказать. ■

Следствие 1. Для произвольной матрицы А выполняется следующее свойство: если р(|А|) < 1, то р(|В|) ^ р(|А|) < 1.

Доказательство. Во-первых, заметим, что при всех г выполняется условие Р^Лц|) < 1, а значит, матричный ряд I+|AJ^г|+1|2+... сходится. Пусть В' — Вз(|А|), тогда все блоки матриц В и В' можно выразить одинаковым образом через блоки матриц А и |А| соответственно, с помощью операций матричного сложения и умножения. А поскольку для любых матриц А и В поэлементно выполняются неравенства А + B| < А + |B| и АВ| < |А||В|, то поэлементно выполняется 0 < В < В', следовательно будем иметь р(|В|) < р(В') < р(|А|) < 1. ■ Заметим, что, вообще говоря, из сходимости процесса — АХ+ Р не следует сходимость процесса Уи+1 — ВУь + Т .В качестве примера можно привести матрицу, рассмотренную Холтоном:

а—(: -

При т — 2 и 0 < : < 1/2 прямыми вычислениями получаем

/О

В=( 'п" у: 0

р(А) — :, р(В) —

л/1 - 2а'

При а, стремящемся к 1/2 слева, получаем р(А) < 1, а р(В) неограниченно растет.

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

1. Взять начальное приближение Y0 = col (Y0(1),..., У"0(т)).

2. Параллельно выполнять внешние итерации по к следующего вида: на i-м из m узле

— собрать посчитанные на предыдущем шаге части вектора Yk;

— каким-либо образом вычислить приближенное решение (3):

потребовав = 6к,г < ~к,г, где через

обозначен i-й блок строк

матрицы В.

3. Для некоторого к0 взять Yko в качестве ответа на задачу.

При этом потребуем выполнения — = 6f¡¿ ^ £k}i- Поскольку для Y¡~, Y¡~+1

и X выполнены равенства

Yk+i = (I - Ai)-1 (a^Y^ + ••• + Ai,i-iYj:i-11+

+ Aii+iYÜ+1) + ••• + AimYkm + F(i)), i = 1,...,m,

= (I - АиГ1 (АаХ{1) + • • • + 1} +

+ Аг,г+1Х(г+1) + ---+АгтХ(т) + Г«

1,...,т,

то, вычитая одно из другого, получаем

3 = 1,...,т

3=г

Таким образом, если рассмотреть матричную норму, согласованную с рассматриваемой векторной, и ввести матрицу I = {¿з }™з=1 такую, что справедливо

= Г П (I — Л^-1^ у, г = 3,

[0, г = 3,

то ошибки на (к + 1)-м шаге можно выразить через ошибки на к-м шаге следующим образом:

У к+1 — Х

<

¿гз

X ¿гз

П = 1,...,т з=г

(6)

У к+1 — У

(г)

к+1

и можно положить £к+1 = Век + где через и обозначены вектора размера то, составленные из компонентов и соответственно.

Пусть С — матрица размера т х т, составленная из норм блоков Агз матрицы А: сгз = \\Aij у. Тогда, если все ее диагональные элементы меньше 1, т.е. сгг < 1, элементы матрицы I можно оценить следующим образом:

сгз/(1 — сгг ), г = 3,

г = 3,

так как (I — Аи)-1 = I + Аи + А2г + ..., \\(I — Аи)-1\\ < 1/(1 — си). Поскольку I = В1,..,т(С) (тривиальное разбиение индексов на блоки: каждый индекс — это отдельный блок), то по теореме 1, если р(С) < 1, то р(1) < 1, а значит, р(1) < 1, так как поэлементно 0 ^ I ^ О.

Особый интерес представляет вычисление приближенных решений методом Монте-Карло. В таком случае будем считать вк^ и в к,г случайными величинами, а £кгг и £к,г —детерминированными. При этом указанные неравенства для погрешностей следует переписать как Евк^ ^ £к,г и Е&к,г ^ ~к,г-

2. Асинхронный подход. Учитывая выведенные неравенства на ошибки, можно уточнить рассмотренную выше схему, использующую метод Монте-Карло вместо внутренних итераций, описав требуемые погрешности ~к,г, необходимые для достижения заданной точности. Однако в этом случае алгоритм будет состоять из чередующихся фаз вычислений и обменов по сети. Поэтому интерес представляет адаптация метода асинхронных итераций [3]. Заменим фазы вычислений и обменов на постоянные вычисления с асинхронной рассылкой текущих посчитанных приближений и

0

использованием информации о чужих приближениях при получении такой информации. Однако, в отличие от [3], будем вычислять сразу блоки неизвестных (т. е. вместо n неизвестных рассмотрим m блоков неизвестных). Также, в отличие от существующих работ (например, [8]), совмещающих метод multi-splitting и асинхронные итерации, вычисления будем производить приближенно с помощью метода Монте-Карло. Этот подход тем не менее требует дополнительного анализа, поскольку, в отличие от метода асинхронных итераций, компоненты решения вычисляются со случайной погрешностью. В частности, нельзя просто отбрасывать результаты предыдущих вычислений при приходе обновлений. Логично предположить, что если предыдущая оценка была посчитана по методу Монте-Карло с объемом выборки, скажем, миллион, то после получения небольшого уточнения «чужой» части вектора решения в момент, когда объем новой выборки равен всего лишь тысячи, старая оценка за счет меньшей дисперсии может быть точнее, несмотря на то что у новой смещение относительно решения исходной задачи будет несколько меньше. Поэтому, в отличие от обычных асинхронных итераций, будем строить оценку по нескольким предыдущим приближениям.

Зафиксируем номер узла i. Пусть Yi,Y2,...,Yk — последовательно приходившие на этот узел уточненные приближения решения исходной задачи X (например, два соседних приближения могут отличаться только в одном блоке индексов, если изменение вызвано обновленным приближением части решения, посчитанным одним из узлов), £\,£2, ■■■,£к —векторы размера m такие, что E6i ^ £i. Для каждого Yi методом Монте-Карло вычислена несмещенная оценка Z(i) = ■ ■ ■ ,zi,i,ni)т вектора

yj^ = B^Yi + J7^) = (уг+1 j 1;..., V[+i i щ)Т■ Для простоты пренебрежем асимптотическим характером нормальности оценок и будем считать, что zi^- имеет распределение No'fi j/Mij), где через М;^ обозначен объем выборки для метода Монте-Карло, и при различных I случайные величины z^ij —ynj независимы. Требуется с учетом значений £i, of,- и Mi ц по оценкам Zf^ получить в каком-то смысле

' ' (i) (i) (i) оптимальную оценку. Будем искать ее в виде Z() = f31Z() + fi2Z2 ) + ••• + в к Z(),

где ei — некоторые неотрицательные числа, такие что в\ + в2 + ••• + вк = 1. Затем

вектор Zбудет разослан i-м узлом в качестве уточненного приближения i-го блока

решения.

Рассмотрим величину —Х^ ^ Ц2 как меру качества приближения и устремим

к минимуму ее оценку сверху:

\i=i

\i=i

J2 A (zf" - в«>и - -F»)) + ((е«У, + - х">)

\

<

<

к а2

Е/^ + Е/^ь (7)

i=1

i=1

где а"2 i = Yln=! а2 i -, а через D (i) обозначена i-я строка какой-либо поэлементной

■yUi

-=1 ai, i, j

E

2

E

2

оценки сверху матрицы D, поскольку

к

Е

1=1

l3i (z^ - B(i)Y - F(i))

к a2-

Для задачи минимизации с данной оценкой выполняются некоторые естественные свойства: пусть имеем Ylo = Ylo+1 (а также £lo = £lo+1)- Зафиксируем все в при l = ¡o,lo + 1 и найдем минимум по /3lo, ^l0+i (fil0 + A0+i = с = const). Вторая сумма в правой части выражения (7) при этом остается постоянной, так как £l0 = £l0+i. Поэтому, не умаляя общности, можно считать, что мы решаем задачу минимизации

в2

Ml м2

при условиях в1 + в2 = 1, вг ^ 0, поскольку аг2о = стг2 +1 г, при этом в¡0 = в1с, в10+1 = в2с. Найдя производную по в1:

Pi , (i-A)2

Ml м2

Л

2 в1

получаем, что минимум достигается при

M1

в1

Mi + M2''

в2

1 1

мТ + М2

M2

1

М2

M1 + M2

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

— Х^ ^ 112 выборки для ^^ и будут объединены.

Пусть теперь для приближений У1о0 и Уь0+1 поэлементно выполняется неравенство

и а10+1 г при вычислении Zl0 и

£l0 ^ £l0+i. Также предположим, что величины af

Zl0+l совпадают. Пусть при этом оптимальной оценкой •)((г1 является ^к=1 вl,lZ'^ 1 Предположим теперь, что уточнение У10+1 было отброшено, а остальные уточнения начинали использоваться в те же моменты времени, что и раньше. По показанному выше можно считать, что вместо У0+1 было использовано У10 (если пренебречь накладными расходами на начало обработки нового приближения). Пусть в этом случае оптимальной оценкой •)((г2 оказалась линейная комбинация ^г= ,2 Zl( 2. Тогда ма-

(i)

тожидание нормы оценки ошибки •)((г1 будет не больше, чем у ^к=1 вl,2ZЦ (в силу

оптимальности выбора вм). В этом случае при переходе к •)(г2 = к=1 в^2Zl('2) первая сумма в правой части формулы (7) останется прежней, а вторая не уменьшится. Таким образом, этот способ нахождения оптимальной оценки обладает также тем свойством, что при некоторых предположениях приход улучшенного приближения решения в любой момент не ухудшает результат.

Формализуем теперь, аналогично описанному в [3], требования к обменам приближениями. Пусть Уц, Уг2 ,Угз,... —приближения, получаемые г-м узлом, при этом все Уг1 совпадают (то есть Уг1 = У1, ег1 = £1 при всех г). Обозначим через Т (Уц) момент времени, когда 3 -й узел посчитал свой компонент приближения, впоследствии

(з)

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

(i)

E

2

при ¡1 < ¡2 частным случаем линейной комбинации ¡2 приближений является комбинация, в которой вь отличны от 0 только для I ^ ¡1). Потребуем отсутствие переупорядочивания сообщений между любыми двумя узлами (то есть если ¡1 < ¡2, то Т (Уц1) < Tj(Уц2)). Тогда покомпонентно выполняется £ц1 ^ £ц2 при ¡1 < ¡2. В [3] было доказано достаточное условие сходимости асинхронных итераций, однако в данном случае дело осложняется наличием случайных ошибок при вычислении компонент приближения, поэтому требуется дополнительное рассмотрение. В следующей теореме доказывается достаточное условие сходимости, при этом формулируются требования, аналогичные предъявляемым в [3] для метода асинхронной итерации.

Теорема 2. Пусть Tj(Уц) монотонно стремится к бесконечности как функция от ¡ для любых г и ]; для любого г сумма^к=1 Мь,г стремится к бесконечности как функция от к; величины а2 г равномерно ограничены сверху значением а2, а также £ц1 ^ £ц2 при ¡1 < ¡2. Тогда если ЦБ||2 < 1, то на любом узле г последовательность приближений Уц сходится в ¡2-норме к решению исходного уравнения.

Доказательство. Не умаляя общности, можно считать, что все компоненты £1 положительны. Докажем по индукции по схеме, аналогичной доказательству теоремы 1 в [3], что для любого к существует такое ¡к, что £цк ^ Шк£1 при всех г, где

Ш

2

(здесь и далее под || • || понимается ¡2-норма). При к = 0 можно взять ¡о = 1. Пусть для к это утверждение уже доказано, докажем для к +1. Рассмотрим момент Ък+1, когда для всех г выполняется неравенство ^1=к+к Мь,г ^ Мк, где величина Мк удовлетворяет условию

• 1 ^

1 -ЦБ ||

Ш к £1

Так как справедливо неравенство Шк£1 ^ ((1 — ЦП11)/2)к£1, все компоненты вектора в правой части положительны, следовательно, при достаточно большом Мк неравенство выполняется. Зафиксируем номер узла г. Тогда, даже если взять вь =0 при ¡ < ¡к и вь = М,г/^Mj,г при ¡к < ¡ < Ьк+1, то будем иметь

Ьк + 1 2 Ьк+1 / Л/Г ~ Ьк +1

'1 Мц ^ V £ щ*)2 (-1

1 — 1к 2

Мк

Мк

Значит, в момент ¡к+1, когда при всех г и ] величина Т (Уцк+1) будет больше момента вычисления

, при любом г можно положить

£ик+1 = ВШк£Х + А • 1 < мгк+1£1.

2

3. Случай разбиения с перекрытием. Как показано в [1], также представляет интерес вариант использования пересекающихся блоков, поскольку такой подход может ускорить сходимость к ответу. Для анализа случая разбиения с перекрытием в [1] используется метод multi-splitting. Рассматривается набор троек (Mi,Ni,Ei) такой, что A = Mi — Ni, и Mi —невырожденная при всех i = 1,...,m, и кроме того, Yim=i Ei = I, причем все Ei —неотрицательные диагональные матрицы. Рассматривается итерационный процесс Хк+1 = ^™1 EiYk,i, MiYk,i = NiXk + F. Ранее рассмотренный итерационный процесс Yk+i = BYk + F, сходящийся к решению системы уравнений (I — A)Y = F, как упоминалось выше, описывается тройками матриц таких, что при всех i на диагонали матрицы Ei стоит либо 0, либо 1. Один из приведенных в [1] способов достичь перекрытия блоков заключается в том, чтобы матрицы Ei оставались блочно-диагональными матрицами с единичным блоком, соответствующим индексам из множества Si, и остальными нулевыми, а при всех i матрица Mi «вырезала» из I — A диагональ и i-й квадратный блок, соответствующий всем индексам из Si и еще некоторым, обеспечивающим перекрытие. Более точно будем рассматривать тройки (Mi, N, Ei), Mi = {(Mi)uv}, Ni = Mi — A, Ei = {(Ei)uv}, где

и для подмножеств Ti индексов матрицы при всех i выполняется условие Si С Ti. Как указывалось выше, матрица B является матрицей итерации метода multi-splitting при всех Mi, равных матрице I — A с обнуленными внедиагональными блоками (в соответствии с разбиением {Si}"=i) и матрицами Ei, описанными в (8). В [1] отмечается, что такую же матрицу B можно получить, рассмотрев метод multi-splitting с тройками (Mi,Ni,Ei), если взять Ti = Si при всех i. Это разбиение в общем случае Si С Ti легко использовать в рассмотренном выше алгоритме, взяв вместо Aii квадратные подматрицы A^ матрицы A, соответствующие индексам Ti (Aii соответствуют Si). При этом на каждом узле будет храниться весь блок строк матрицы A, соответствующий множеству индексов Ti. Заметим, что i-й узел все равно должен вычислять только компоненты вектора приближения, соответствующие индексам из Si. Это в случае использования метода Монте-Карло означает, что для компонентов из Ti\Si оценки строить не требуется (в силу того, что каждый компонент оценивается отдельно). При этом вместо матрицы B следует рассматривать B = ™i EiM-1 Ni.

Из результатов, изложенных в [1], очевидно выводится следующее утверждение.

Утверждение 2. Если матрица A неотрицательная и р(А) < 1, то справедливо неравенство p(B) ^ p(B).

Доказательство. Очевидно, что матрица (I — A) является M-матрицей, поэтому требуемое утверждение следует из теоремы 2 [1]. ■

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

если u € Ti и v € Ti, ( если u = v,

0, иначе,

1, если u = v € Si,

(8)

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

когда имеется два процесса. Во время между т-м и (т +1)-обменами первый процесс имеет в памяти о,1 и Ь1ут-1 + fl, а второй соответственно 62 и а,2Хт-1 + /2 (хо и уо заданы). За это время находятся хт и ут из условий

Если |ai| < 1 и |&21 < 1, то каждое из уравнений можно решать методом итераций. Если мы успеем сделать достаточно много этих итераций, можно считать, что выполняется равенство

Для m = 1, 2,... имеем итерационный процесс, который может сходиться при lai| + |bi| < 1 и |«21 + |&21 < 1, но может и расходиться, например, при ai=&2=0, 5; bi=1,1; «2=—1,1. Более общие случаи исследуются достаточно сложно, кое-что в этом направлении сделано в данной статье.

Литература

1. Frommer A., Pohl B. A comparison result for multisplittings based on overlapping blocks and its application to waveform relaxation methods. Zurich: ETH, 1993.

2. Szyld D. B, Jones M. T. Two-stage and multisplitting methods for the parallel solution of linear systems // SIAM J. Matrix Anal. 1992. Vol. 13, issue 2. P. 671-679.

3. Bo,ud,et G. M. Asynchronous Iterative Methods for Multiprocessors // Journal of the Association for Computing Machinery. 1978. Vol.25, N2. P. 226-244. URL: http://doi.acm.org/10.1145/322063.322067 (дата обращения: 30.08.2016).

4. Strikwerdo J. C. A probabilistic analysis of asynchronous iteration // Linear Algebra and its Applications. 2002. Vol.349, issues 1-3. P. 125-154.

5. Ермаков С.М. Метод Монте-Карло в вычислительной математике (вводный курс). Санкт-Петербург: Невский диалект, Бином, 2009.

6. O'Leory D. P., White R. E. Multi-splittings of matrices and parallel solution of linear systems // SIAM Journal on Algebraic and Discrete Methods. 1985. Vol.6, N4.

7. Elsner L. Comparisons of weak regular splittings and multisplitting methods // Numerische Mathematik. 1989. Vol. 56. P. 283-289.

8. Bru R., Eisner L., Neumann M. Models of parallel chaotic iteration methods // Linear Algebra and its Applications. 1988. Vol.103. P. 175-192. URL: http://www.sciencedirect.com/science/article/pii /0024379588902273 (дата обращения: 30.08.2016).

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

в [1].

ПРИЛОЖЕНИЕ

Поясним ситуацию на простейшем случае двух уравнений

x = a\x + biy + fi, y = a2x + b2y + ¡2,

Xm = aixm + biym-i + ¡1, Ут = a2xm-i + Ь2Ут + f2■

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

Ермаков Сергей Михайлович — доктор физико-математических наук, профессор; [email protected]

Тросиненко Анатолий Владимирович — инженер-программист; [email protected]

ON MONTE CARLO METHOD IN DISTRIBUTED MEMORY SYSTEMS

Sergey M. Ermakov1, Anatolii V. Trosinenko2

1 St. Petersburg State University, Universitetskaya nab., 7—9, St. Petersburg, 199034, Russian Federation; [email protected]

2 JSC "Technologies for aviation", Maliy pr. V.O., 54—4, St. Petersburg, 199178, Russian Federation; [email protected]

The work is devoted to solving systems of linear algebraic equations on distributed memory computers. It is supposed that there are M computing nodes, each of which has a limited fast memory, and communication between nodes takes considerable time.

Provided that the matrix elements and vector of the right side can not be placed in its entirety in the memory of one node, the problem of effective use of equipment in between exchanges, i.e. whether each node to use the data available to him to be able to reduce the overall residual. In general assumptions about the matrix of the answer to this question is negative. An example is placed in the Appendix.

We consider the case when the system has a large enough order and appropriate to apply the Monte Carlo method, the matrix is divided between computing nodes on non-overlapping blocks of rows with the same partition into blocks of indices of rows and columns. We also consider a modification of the method of simple iteration, based on this partition, consisting of two nested iterative process, so that only the outer iteration involve the exchange of messages between nodes.

This iterative process leads naturally to a similar process using the Monte Carlo method that does not require a full copy of the storage matrix of the system on each compute node.

The paper constructed and investigated for solving the linear estimates for the case and under certain additional conditions on the matrix system proved sufficient conditions for convergence. Refs 8. Keywords: Monte Carlo method, parallel computation, asynchronous iterative method.

References

1. Frommer A., Pohl B., A comparison result for multisplittings based on overlapping blocks and its application to waveform relaxation methods (ETH, Zürich, 1993).

2. Szyld D.B., Jones M.T., "Two-stage and multisplitting methods for the parallel solution of linear systems", SIAM J. Matrix Anal. 13, Issue 2, 671-679 (1992).

3. Baudet G.M., "Asynchronous Iterative Methods for Multiprocessors", Journal of the Association for Computing Machinery April. 25(2), 226-244 (1978). Available at: http://doi.acm.org/10.1145 /322063.322067 (accessed 30.08.2016).

4. Strikwerda J.C., "A probabilistic analysis of asynchronous iteration", Linear Algebra and its Applications 349, Issues 1-3, 125-154 (2002).

5. Ermakov E., Monte Carlo method in computational mathematics (introductory course) (Saint Petersburg, Nevskiy dialekt, Binom, 2009) [in Russian].

6. O'Leary D.P., White R. E., "Multi-splittings of matrices and parallel solution of linear systems", SIAM Journal on Algebraic and Discrete Methods 6(4) (1985).

7. Elsner L., "Comparisons of weak regular splittings and multisplitting methods", Numerische Mathematik 56, 283-289 (1989).

8. Bru R., Elsner L., Neumann M., "Models of parallel chaotic iteration methods", Linear Algebra and its Applications 103, 175-192 (1988). Available at: http://www.sciencedirect.com/science/article/pii /0024379588902273 (accessed 30.08.2016).

Для цитирования: Ермаков С. М., Тросиненко А. В. О методе Монте-Карло в системах с распределенной памятью // Вестник Санкт-Петербургского университета. Серия 1. Математика. Механика. Астрономия. 2016. Т.3(61). Вып. 4. С. 558-569. DOI: 10.21638/11701/spbu01.2016.405

For citation: Ermakov S. M., Trosinenko A. V. On Monte Carlo method in distributed memory systems. Vestnik of Saint Petersburg University. Series 1. Mathematics. Mechanics. Astronomy, 2016, vol. 3(61), issue 4, pp. 558-569. DOI: 10.21638/11701/spbu01.2016.405