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

К РЕШЕНИЮ СИСТЕМ

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

Т. М. Товстик

С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, Peter.Tovstik@mail.ru

В статье представлен алгоритм приближенного решения системы линейных алгебраических уравнений методом Монте-Карло в сочетании с идеями моделирования полей Гиббса и Метрополиса [1]. Оценивается решение в виде ряда Неймана. Находится сразу весь вектор решений. Размерность системы может быть большой. Приводятся формулы для вычисления ковариационной матрицы отдельной имитации. Метод решения перекликается с методом, предложенным в статье Ермакова С. М. и Рукавишниковой А. И. [2]. На примерах систем третьего и сотого порядков производится сравнение точности приближения предложенного метода с методом из [2] и с классическим методом Монте-Карло [3], заключающемся в последовательном оценивании компонент неизвестного вектора.

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

здесь А = \\Aij\\І^=1п —квадратная матрица, а X = (хі,...,хп)Т и / =

(/і, . ..,/п)Т — векторы размерности п.

Представим решение X уравнения (1) в виде ряда Неймана, используя в качестве начального приближения вектор /,

Известно несколько достаточных условий сходимости этого метода [4]. Приведем два из них. В первом случае предполагается выполнение неравенства для т-нормы матрицы А,

X = АХ + /,

(1)

(2)

п

(3)

а во втором — для /-нормы матрицы А,

п

\\А\\і = тах ^\Акі\ < 1

1<І<П

(4)

к=1

Рассмотрим приближение к решению X вида

м

к=0

которое имеет точность \\— X(м)\\ < \\А\\м+1\\/\\/(1 -|\А\\) (см. [4]).

© Т.М.Товстик, 2011

Величину X(м) будем оценивать стохастической суммой

^(М) _ ^(0) + ^(1) + ... + £(м), (6)

в которой £(0) _ /, а последующие £(т), при т _ 1,..., М, будут оценками векторов Ат/, получаемых с помощью метода Монте-Карло.

С матрицей А свяжем стохастическую матрицу

П

Р _ \\PijII, И ^ 0, 1 - *,3 — п, _ 1, 1 - * - п.

^=1

Матрицу Р будем считать допустимой по отношению к матрице А, если р^ > 0 для тех пар (1,3), для которых А^ _ 0 или Aji _ 0, если выполнены, соответственно, неравенства (3) или (4).

В зависимости от того, каким неравенствам, (3) или (4), удовлетворяют элементы

матрицы А, подбираем стохастические матрицы Р или Q _ \ \qij \\, для которых будут

выполнены, соответственно, неравенства

\Аю\lpij < 1, *,3 _ 1 ••• ,п (7)

или

\Aji\/qij < 1, *,з _ 1 •••,п. (8)

В алгоритмах I и III, которые будет описаны ниже, должны быть выполнены

неравенства (7), а в алгоритме II — неравенства (8). Элементы соответствующих стохастических матриц Р и Q можно определить, например, следующим образом:

П П

Pik _ \/^2\А^\, qik _ \А^\/^2 . (9)

j=l j=l

Приведенные ниже примеры показывают, что неравенства (7) и (8) не являются необходимыми для применения соответствующих алгоритмов.

На матрицы Р и Q, если они допустимы для А, никаких дополнительных ограничений не накладывается. В частности, при соответствующей нумерации состояний они могут состоять из отдельных ненулевых блоков вдоль главной диагонали.

Рассмотрим два способа приближенного решения системы линейных уравнений методом Монте-Карло, на который накладываются идеи моделирования полей методом Гиббса и Метрополиса. Сокращенно будем их называть методами Гиббса и Гиббса—Метрополиса. Способы II и III приводятся для того, чтобы на примерах провести сравнения различных методов.

Случайное поле Гиббса определяется следующим образом. Имеется конечное число упорядоченных узлов 5, и в каждом узле в € 6 имеется конечное пространство состояний и8, так что и _ Пяеди8. Точки х, х € и называются конфигурациями. На и задается вероятностная мера П(х),х € и, выражающаяся через энергетическую функцию Н(х) в форме

