Алгоритм Монте-Карло для решения систем линейных алгебраических уравнений методом Зейделя∗ Текст научной статьи по специальности «Математика»

УДК 519.2

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

АЛГОРИТМ МОНТЕ-КАРЛО

ДЛЯ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ

АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ МЕТОДОМ ЗЕЙДЕЛЯ*

Т. М. Товстик, К. С. Волосенко

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

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

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

Ключевые слова: алгоритм Монте-Карло, метод Зейделя, система линейных алгебраических уравнений.

1. Введение. Метод Монте-Карло для решения системы линейных алгебраических уравнений рекомендуется использовать в том случае, если ее порядок достаточно велик [1]. В монографиях [1, 2] с помощью метода Монте-Карло вычисляется одна компонента вектора решения или скалярное произведение вектора решения и произвольного вектора. В работах [3-5] оценивается сразу весь вектор решений. В работе [3] используется стохастический метод итераций, причем математическое ожидание случайного вектора очередной итерации совпадает с суммой соответствующего отрезка ряда Неймана. В статьях [4, 5] представлены алгоритмы Монте-Карло, позволяющие получать итерационное решение системы линейных уравнений при условиях более слабых, чем условия обычного метода Монте-Карло. В данной работе, в отличие от алгоритма в [3], особенность построения случайного вектора заключается в том, что математическое ожидание очередной итерации совпадает с итерацией Зейделя [6]. В дополнение к работе [7] доказывается существование и конечность предельных дисперсий случайного вектора решений системы.

Пусть задана система линейных алгебраических уравнений вида

X = A ■ X + f, (1)

где А = — квадратная матрица, а X = (Хь ..., Хп)т и / = (/ь ..., fn)T —

векторы.

В работе рассматривается алгоритм Монте-Карло, в основе которого лежит метод Зейделя. Будем предполагать, что для нормы матрицы A выполняется условие

n

IIAll = max ]Т\Aik| < 1, (2)

Кi<n *—' --k=1

которое, в частности, обеспечивает сходимость метода Зейделя к единственному решению системы (1) при любом выборе начального вектора [6].

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

Пусть задан начальный вектор X(0) = /, тогда на то-й итерации по методу Зей-деля [6] компоненты вектора X(т) вычисляются по формулам

г-1

= £ х(т) А. х(т-1) + /г, г = 1,...,п, то > 1.

(3)

3 = 1

2. Алгоритм Монте-Карло для метода Зейделя. Рассмотрим метод, который имитирует алгоритм Зейделя. Будем моделировать случайные векторы то > 0 такие, что выполняется £(0) = /, а при то > 1 справедливо Е £(т) = X(т), причем компоненты X(т) удовлетворяют равенствам (3). В основе моделирования лежит имитация цепи Маркова, которая задается вектором начального распределения п и стохастической матрицей Р вероятностей перехода:

\\Pij||, Р. > ° 1 < г,] < П

3=1

1 < г < п.

(4)

Матрицу Р будем считать допустимой по отношению к матрице А, если выполняется р. > 0 для тех пар (г,]), для которых справедливо Аг. = 0. На матрицу Р не накладывается никаких ограничений. В частности, она может состоять из отдельных блоков.

Моделирование вектора С(т) осуществляется следующим способом. Выбирается начальное приближение С(0) = /, далее при то > 1 компоненты случайного вектора обновляются последовательно, начиная с г = 1, причем при моделировании г-й компоненты учитывается результат обновления всех предыдущих компонент ((т\ ] < г. При каждом г согласно вероятностям (4) разыгрываем номер состояния ]г как дискретную случайную величину с распределением

1 2

Рг1 Рг2

п

ргп

Тогда вектор С(т) можно записать в виде

с(т) = • (гп-1 • (... • • (г1 + /1) + /2) +...

1)

(5)

причем

Zi

1

\

0

\

0

1

.0

0

/г 0

.0

г = 1,2,... ,п.

Матрица в каждой строке имеет только один ненулевой элемент, при этом в г-й строке на ]г-м месте стоит элемент А. /р., а в остальных строках на диагоналях стоят единицы.

1

0

1

Итак, имеем = /, г = 1, 2,...,п, компоненты вектора С(т) при т > 1 моделируются последовательно, начиная от г = 1 до п по формулам

(т) , 1 ^АЪ /РЪ пРи ^ < ^

СГ = fi + { :m ' (6)

Qi 'Aiji / Piji при ji > г,

где случайные состояния ^ разыгрываются в соответствии с распределением (5). Формулы (6) допускают представление

i—1 A' ' / \ n А- ■

Ai,j ..im) X ^ Aij,i

