close

Вход

Забыли?

вход по аккаунту

?

К решению систем линейных алгебраических уравнений методом Гиббса.

код для вставкиСкачать
УДК 519.2
Вестник СПбГУ. Сер. 1. 2011. Вып. 4
К РЕШЕНИЮ СИСТЕМ
ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
МЕТОДОМ ГИББСА
Т. М. Товстик
С.-Петербургский государственный университет,
канд. физ.-мат. наук, доцент, Peter.Tovstik@mail.ru
В статье представлен алгоритм приближенного решения системы линейных алгебраических уравнений методом Монте-Карло в сочетании с идеями моделирования
полей Гиббса и Метрополиса [1]. Оценивается решение в виде ряда Неймана. Находится сразу весь вектор решений. Размерность системы может быть большой. Приводятся формулы для вычисления ковариационной матрицы отдельной имитации. Метод
решения перекликается с методом, предложенным в статье Ермакова С. М. и Рукавишниковой А. И. [2]. На примерах систем третьего и сотого порядков производится
сравнение точности приближения предложенного метода с методом из [2] и c классическим методом Монте-Карло [3], заключающемся в последовательном оценивании
компонент неизвестного вектора.
Рассмотрим систему линейных алгебраических уравнений вида
X = AX + f,
(1)
здесь A = ||Aij ||i,j=1,n — квадратная матрица, а X = (x1 , . . . , xn )T и f =
(f1 , . . . , fn )T — векторы размерности n.
Представим решение X̄ уравнения (1) в виде ряда Неймана, используя в качестве
начального приближения вектор f,
X̄ =
∞
Ak f.
(2)
k=0
Известно несколько достаточных условий сходимости этого метода [4]. Приведем
два из них. В первом случае предполагается выполнение неравенства для m-нормы
матрицы A,
n
||A||m = max
|Aik | < 1,
(3)
1≤i≤n
k=1
а во втором — для l-нормы матрицы A,
||A||l = max
1≤i≤n
n
|Aki | < 1.
(4)
k=1
Рассмотрим приближение к решению X̄ вида
X (M) =
M
Ak f,
k=0
которое имеет точность ||X̄ − X (M) || ≤ ||A||M+1 ||f ||/(1 − ||A||) (см. [4]).
c
90
Т. М. Товстик, 2011
(5)
Величину X (M) будем оценивать стохастической суммой
ζ̂ (M) = ζ (0) + ζ (1) + · · · + ζ (M) ,
(6)
в которой ζ (0) = f , а последующие ζ (m) , при m = 1, . . . , M , будут оценками векторов
Am f , получаемых с помощью метода Монте-Карло.
С матрицей A свяжем стохастическую матрицу
P = ||pij ||,
pij ≥ 0,
1 ≤ i, j ≤ n,
n
pij = 1,
1 ≤ i ≤ n.
j=1
Матрицу P будем считать допустимой по отношению к матрице A, если pij > 0
для тех пар (i, j), для которых Aij = 0 или Aji = 0, если выполнены, соответственно,
неравенства (3) или (4).
В зависимости от того, каким неравенствам, (3) или (4), удовлетворяют элементы
матрицы A, подбираем стохастические матрицы P или Q = ||qij ||, для которых будут
выполнены, соответственно, неравенства
|Aij |/pij < 1,
i, j = 1, · · · , n,
(7)
|Aji |/qij < 1,
i, j = 1, · · · , n.
(8)
или
В алгоритмах I и III, которые будет описаны ниже, должны быть выполнены
неравенства (7), а в алгоритме II — неравенства (8). Элементы соответствующих стохастических матриц P и Q можно определить, например, следующим образом:
pik = |Aik |/
n
|Aij |,
qik = |Aki |/
j=1
n
|Aji |.
(9)
j=1
Приведенные ниже примеры показывают, что неравенства (7) и (8) не являются
необходимыми для применения соответствующих алгоритмов.
На матрицы P и Q, если они допустимы для A, никаких дополнительных ограничений не накладывается. В частности, при соответствующей нумерации состояний
они могут состоять из отдельных ненулевых блоков вдоль главной диагонали.
Рассмотрим два способа приближенного решения системы линейных уравнений
методом Монте-Карло, на который накладываются идеи моделирования полей методом Гиббса и Метрополиса. Сокращенно будем их называть методами Гиббса и
Гиббса—Метрополиса. Способы II и III приводятся для того, чтобы на примерах
провести сравнения различных методов.
Случайное поле Гиббса определяется следующим образом. Имеется конечное число упорядоченных узлов S, и в каждом узле s ∈ S имеется конечное пространство
состояний Us , так что U = Πs∈S Us . Точки x, x ∈ U называются конфигурациями. На
U задается вероятностная мера Π(x), x ∈ U , выражающаяся через энергетическую
функцию H(x) в форме
Π(x) = exp(−H(x))/
exp(−H(u)).
u∈U
91
Если нужно найти конфигурацию, на которой достигается минимум функции H(x),
то подбирают матрицу P (x, y) условных вероятностей перехода от конфигурации x
к конфигурации y, предельное распределение у которой было бы Π(x). Тогда при
разыгрывании любого начального распределения ν(x) и последующего многократного моделирования цепи Маркова с матрицей вероятностей перехода P (x, y) система
приходит к конфигурации u0 , для которой величина H(u0 ) близка к минимальной.
При моделировании случайного поля методом Гиббса матрица P (x, y) представляется в виде произведения матриц, каждая из которых обновляет конфигурацию в
одном узле. При этом обновление происходит во всех узлах в определенном порядке.
Такое моделирование по всем узлам называется сканом. Многократное сканирование
приводит поле к конфигурации u0 со свойствами, описанными выше.
Способ I. Метод Гиббса. Воспользуемся моделированием методом Гиббса, при
котором на шаге m полностью обновляются все компоненты вектора ζ (m) в (6). Векторы ζ (m) , m ≥ 1, будут моделироваться последовательно. При каждом m для всех i
(i = 1, . . . , n) согласно распределениям
pi1 , pi2 , . . . , pin
(10)
разыгрывается номер состояния; пусть это будет im . Компоненты вектора ζ (m) находим по формулам
Aiim (m−1)
ζ
(im ).
ζ (m) (i) =
piim
В результате, ряд (6) принимает вид
ζ̂ (M) = f +
M n
Aiim (m−1)
ζ
(im )e(i).
p
m=1 i=1 iim
(11)
Здесь e(i) — вектор, у которого на i-м месте стоит единица, а на остальных — нули.
Вектор ζ (m) можно записать в виде
ζ
(m)
= Lm ζ
(m−1)
,
ζ
(m)
=
m
+
Lk f,
m ≥ 1,
ζ (0) = f.
(12)
k=1
Матрица Lm в каждой строке имеет только один ненулевой элемент. В i-й строке
матрицы Lm ненулевой элемент стоит на im -м месте и равен Aiim /piim . Равенство
(11) в обозначениях (12) записывается в виде
ζ̂ (M) = f +
M +
m
Lk f.
(13)
m=1 k=1
Заметим, что условное математическое ожидание ζ (m) (i), если ζ (m−1) известно,
равно
E(ζ (m) (i)|ζ (m−1) ) =
n
Aij ζ (m−1) (j).
j=1
Так как E(ζ (0) ) = ζ (0) = f , выполняется E(ζ (1) ) = Af и E(ζ (m) ) = Am f, m =
2, . . . , M и, следовательно, E(ζ̂ (M) )) = X (M) .
92
При достаточно большом M можно оценить решение системы (1) в форме (2)
усреднением N реализаций случайных векторов ζ̂ (M) , а именно,
(M)
X̄ ≈ X (M) ≈ ζN
=
N
1 (M)
(ζ̂
)s ,
N s=1
где (ζ̂ (M) )s — результат s-й имитации.
Найдем ковариационную матрицу R(M)
(M)
||Rij || вектора ζ̂ (M) . Пусть
=
(m)
(Lm f )(i) — i-я компонента вектора Lm f , ее дисперсия Di = D(Lm f )(i) равна
⎛
⎞2
n
n
2
2
2
Aii Aim im−1
A
Aiim Aim im−1
Ai i
⎜ ⎟
(m)
m
Di =
· · · i2 i1 f 2 (i1 )−⎝
· · · 2 1 f (i1 )⎠ .
p
p
p
p
p
p
ii
i
i
i
i
ii
i
i
i
i
m
m
m−1
2
1
m
m
m−1
2
1
i = 1
i = 1
r
1 ≤ r ≤ m
r
1 ≤ r ≤ m
(M)
В наших обозначениях Rii
формулу для ее вычисления:
(M)
Rii
M
=E
k
(L f )(i)
=E
(L f )(i) −
k
M
M
2
k
(L f )(i)
=
k=1
2
((L f )(i)) + 2E
M M−k
k
(L f )(i)(L
k+s
f )(i) −
k=1 s=1
k=1
=
E
m
m=1
k=1
M
M
— дисперсия i-й компоненты вектора ζ̂ (M) . Найдем
D(Lk f )(i) + 2
k=1
E
M
2
k
(L f )(i)
=
k=1
M M−k
n
k=1 s=1
ir = 1
k+1 ≤ r ≤ k+s
(k)
Aiik+s · · · Aik+2 ik+1 Jik+1 i ,
где при k = 1
(1)
(1)
Ji2 j2 = δ(i2 , j2 )Di2 ,
а при k ≥ 2
(k)
(k)
Jik+1 jk+1 = δ(ik+1 , jk+1 )Dik+1 +
+
k−2
n
l
+
l=0
ir = 1
1 ≤ r ≤ k
t=0
(k−t−1)
(1 − δ(ik+1−t , jk+1−t ))Aik+1−t ik−t Ajk+1−t jk−t δ(ik−t , jk−t )Dik−t
.
Здесь δ(i, j) — символ Кронекера.
(m,s)
Введем обозначение Kij
для ковариаций Lm f (i) и Ls f (j) элементов соответствующих векторов
(m,s)
Kij
= E((Lm f )(i) − E(Lm f )(i))((Ls f )(j) − E(Ls f )(j));
(M)
тогда ковариации Rij
(M)
Rij
=
M
m=1
(14)
можно записать в виде
(m,m)
Kij
+
M M−m
m=1 s=1
(m,m+s)
Kij
+
M M−m
(m,m+s)
Kji
.
m=1 s=1
93
При i = j ковариации (14) находим по формулам
(m,m)
Kij
n
=
(m−1)
Aiim Ajjm (δ(im , jm )Dim
(m−1)
+ Jim jm ),
i = j,
im ,jm =1
(m,k)
Kij
=
n
(k)
Aiim · · · Aik+2 ik+1 (δ(ik+1 , j)Dj
(k)
+ Jik+1 j ),
i = j,
m > k.
ir = 1
k+1 ≤ r ≤ m
Переходя к пределу при M → ∞ в равенствах (5) и (13), получаем для оценки
X̄ стохастический вектор ζ, имеющий ковариационную матрицу R, где
ζ =f+
∞ +
m
Lk f,
Eζ = X̄,
R = cov(ζ, ζ) = lim R(M) .
M→∞
m=1 k=1
Способ II. Метод Ермакова С. М. и Рукавишниковой А. И. Данный алгоритм описан в работе [2] и, фактически, заключается в следующем. Вместо уравнения
(2) рассматривается его транспонированный вариант
(X̄)T =
∞
f T (AT )k .
(15)
k=0
Слагаемые суммы (15) оцениваются при значениях k ≥ 1, а нулевое слагаемое
равно f. Матрица Q = ||qik || вероятностей перехода в цепи Маркова определяется,
например, по формулам (9), и для нее должны выполняться неравенства (8). Случайные состояния
(16)
i1 → i2 → · · · → iM
являются результатом моделирования соответствующей цепи Маркова со стохастической матрицей Q. Выбор начального
распределения для моделирования i1 в (16),
n
например, таков: π(i) = |f (i)|/ j=1 |f (j)|, 1 ≤ i ≤ n.
В статье [2] оцениваются члены ряда (15) при k ≥ 2, а выбор начального состояния i2 в (16) происходит согласно распределению
n
|f (i)| j=1 |Aji |
,
1 ≤ i ≤ n.
μ(i) = n
k,j=1 |Ajk f (k)|
Для оценки X (M) имитируется ряд
AiM iM −1
f (i1 ) Ai2 i1
Ai3 i2
(M)
η
=f+
e(iM ) . . . .
e(i2 ) +
e(i3 ) + · · · +
π(i1 ) qi1 i2
qi2 i3
qiM −1 iM
(17)
В правой части равенства (17) первое слагаемое f не случайное, математическое
ожидание второго слагаемого равно
f (i1 ) Ai2 i1
e(i2 ) =
Aij f (j)e(i) = Af,
E
π(i1 ) qi1 i2
i=1 j=1
n
n
а остальных — Ak f, 2 ≤ k ≤ M, поэтому Eη (M) = X (M) .
94
Если η (M) при одном и том же M промоделировать N раз, то
X
(M)
≈
(M)
ηN
N
1 (M)
=
(η
)s ,
N s=1
(18)
где (η (M) )s — результат s-й имитации. Для выборочного среднего используем, как и
(M)
в первом способе, обозначение ηN .
Если цепь Маркова со стохастической матрицей Q, являющейся допустимой по
отношению к матрице A, имеет несколько эргодических классов, то должно выполняться условие N $ n. Действительно, из формулы (17) следует, что при вычислении
одной имитации η (M) , если начальное состояние i1 попадает в какой-либо эргодический класс, система не выходит из него. В то же время, для получения хорошей
оценки X (M) нужно каждый эргодический класс посетить достаточно большое число
раз.
Способ III. Классический метод Монте-Карло. Рассмотрим классический
метод, не использующий поглощение. (Описание метода см. в книге Соболя И. М. [3].)
Для матрицы A выполнены неравенства (3), с ней связываем допустимую стохастическую матрицу P . Каждая из компонент вектора X (M) оценивается независимо
от других. Для оценки, например, X (M) (j) полагается i0 = j, имитируется случайная последовательность (16), в которой переход ik → ik+1 , k = 1, . . . n, происходит
согласно распределению (10) с i = ik . После реализации ряда (16) получаем оценку
ξ (M) (j) = f (j) +
M
Ai
Aji1 Ai1 i2
i
·
· · · m−1 m f (im )
p
p
p
ji
i
i
i
i
1
1
2
m−1
m
m=1
компоненты X (M) (j) вектора X (M) .
Если сделать N независимых реализаций ξ (M) (j), то при больших M и N получаем приближенное равенство
(M)
X̄(j) ≈ ξN (j) =
N
1 (M)
(ξ
(j))s ,
N s=1
где (ξ (M) (j))s — результат s-й имитации. В примерах элементы матицы P вычисляются по формулам (9).
Способ IV. Метод Гиббса—Метрополиса. Для того, чтобы воспользоваться методом Метрополиса нужно ввести энергетическую функцию H(Y ) вектора Y.
Определяем ее выражением
n
H(Y ) = max Aij Y (j) + f (i) − Y (i) .
1≤i≤n j=1
Заметим, что H(X (M) ) — норма погрешности приближения X (M) . Находим энергетическую функцию нулевого приближения H0 = H(X (0) ) = H(f ).
95
На каждом шаге сначала моделируется пробный вариант. Будем предполагать,
что он основан на способе I. На m-м шаге моделируется вектор ζ (m) , вычисляется
вектор ζ̂ (m) , который считается пробной конфигурацией и должен пройти проверку.
Если ζ̂ (m) не проходит проверку, то моделируется новый вариант ζ (m) . Пусть V =
ζ̂ (m−1) , W = ζ̂ (m) .
Правила принятия пробной конфигурации:
а) если H(W ) ≤ H(V ), то на m-м шаге принимаем пробную комбинацию ζ̂ (m) , а
с ней и ζ (m) ;
б) если H(W ) > H(V ), то с вероятностью exp(−(H(W ) − H(V ))) все равно принимаем пробную конфигурацию ζ̂ (m) ;
в) если условия а) и б) не проходят, то пробная конфигурация ζ̂ (m) отвергается
и моделируется новый вариант ζ (m) .
Если в качестве пробной конфигурации использовать вектор ζ̂ (m) , полученный
с помощью способа I, то вычисление H(ζ̂ (m) ) на каждом шаге увеличивает общее
время счета, а точность, как увидим на примере 1, существенно не меняется.
Если пробную конфигурацию получать с помощью способа II, то вычисление
H(W ) = H(η̂ (m) ) можно сократить. Пусть
H(η̂ (m) ) = max |B (m) (i)|,
i
B (m) (i) =
n
Aij η̂ (m) (j) + f (i) − η̂ (m) (i).
j=1
Так как
η̂ (m) = η̂ (m−1) + h(im )e(im ),
h(im ) =
m
f (i0 ) + Aik ik−1
,
π(i0 )
qik−1 ik
k=1
вычисления H(η̂ (m) ) упрощаются ввиду рекуррентного соотношения
B (m) (i) = B (m−1) (i) + Aiim h(im ) − δ(i, im )h(im ).
Приведем результаты вычислений по описанным алгоритмам. В первом примере
n = 3, а во втором и третьем примерах n = 100.
Пример 1. Матрицу A возьмем симметричной, так что матрицы P и Q с элементами (9) будут совпадать. Исходные данные системы (1) таковы:
⎛
⎞
⎛
⎞
0.1 −0.5 0.3
0.5
A = ⎝ −0.5 −0.2 0.1 ⎠ ,
f = ⎝ 0.2 ⎠ .
0.3
0.1 0.2
0.3
Обозначим через ρ(M, N ) максимальное расхождение между компонентами точного
решения X̄ и найденного методом Монте-Карло приближенного решения. В обозначениях первого способа имеем
(M)
ρ(M, N ) = max X̄(i) − ζN (i) .
1≤i≤n
96
Таблица 1. Значения ρ(M, N ) при n = 3
I
II
III
IV
M
N
40
104
0.00290
0.01357
0.00861
0.03792
M
N
40
106
0.00064
0.00074
0.00026
0.04070
M
N
60
103
0.01249
0.04239
0.03017
0.03429
M
N
60
105
0.00420
0.00453
0.00167
0.04333
M
N
100
104
0.00224
0.00436
0.00100
0.04060
M
N
100
106
0.00115
0.00090
0.00089
0.04122
В таблице для четырех описанных выше способов и шести пар значений M и
N приводятся величины ρ(M, N ). В данном примере H0 = 0.26. В методе Гиббса—
Метрополиса пробный вектор ζ (m) вычислялся способом I.
Пример 2. Размерность системы n = 100. В вектор f и во вспомогательную симметричную матрицу B размерности 100*100 засылаем независимые (псевдослучайные) равномерно распределенные на (0,1) числа. Элементы матрицы P и
Q = P получаем по формулам (9), в которых Aij следует заменить на Bij .
Для способов I и III элементы матрицы A при всех i и j сначала получаем в
виде Aij = 0.9pij , а затем примерно у 800 из них знак меняем на противоположный.
У этих способов H0 (I, III) = 0.52609. У способа II матрица в уравнении (1) будет
транспонированной по отношению к только что найденной матрице A, и в данном
случае H0 (II) = 0.57133.
Пример 3. Размерность системы n = 100. Матрицы P и Q находим так же,
как и в примере 2. Для способов I и III элементы матрицы A с точностью до знака определяем
r следующим способом: Aij = 1.3pij при 1 ≤ i ≤ 10, 1 ≤ j ≤ n0 ,
n0 = {minr k=1 pi,k > 0.5}, при всех остальных значениях i и j полагаем Aij = 0.9pij .
Примерно 800-м элементам Aij присваиваем отрицательный знак. В данном примере не выполняются условия (3) и (4). Нормы матрицы A таковы: ||A||m = 1.1048,
||A||l = 1.08744. Из всех 104 отношений Aij /pij примерно в 500 случаях Aij /pij = 1.3,
т. е. не выполнено условие (7). Пример 3 показывает, что скорость сходимости к решению по сравнению с примером 2 уменьшается, однако методом Монте-Карло можно
воспользоваться и в этом случае.
В таблице 2 для описанных выше трех способов при разных комбинациях M
и N приводится величина расхождения ρ(M, N ) между точным и приближенным
решениями. Первые четыре столбца относятся к примеру 2, остальные — к примеру 3.
Таблица 2. Значения ρ(M, N ) при n = 100
I
II
III
M
N
60
104
0.02523
0.42698
0.07202
M
N
60
105
0.00788
0.09447
0.02996
M
N
60
106
0.00312
0.03515
0.01147
M
N
60
104
0.03918
0.45595
0.08299
M
N
60
105
0.01161
0.12913
0.03696
M
N
60
106
0.00362
0.03693
0.01290
(M)
При вычислении ζN способом II число операций примерно в n раз меньше, чем
при других, хотя сравнение времени счета этот факт не подтверждает. Так, если t1 —
время счета способом I при M = 60, N = 105 , а t2 — время счета способом II при
M = 60, N = 107 , то отношение их примерно таково t1 : t2 = 4.5 : 17. При подсчете
способом II получены следующие результаты: в примере 2 расхождение ρ(60, 107 ) =
0.00956, а в примере 3 оно равно ρ(60, 108) = 0.00506.
97
В примерах с размерностью системы n = 100, если у матрицы A все элементы
больше нуля, то точность решения ρ(M, N ) у способов I и III примерно одинаковая.
Если имеются отрицательные элементы, то, как показывают примеры 2 и 3, точность
решения ρ(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 MonteCarlo 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 г.
98
Документ
Категория
Без категории
Просмотров
6
Размер файла
294 Кб
Теги
решение, методов, уравнения, гиббса, система, линейный, алгебраический
1/--страниц
Пожаловаться на содержимое документа