П(х) _ ехр(—Н(х))/ £ ехр(—Н(и)).

пей

Если нужно найти конфигурацию, на которой достигается минимум функции Н(ж), то подбирают матрицу Р(ж, у) условных вероятностей перехода от конфигурации х к конфигурации у, предельное распределение у которой было бы П(ж). Тогда при разыгрывании любого начального распределения V(ж) и последующего многократного моделирования цепи Маркова с матрицей вероятностей перехода Р(ж, у) система приходит к конфигурации ио, для которой величина Н(ио) близка к минимальной.

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

Способ I. Метод Гиббса. Воспользуемся моделированием методом Гиббса, при котором на шаге т полностью обновляются все компоненты вектора С(т) в (6). Векторы ((т), т > 1, будут моделироваться последовательно. При каждом т для всех г (г = 1,...,п) согласно распределениям

Ра,Р12,---,Рт (10)

разыгрывается номер состояния; пусть это будет гт. Компоненты вектора С(т) находим по формулам

С(-)(г) = ^С(™-1)(гт). рит

В результате, ряд (6) принимает вид

М п Л-

с(м) = /+ЕЕ^с(т-1)Ые(г). (п)

1 - 1 РИт

т=1 г=1 т

Здесь е(г) —вектор, у которого на г-м месте стоит единица, а на остальных — нули. Вектор ((т) можно записать в виде

т

с(т) = ЬтС(т-1), С(т) = П Ьк I, т > 1, с(0) = I. (12)

к=1

Матрица Ьт в каждой строке имеет только один ненулевой элемент. В г-й строке

матрицы Ьт ненулевой элемент стоит на гт-м месте и равен Лцт /рцт. Равенство

(11) в обозначениях (12) записывается в виде

Мт

с(М) = I+Е Иьк I. (13)

т=1 к=1

Заметим, что условное математическое ожидание £(т)(г), если £(т-1) известно, равно

п

Е(С (т)(*)|С (т-1)) = Е Ап с (т-1)(з).

3=1

Так как Е(£(0)) = С(0) = I, выполняется Е(£(1)) = Л1 и Е(£(т)) = ЛтI, т = 2,...,М и, следовательно, Е(^(М))) = X(М).

При достаточно большом М можно оценить решение системы (1) в форме (2) усреднением N реализаций случайных векторов С,(М), а именно,

где (С(м) )8 —результат в-й имитации.

Найдем ковариационную матрицу Д(м) = ||Д(м)|| вектора £(м). Пусть

(^т/ )(і) — і-я компонента вектора Ьт/, ее дисперсия Д(т) = &(Ьт/)(і) равна

/ Ч 2

п л2 Л2 Л2

1)(т) — \ ' »то»то-1 »2»1 ^2^^

і = г РИт Pi-m.i-m.-l РІ2І1

1 < г < т. \ 1 < г < т.

” д.. д. . Л- •

Л**т Л*т*т-1 Л*2*1

^ = 1 Літ Pim.im.-1 РІ2І1

/(і1)

В наших обозначениях К(м) —дисперсия і-й компоненты вектора С,(м^ Найдем формулу для ее вычисления:

м м / м

Д(м} = Е]Т(Ьк/)(і)^(ЬШ/)(і) - (Ьк/)(і)} =

к=1 т=1 \ к=1 )

м м м-к / м ^

= Е]Т((Ьк/)(і))2 +2Е^ ]Г (Ьк/)(і)(Ьк+ /)(і) - е^2(ьк/)(і)

к=1 к=1 з = 1 \ к=1