Лг) Ai,ji Am) i Ai,j Am— 1) , x 1 г»

Q =2^X(k,ji)-j2-Cji '+2^X(k,ji)-j2-Cji > + fi, г = 1,2,...,п, (7)

k=1 Pij k=i Pij

где X(k,p) — символ Кронекера.

Вычисление математического ожидания от обеих частей равенства (7) приводит к равенствам

i— 1 n

ЕСГ = £ Aij EZ(m) + ^ ECjm—1) + fi, г = 1, 2,...,n, m > 1. (8) j=1 j=i

Отсюда при m =1 получаем

n i— 1 n

EZ(1) = £ A1j fj + f1, EZ(1) = £ Aij Ej + E Aij fj + fi, г = 2,...,n. (9)

j=1 j=1 j=i

Из формул (9) следуют равенства EZ(1) = X(1), EZ(1) = X(1), в которых X(1) те же, что и в (3) метода Зейделя, если начальный вектор удовлетворяет X(0) = f. Нетрудно проверить, что при всех m > 1 выполнены равенства EZ(m) = X\m\ Теорема 1. Пусть для системы линейных алгебраических уравнений вида

X = A ■ X + f

выполнено условие

n

\\A\\ = max ]Т\Aik\ < 1,

1<i<n ^—' --k=1

и матрица P является допустимой для A. Пусть случайный вектор Z(m) имеет компоненты

i—1 A n A

ä0) = ь, i < < < n, cf° = Е^^сГ+Е^^г1^/, rn>i.

k=1 Piji k=i Piji

Тогда справедливы равенства

i— 1 n

ECr = ХГ = £ AijX(m) + ^ Aijx(m—1) + fi, 1 < г < n, (10) j=1 j=i

lim EZ(m) = Xi, 1 < i < n. (11)

Доказательство. Равенство (10) доказано выше. В [6] для то-го приближения Х(т) метода Зейделя доказано \Хг — Х(т^\ —> 0 при то ^ ж для всех г. Так как имеет место равенство Е(г(т) = Х(т), выполняется \Хг — Е^™^ —> 0 при то ^ ж, тем самым теорема доказана.

Для оценки решения уравнения (1) следует достаточно большое число N

(м)

N

раз моделировать случайный вектор С(М), тогда будем иметь ес(м) «< С(М)(S)/N, где в — номер серии.

3. Оценка числа итераций. Чтобы оценить номер итерации М, на котором следует остановить моделирование для получения нужной точности, воспользуемся соответствующими неравенствами, доказанными в [6] для метода Зейделя при норме матрицы, имеющей вид (2). Для компонент вектора С(М) при всех г справедливы неравенства

м

(12)

1 — ^

где ц и 5 определяются выражениями

V" 1А--1

ц = тах < ||А|| < 1, <5 = тах - Е^>|.

4 1-Е}:! \м <

Пусть о"г(М) = 1 < % < п, и <т(м) = тах^ стДМ), тогда согласно центральной

предельной теореме с вероятностью 7 близкой к единице при любом г выполняется

соотношение

1 N

в=1

г а(м)

< (13)

Например, если 7 = 0.95, имеем = 1.96.

Неравенства (12) и (13) служат рекомендацией для оценки величины М. Принимая

а(м) Н

с той же вероятностью 7 получаем оценку

1, (14)

1 N

в=1

а(м )(г7 + Н) 1 .

Величину М следует находить из соотношения (14). Величина а(м) неизвестна, но если нормы матрицы А и вектора / меньше единицы, будем иметь а(м) « 0(1).

Моделирование вектора £(т) дает возможность найти аг(М) —выборочные оценки среднеквадратических отклонений аг (М) и построить доверительные интервалы для величин Хг.

4. Корреляции предельного вектора. Введем вектор X(т) с компонентами

2(т) = п((Ш) = Е М-А2 , 1 < г < п, то > 0.

Теорема 2. Пусть выполнены условия теоремы 1 и норма матрицы В = [В. ] = [А. /ргз] удовлетворяет неравенству

п

\\В\\ = та? £Вш < 1,

1<г<п --к=1

тогда существуют конечные пределы

Ег] = lim E (c(m)c(m)) , Кг] = lim E (dC(m-1)) , 1 < г,3 < n, (15)

т^ж \ J / т^ж \ J /

причем корреляции Rii удовлетворяют системе n линейных уравнений