м м м-к п

][>(ікЯ« + 2££ £ А„к+,-■ ■ А,„„+14к1,,

к=1 к=1 8 = 1 іг = 1

к + 1 < г < к + б

где при к =1

'І2 32 ~ ” У12^2^ І2

а при к > 2

^І13к + 1 = 5(ік+1 ,^к+1 )0ік:11 +

к —2 п I

'і(к — і—1)

1=0 іг = 1 4=0

+ ^ у ^ V РТ(1 ^(ік+1 — І,Зк+1 — І ))Дік+1-і Ік-ь ДЗк+1-ЬЗк_Ь З(ік — І ,Зк — і)Рік_і

1<г<к

Здесь 6(г^) — символ Кронекера.

Введем обозначение К(т,3') для ковариаций ЬтI(г) и ЬвI(?) элементов соответствующих векторов

3 = E((LmI)(г) - Е(Ьт!Х*))^)(,’) - E(LsI)(3)); (14)

(М)

тогда ковариации К\з можно записать в виде

м м м—т м м—т

*3) = Е Кгт) + Е Е к(т’т+з) + е Е к(т,т+8).

= 1 т=1 в = 1 т=1 в=1

2

2

При г = ] ковариации (14) находим по формулам

П

КТт) = Е А«-А«-дат^г1’+^т-1’)- *=^

Ът ?— т = 1

П

К—'к) = Е АЪЪ: • • • А*^+2^+1 (^(гк+ъ.Л' + ^1-)- * = ^ т > к.

к + 1 < г < т

Переходя к пределу при М — ж в равенствах (5) и (13), получаем для оценки X стохастический вектор £, имеющий ковариационную матрицу К, где

оо т

С = I + Е П Ьк I- ЕС = X, К = еоу(С, о = Иш Км’.

т=1 к=1

Способ II. Метод Ермакова С. М. и Рукавишниковой А. И. Данный алгоритм описан в работе [2] и, фактически, заключается в следующем. Вместо уравнения (2) рассматривается его транспонированный вариант

(Х)т = Ё ІТ (АТ )к.

(15)

к=0

Слагаемые суммы (15) оцениваются при значениях к > 1, а нулевое слагаемое равно I. Матрица Q = \\qik || вероятностей перехода в цепи Маркова определяется, например, по формулам (9), и для нее должны выполняться неравенства (8). Случайные состояния

*1 —— *2 —— • • • —— гм (16)

являются результатом моделирования соответствующей цепи Маркова со стохастической матрицей Q. Выбор начального распределения для моделирования :ч в (16), например, таков: п(г) = \/(г^/Е-^ \1 С?)\, 1 < г < п.

В статье [2] оцениваются члены ряда (15) при к > 2, а выбор начального состояния *2 в (16) происходит согласно распределению

М(і)

1/(*)1Е"=1 ІЛ

£м=і Азк І (*)!’ Для оценки X(м) имитируется ряд

1 < і < п.

П

(м)

І +

Ян) АІ2, тг(н) Япг

е(І2)

Аі

Яі2

е(*з) Н--------Ь

Аі

Чім-гі

-є(ім)

(17)

В правой части равенства (17) первое слагаемое I не случайное, математическое ожидание второго слагаемого равно

е(*2) = ЕЕ^-^>« = А/>

п(іі) ЯігІ2 І=1 Р1

а остальных — Ак І, 2 < к < М, поэтому Еп(м) = X(м). 94

Если п(М’ при одном и том же М промоделировать N раз, то

____ N

Х(м)«^М) = -Е(^м))8, (18)

в=1