n A2

/■':: Y.—■ 2f:X: 1 < * < (16)

j=1 Pij

а взаимные корреляции Я^. = Я.г, г < ], и корреляции Квг при 1 < < п связаны системой г = (3п2 — п)/2 линейных алгебраических уравнений

г-1 п

Кгк = £ Ац К.к + £ Ац Кк. + /г Хк, к < г,

— 3=пг (17)

Каг = £ АцК.г + £ Ац. + /Х 1 < в, 4 < п,

3=1 3 = *

которая имеет единственное решение.

Доказательство. Так как имеем £(0) = /, вторые моменты будут рассмотрены только при то > 1. Пользуясь равенствами (6), находим второй момент разности

/ Ат) с \ тл I ^(т)\ Лг(т)

(С /г) и с учетом Е I Сг ) = Х> получаем

2 г — 1 Л2 2 п л2 2

Е (с(т)) =Е—Е(Сзт)) +Е—Е(^1)) то > 1. (18)

3=1 ц=г

Равенство (18) дает

г— 1 п

КГ] — КТ—1) = £ Вгз (Ц) — К.— 1)) + £ Вгз — 1) — Я.^ +

3=1 3=г

г(т) ~,г(т—1)л

+ 2fi (X(m) - X(m-1)) . (19)

Покажем, что для всех г будет выполняться \Я(т) — Я™ ^^ 0 при то ^ ж.

В [6] доказано, что в методе Зейделя для разностей итераций X(т), X(т—1) справедливы неравенства

||Х(т) _х(т-1)|| <мт-1||Х(1) _Х(0)||; ^ = тах ^ -Ы-> М<||А||<1.

г 3=1 1 — £ А \

3=г

Используем аналогичный способ. Пусть имеем

г— 1 п

^ О ^ О тахг Vг , • п ч

3=1 3=г г

Можно показать, что справедливо неравенство V < \\В\\ < 1. Из (19) следует, что для любого г выполняется

Е(С(ТО)}2 - Е(с(т—1))2

Е(т) - Е(т-1)

<

<

%(т) — %(т—1) < у X(т—1) — X(т—2)

\х(т) - X(т—1) I

Итерируя неравенство (20), получаем

X(т) — X(т—1)

<

X(1) - X(0)

т — 1

= Е ^ »т—1—.

К

(21)

Так как имеем 0 < и,^ < 1, несложно показать, что существуют константы 0 < в < 1 и С > 0 такие, что выполняется ат < С вт—1. Тогда получим

X(т) — X(т—1)

< С1 вт—1, С1 =

X(1)- X(0)

211/11

К '

что гарантирует существование конечного предела

Иш Е(т) = Иш Е (с(т)6т)) = Ягг-

(22)

Переходя при т ^то к пределу в равенстве (18) и учитывая предельные равенства (11), убеждаемся, что корреляции Нгг удовлетворяют системе линейных уравнений (16).

Матрица Д симметричная и ее диагональные элементы можно считать найденными, поэтому для корреляционной матрицы Д достаточно найти элементы Н^ при к < г. Чтобы вычислить корреляции векторов С(т) и С(т—1) и взаимные корреляции компонент, обратимся к формулам (6) или (7) и найдем при 1 < г, к < п следующие математические ожидания:

1—1 п

Е ((Ст) - /гК^) = £ А^Е (с(т)ст) +53 АъЕ (с(т—1)с(т)), к < г,

3=1 з=г

Е ((С(т) - /*К1т—1)) = Е АъЕ (Ът)С(т—1)) +

3=1

п ( )

+ £ А(3Е (с3т—1)С(т—, 1 < М < п. (23) 3=(

Перейдем к доказательству равенств (15). Введем обозначения:

г—1

К = ш1п | 1 - ^ 1Аз I I , Н(т) = Е (Сг(т)С(тП , к(т) = Е и(т)С(т—1)

3=1

Ф(т)

шах

3

н(т) - н(т—1)

Ф(т)

шах

3

К(т) - к(т—1)

а

т

т

Пользуясь равенствами (23) и симметрией ковариаций Д™1^, получаем неравенства, аналогичные (20),

|Д(») _ д(--1>| < ф(т) < ф(ш) + WfW ^Н .

1 ij ij 1 --- Г ^ ' / -Л

\K(f _ < ф(т) < мф(т-1) + II/II •Il^(m-_1)-^(m-2)ll) 1 < г, j < П.

(24)

Итерирование (24) приводит к неравенствам, аналогичным (21), что с учетом доказанных пределов (22) ведет к доказательству предельных равенств (15).

Если в равенствах (23) перейти к пределам при m ^ ж, получим равенства (17).

Убедимся, что система уравнений (17) имеет единственное решение. В каждой строке системы (17) ненулевыми элементами являются (n — 1) элемент одной строки матрицы A, так как все Rii, 1 < i < n, не входят в качестве неизвестных в эту систему. Отсюда с учетом неравенства (2) следует, что норма матрицы системы (17) меньше единицы.

Замечание. Интересно отметить, что элементы корреляционной матрицы R = \\Hij || предельного вектора lim Z(m) тесно связаны с предельными корреляциями K = \\Kij\\ двух последовательных приближений Z(m) и Z(m-1). Те и другие корреляции входят в одну систему линейных алгебраических уравнений и могут быть вычислены только одновременно. Матрица R является симметричной, матрица K — нет. Отсутствие симметрии предельной матрицы K является следствием того, что

х , , „ „ - х х - (Zi Zj ) и E(Zj 'Q )

вычисляются по разным формулам (см. (23)). В частности, будем иметь

Kg) = E (сРZT) = X^fj = j = E = X(1>U

Пример 1. Рассмотрим систему линейных алгебраических уравнений вида (1), где матрица А и вектор f определяются выражениями

0.3 —0.5 0.1 A = I —0.2 0.3 0.4 0.4 -0.3 0.2

0.1 —0.5 0.4

Норма матрицы А равна ||А|| = ш&х^ | = 0.9. Элементы матрицы Р вычис-

лялись по формулам

Ра =

^к = l \Aik \

1 < i,j < n.

(25)

Решением системы (1) является X = (0.5226, —0.3529, 0.8937)т. Предельный вектор lim Z(m) имеет среднеквадратические отклонения а = (0.8553, 0.9988, 0.8298)т и коэффициенты корреляций между компонентами pi2 = —0.2041,pi3 = 0.4902,p23 = —0.4219. Приведем корреляционные матрицы R и K с элементами (17):

1.0046 —0.3588 0.8858 R = I —0.3588 1.1222 —0.6651 0.8858 0.6651 1.4873

0.6216 —0.7705 0.8364 K = I —0.1390 0.4012 —0.2187 0.6766 0.7028 1.0551

f

При N = 106, М = 90 максимальные отклонения выборочных ковариаций См = Соу{См, ^м} от предельных ковариаций С. = См = Я. — ХгХ.

удовлетворяют неравенству шах. \См — Сц \ < 0.0027, а для соответствующих сред-неквадратических отклонений имеет место неравенство шахг \агг(М) — аг \< 0.0006.

Обозначим через р(М, N максимальное расхождение между компонентами точного решения X и найденного приближенного: р(М, N) = шахг

Х Лм)

Хг — ъ N

(г)

. В таб-

лице приведены значения этой функции для разных М и N. Жирным шрифтом выделены значения М, которые получены из равенства (14) при Н = 1.

N 10000 1000000

М 50 58 70 60 80 100

р(м, ло 0.0155 0.0071 0.0065 0.0013 0.0012 0.0011

Результаты, приведенные в таблице, показывают, что для достижения нужной точности при небольших N N = 104) необходимо выбирать М > М*, а при больших N N = 106 и больше), наоборот, брать М < М*, где М* — граница, полученная из равенства (14).

Пример 2. Пусть размерность системы равна п = 100. Пусть компонентами вектора / и вспомогательной матрицы размерности 100 х 100 являются независимые (псевдослучайные) равномерно распределенные на (0,1) числа. Элементы вспомогательной матрицы преобразовываем и отправляем в матрицу А, поменяв знаки части элементов на противоположные. Норма матрицы А равна \\А\\ = 0.9. Матрицу Р строим по схеме (25).

При моделировании процесса См, М = 80, оценки Х(м), 1 < г < п решений Хг были вычислены N = 106 раз. Кроме того, были найдены теоретические а2 и выборочные а2 дисперсии компонент вектора С. Для данного примера формула (14) дает М = 80.

Приведем результаты вычислений: шахг \Хг\ = 3.3839, шахг \Х(м) — Хг \ = 0.0066, шахг \а2 — а\ \ = 0.0212, шахг \аг — аг \ = 0.0045, шахг аг = 2.495.

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

Литература

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

2. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. М.: Академия, 2006. 367 с.

3. Товстик Т. М. К решению систем линейных алгебраических уравнений методом Гиббса // Вестн. С.-Петерб. ун-та. Сер. 1. Математика. Механика. Астрономия. 2011. Вып. 4. С. 90—98.

4. Беляева А. А., Ермаков С.М. О методе Монте-Карло с запоминанием промежуточных результатов // Вестн. С.-Петерб. ун-та. Сер. 1. Математика. Механика. Астрономия. 1996. Вып.3. С. 8-11.

5. Дмитриев А. В., Ермаков С.М. Метод Монте-Карло и асинхронные итерации // Вестн. С.-Петерб. ун-та. Сер. 1. Математика. Механика. Астрономия. 2014. Т. 1(59), вып. 4. С. 517-528.

6. Демидович Б. П., Марон И. А. Основы вычислительной математики. М.: Наука, 1981. 664 с.

7. Товстик Т. М., Волосенко К. С. Алгоритм Монте-Карло для решения систем алгебраических уравнений методом Зейделя // Труды межд. конф. «Актуальные проблемы вычислительной и прикладной математики-2015». Новосибирск, 2015. С. 771-776 (2015). ISBN 978-5-9905347-2-8.

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

Товстик Татьяна Михайловна — кандидат физико-математических наук, доцент; peter. tovst ik@mail.ru

Волосенко Ксения Сергеевна — студент; ksuhasvb@mail.ru

MONTE CARLO ALGORITHM FOR SOLUTION OF A SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS BY THE SEIDEL'S METHOD

Tatiana M. Tovstik, Kseniya S. Volosenko

St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation; peter.tovstik@mail.ru, ksuhasvb@mail.ru

For solution of a system of linear algebraic equations by the Monte Carlo method a method of the successive approximations is used. The next in turn iteration is modeled as a random vector, the expectation of which coincides with the corresponding Seidel's iteration. The existence and the boundedness of dispersions of the limiting vector is proved. The system of linear equations, to which satisfy the cross-correlations of the limiting vector, and the limiting correlations of two successive approximations is delivered. All these correlations are unknowns of the common linear system and may be found only simultaneously. The boundedness of dispersions of the limiting vector leads to the stochastic stability of algorithm. The number of iterations giving the desirable exactness is estimated. As an examples, the systems of the orders of n = 3 and n = 100 are studied. Refs 7. Tables 1.

Keywords: method Monte Carlo, Seidel's algorithm, system of linear algebraic equations.

References

1. Ermakov S.M., The Monte Carlo Method and Related Questions (Nauka, Moscow, 1975) [in Russian].

2. Mikhailov G. A., Voytishek A. V., Numerical Statistical Simulation. Monte Carlo methods. (Publ. Center Academy, Moscow, 2006) [in Russian].

3. Tovstik T.M., "On the solution of systems of linear algebraic equations by Gibbs's method", Vestnik St. Peterburg Univ.: Math. 44(4), 317-323 (2011).

4. Belyaeva A. A., Ermakov S. M., "On the Monte Carlo Method with remembering of the intermediate results", Vestnik St. Petersburg University. Ser.1. Mathematics. Mechanics. Astronomy Issue 3, 8-11 (1996) [in Russian].

5. Dmitriev A. V., Ermakov S.M., "Method Monte Carlo and non-synchronic iterations", Vestnik St. Petersburg University. Ser. 1. Mathematics. Mechanics. Astronomy 1(59), issue 4, 517-528 (2014) [in Russian].

6. Demidovich B.P., Maron I. A., Computational mathematics (Nauka, Moscow, 1981) [in Russian].

7. Tovstik T. M., Volosenko K. S., "Monte Carlo algorithm for solution of systems of linear algebraic equations by the Seidel's method", Actual problems of calculate and applied mathematics, 771-776 (Novosibirsk, 2015, ISBN 978-5-9905347-2-8).

Для цитирования: Товстик Т. М., Волосенко К. С. Алгоритм Монте-Карло для решения систем линейных алгебраических уравнений методом Зейделя // Вестник Санкт-Петербургского университета. Серия 1. Математика. Механика. Астрономия. 2016. Т. 3(61). Вып. 3. С. 440-448. DOI: 10.21638/11701/spbu01.2016.312

For citation: Tovstik T. M., Volosenko K. S. Monte-Carlo algorithm for solution of a systems of linear algebraic equations by the Seidel's method. Vestnik of Saint Petersburg University. Series 1. Mathematics. Mechanics. Astronomy, 2016, vol. 3(61), issue 3, pp. 440-448. DOI: 10.21638/11701/spbu01.2016.312