где (п(м’ )я —результат в-й имитации. Для выборочного среднего используем, как и

в первом способе, обозначение пN’.

Если цепь Маркова со стохастической матрицей Q, являющейся допустимой по отношению к матрице А, имеет несколько эргодических классов, то должно выполняться условие N ^ п. Действительно, из формулы (17) следует, что при вычислении одной имитации г\(м’, если начальное состояние *1 попадает в какой-либо эргодиче-ский класс, система не выходит из него. В то же время, для получения хорошей оценки X(м’ нужно каждый эргодический класс посетить достаточно большое число раз.

Способ III. Классический метод Монте-Карло. Рассмотрим классический метод, не использующий поглощение. (Описание метода см. в книге Соболя И. М. [3].)

Для матрицы А выполнены неравенства (3), с ней связываем допустимую стохастическую матрицу Р. Каждая из компонент вектора X(м’ оценивается независимо от других. Для оценки, например, X(м’(2) полагается *о = 2, имитируется случайная последовательность (16), в которой переход *к — *к+1, к = 1,...п, происходит согласно распределению (10) с г = *к. После реализации ряда (16) получаем оценку

м А А А

£(м)сл = т + е —• —••• 1т-11т/(гт)

т=1 Рзч Рг1%2 Ргт-1Ът

компоненты X(м’(2) вектора X(м’.

Если сделать N независимых реализаций £(м’(2), то при больших М и N получаем приближенное равенство

____ N

вд«4м)(л = ^£(£(м)оо)«,

в=1

где (£(м’(2))8 —результат в-й имитации. В примерах элементы матицы Р вычисляются по формулам (9).

Способ IV. Метод Гиббса—Метрополиса. Для того, чтобы воспользоваться методом Метрополиса нужно ввести энергетическую функцию Н(У) вектора У. Определяем ее выражением

Н (У) = шах

1<<и

£А- У (2 ) + / (г) - У (г)

-=1

Заметим, что Н(X(м’) — норма погрешности приближения X(м’. Находим энергетическую функцию нулевого приближения Но = Н(X(0)) = Н(/).

На каждом шаге сначала моделируется пробный вариант. Будем предполагать, что он основан на способе I. На m-м шаге моделируется вектор Z(m), вычисляется вектор C(m), который считается пробной конфигурацией и должен пройти проверку. Если C(m) не проходит проверку, то моделируется новый вариант Z(m). Пусть V =

C(m-l), W = c(m).

Правила принятия пробной конфигурации:

а) если H(W) < H(V), то на m-м шаге принимаем пробную комбинацию C(m), а с ней и Z(m);

б) если H(W) > H(V), то с вероятностью exp( —(H(W) — H(V))) все равно принимаем пробную конфигурацию C(m);

в) если условия а) и б) не проходят, то пробная конфигурация Z(m) отвергается и моделируется новый вариант Z(m).

Если в качестве пробной конфигурации использовать вектор Z(m), полученный с помощью способа I, то вычисление H(C(m)) на каждом шаге увеличивает общее время счета, а точность, как увидим на примере 1, существенно не меняется.

Если пробную конфигурацию получать с помощью способа II, то вычисление H(W) = H(n(m)) можно сократить. Пусть

n

H(n(m)) = max \B(m)(i)\, B(m)(i) = V A3n(m)(j) + f (i) — f,(m)(i). i

3=1

Так как

'x(m) n(m-1) 1

\ ^ _ f(i0) m Aikik-1

Г-' = + h(im)e(im), h(im) = TT

n(i0) fc=i qik-iik

вычисления Н (г](т)) упрощаются ввиду рекуррентного соотношения В(т)(і) = Б(т-1)(і)+АіітН(іт) - 5(і, Іт)Н(Іт).

Приведем результаты вычислений по описанным алгоритмам. В первом примере п = 3, а во втором и третьем примерах п = 100.

Пример 1. Матрицу А возьмем симметричной, так что матрицы Р и Q с элементами (9) будут совпадать. Исходные данные системы (1) таковы:

A

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

0.І —O.B 0.З \ O.B

O.B —O.2 0.І I , f = 1 O.2

0.3 0.І 0.2 1 l 0.З

p(M, N) = max -X(i) — Z^^) (i)

l< i<r,

м N м N м N м N м N м N

40 104 40 10ь 60 10^ 60 10ь 100 104 100 10ь

I 0.00290 0.00064 0.01249 0.00420 0.00224 0.00115

II 0.01357 0.00074 0.04239 0.00453 0.00436 0.00090

III 0.00861 0.00026 0.03017 0.00167 0.00100 0.00089

IV 0.03792 0.04070 0.03429 0.04333 0.04060 0.04122

В таблице для четырех описанных выше способов и шести пар значений М и N приводятся величины р(М^). В данном примере Но = 0.26. В методе Гиббса— Метрополиса пробный вектор £(т) вычислялся способом I.

Пример 2. Размерность системы п = 100. В вектор / и во вспомогательную симметричную матрицу В размерности 100*100 засылаем независимые (псевдослучайные) равномерно распределенные на (0,1) числа. Элементы матрицы Р и Q = Р получаем по формулам (9), в которых АЪ- следует заменить на ВЪ-.

Для способов I и III элементы матрицы А при всех г и 2 сначала получаем в виде А— = 0.9р—, а затем примерно у 800 из них знак меняем на противоположный. У этих способов Но(I, III) = 0.52609. У способа II матрица в уравнении (1) будет транспонированной по отношению к только что найденной матрице А, и в данном случае Н0(П) = 0.57133.

Пример 3. Размерность системы п = 100. Матрицы Р и Q находим так же, как и в примере 2. Для способов I и III элементы матрицы А с точностью до знака определяем следующим способом: А— = 1.3р— при 1 < г < 10, 1 < 2 < по, по = {штг£к=1 Рг,к > 0.5}, при всех остальных значениях г и 2 полагаем А— = 0.9р—. Примерно 800-м элементам А— присваиваем отрицательный знак. В данном примере не выполняются условия (3) и (4). Нормы матрицы А таковы: \\А\\т = 1.1048, \\А\\; = 1.08744. Из всех 104 отношений А— /р— примерно в 500 случаях А—/р— = 1.3, т. е. не выполнено условие (7). Пример 3 показывает, что скорость сходимости к решению по сравнению с примером 2 уменьшается, однако методом Монте-Карло можно воспользоваться и в этом случае.

В таблице 2 для описанных выше трех способов при разных комбинациях М и N приводится величина расхождения р(М, N) между точным и приближенным решениями. Первые четыре столбца относятся к примеру 2, остальные — к примеру 3.

Таблица 2. Значения р(М, М) при п = 100

м N м N м N м N м N м N

60 104 60 10ь 60 10ь 60 104 60 10ь 60 10ь

I 0.02523 0.00788 0.00312 0.03918 0.01161 0.00362

II 0.42698 0.09447 0.03515 0.45595 0.12913 0.03693

III 0.07202 0.02996 0.01147 0.08299 0.03696 0.01290

При вычислении ^^М) способом II число операций примерно в п раз меньше, чем при других, хотя сравнение времени счета этот факт не подтверждает. Так, если і і — время счета способом I при М = 60, N = 105, а І2 —время счета способом II при М = 60, N = 107, то отношение их примерно таково іі : І2 = 4.5 : 17. При подсчете способом II получены следующие результаты: в примере 2 расхождение р(60,107) =

0.00956, а в примере 3 оно равно р(60,108) = 0.00506.

В примерах с размерностью системы n = 100, если у матрицы A все элементы больше нуля, то точность решения p(M, N) у способов I и III примерно одинаковая. Если имеются отрицательные элементы, то, как показывают примеры 2 и З, точность решения p(M,N) метода с элементами Гиббса лучше, чем у других методов. Метод с элементами Метрополиса увеличивает время счета, но не улучшает сходимости по сравнению с методом Гиббса. Метод Метрополиса можно также использовать при моделировании способами II и III, проверяя на каждом шаге точность счета, но и в этом случае проявляются те же недостатки.

Литература

1. Винклер Г. Анализ изображений, случайные поля и динамические методы Монте-Карло. Изд. СО РАН. 2002. (Winkler G. Image, Analysis, Random Fields and Dynamic Monte Carlo Methods. Berlin, Heidelberg: Springer-Verlag, 1995.)

2. Ermakov S. M., Rukavishnikova A. I. Application of the Monte Carlo and Quasi Monte-Carlo methods to solving systems of linear equations // Proc. 6th St.Petersburg Workshop on Simulation. St. Petersburg, 2009. P. 929-934.

3. Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1973.

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

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