Улучшение одной из оценок скорости сходимости метода Зейделя Текст научной статьи по специальности «Математика»

2019

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Т. 6 (64). Вып. 2

МАТЕМАТИКА. МЕХАНИКА. АСТРОНОМИЯ

МАТЕМАТИКА

УДК 519.612 МБС 65Р10

Улучшение одной из оценок скорости сходимости метода Зейделя

А. Н. Борзых

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

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

Для цитирования: Борзых А. Н. Улучшение одной из оценок скорости сходимости метода Зейделя // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2019. Т. 6(64). Вып. 2. С. 185-195. https://doi.org/10.21638/11701/spbu01.2019.201

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

Ключевые слова: метод Зейделя, одношаговый циклический процесс, СЛАУ, система линейных алгебраических уравнений, итерационные методы решения, сходимость метода Зейделя, оценка скорости сходимости метода Зейделя.

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

х = Ах + ],

(1)

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

выборе произвольного xo и последующем вычислении = Axk + f. Известно (см., например, [1]), что метод простой итерации сходится, если спектральный радиус матрицы A меньше единицы: р(А) < 1.

Одним из вариантов метода простой итерации является метод Зейделя (одно-шаговый циклический процесс). Он отличен от метода простой итерации тем, что переход от вектора xk к вектору xk+1 происходит не сразу, а за n «простых» шагов. Каждый «простой» шаг заключается в вычислении i-го элемента вектора x^+i, где i меняется от 1 до n. Для определения г-го элемента вектора x^+i используются уже вычисленные элементы с 1 по i — 1 вектора x^+i и элементы с i по n вектора x^. Формула перехода такова:

i— 1 n

xfc+i = ^3 aij xk+iaij xk+f • j=i j=i

Известно (см., например, [1]), что метод Зейделя эквивалентен методу простой итерации для равносильной системы x = Bx + g, где В = (E — L) — i (D + R), g = (E — L) — if, E — единичная матрица, L, D, R — нижнетреугольная (Left), диагональная (Diag) и верхнетреугольная (Right) матрицы, образующие матрицу А: А = L+D+R.

Следовательно, метод Зейделя сходится, если спектральный радиус B меньше единицы: р(В) < 1. Величина р(В), кроме того, является коэффициентом скорости сходимости.

На практике вычисление р(В) является столь же сложной задачей, как и решение СЛАУ: требуется матричное обращение, матричное умножение и вычисление наибольшего по модулю собственного числа. В силу этого для определения факта сходимости (возможности применения метода Зейделя) и примерной оценки скорости сходимости (количества шагов для достижения заданной точности) пользуются какими-либо достаточными условиями сходимости. Для метода Зейделя удобно вычисляемым достаточным условием является условие

i— i n

ц < 1, где ц = max= Ъ , /?¿ = ^ |a¿¿|, 7» = У^ \aij\-i 1 Pi . . .

3=i 0=1

Недостаток этой оценки — слишком «грубый» характер величины Ставится вопрос о способах «улучшения» которые не требовали бы существенных вычислительных затрат и давали бы меньшее значение более близкое к р(В).

Один из способов был предложен в работе [2]. Идея заключалась в домножении (1) слева на некоторую матрицу D:

Dx = DA(D—iD)x + Df,

и получении эквивалентной системы x = Ax + f, где x = Dx, A = DAD—i, f = Df.

Матрица D выбиралась из класса матриц перестановок (в каждой строке и каждом столбце все нули и одна единица) и «подбиралась» таким образом, чтобы значение ^(A) было наименьшим. Несмотря на факториальное количество различных вариантов матриц D, был предложен и доказан алгоритм, определяющий оптимальную матрицу D за время O(n2).

В настоящей работе мы продолжим исследование «оптимизации» используем тот же прием, но матрицу D выбираем из класса диагональных матриц, т. е. ставим

задачу поиска диагональной матрицы D, обеспечивающей минимум величины p(A), где Ä = DAD-1.

Отметим, что преобразование Ä = DÄD-1 изменяет p, но никоим образом не влияет на p(B), т.е. на «действительную» скорость сходимости метода Зейделя:

Ä = DÄD-1 = D(L + D + R)D-1 = DLD-1 + DD D-1 + DRD-1, B = (E — DLD-1)-1(DDD-1 + DRD-1) = D(E — L)-1(D + R)D-1 = DBD-1,

B - в ^ P(B) = P(B).

В этом смысле целью работы является получение более качественной оценки, а не более быстрого вычислительного алгоритма.

2. Предлагаемый алгоритм. Матрицу Ä строим как предел итерационного процесса:

Afc+1 = DkÄkD-1, где Ao = A, Dk(i, а) = diag(1,..., а,..., 1).

i

Тогда D = ... Dk ... D2D1D0. Параметры i и а на каждом шаге выбираем следующим образом: i = argminpi(Ak), то есть i определяет минимальное значение pi

i

матрицы Afc; а — решение задачи на минимум p(Ak+1):

p(Ak+1) = max pj (Ak+1) -—min . j

В этом смысле описанный алгоритм эквивалентен методу поиска экстремума

f (¿1,..., dn) = p(DAD-1) f-1'"''d"> min, D = diag(d1,..., dn),

использующему покоординатный спуск.

3. Практическая реализация предлагаемого алгоритма. Здесь и далее считаем, что в начальной матрице A нет нулевых элементов. Если это не так, то построение эквивалентной СЛАУ x = Ax + f, не содержащей нулевых элементов в A, является отдельной задачей, которую в данной работе не рассматриваем. Отсутствие нулевых элементов в A гарантирует отсутствие нулевых элементов во всех матрицах Ak .

Утверждение 1. Если а > 1, то pi(Ak+1) > pi(Ak) и Vj = i pj(Ak+1) < Pj (Ak), то есть осуществляемое преобразование при а > 1 увеличивает значение pi (в матрице Ak оно было минимальным) и уменьшает все остальные pj.

Доказательство. Обозначим A/i = yí — |aii| = Y^n=i+1 laij|. В таких обозначениях имеем pi{Ak) = ^ "^ . Преобразование Ак+\ = DkAkD¿Г1 умножает эле-

1 Pi

менты строки i на а, делит элементы столбца i на а, оставляет без изменения aii и все остальные элементы. Тогда получаем

^Ak+Í) = (2)

1 — а/ji

Рассмотрим три случая:

1) если г = 1, то вг = 0, 7г > 0, |агг| + а7 > |агг| + 7г, 1 - авг = 1 -вг,

2) если 1 < г < п, то вг > 0, 7® > 0, |«гг| + а7г > |агг| + 7г, 1 - авг < 1 - вг,

3) если г = п, то вг > 0, 7г = 0, |а«г| + а7г = |ай| + 7г, 1 - авг < 1 - вг-Во всех трех случаях получаем

(л N Ы + «7» ^ \ад\+1ъ (лл

= > т^аГ =

Первая часть утверждения 1 доказана.

Поскольку преобразование Ак+1 = ^АдО-1 изменяет не только строку, но и столбец, то меняются все остальные ^. Покажем, что они уменьшаются. Возможны два случая.

1. Случай ] < г. В строке ] изменяется элемент а^г, лежащий правее диагонали и влияющий на величину параметра 7], который уменьшается:

Ъ(Ак+1) = (ъ(Ак) ~ М) + -\ajil < ъ(Ак),

а

а параметр /3^ остается без изменения. Числитель дроби ^ ^ ^ уменыпает-

1 - вз

ся, следовательно дробь уменьшается: ^ (Ад+1) < ^(Ад).

2. Случай ] > г. В строке ] изменяется элемент а^г, лежащий левее диагонали и влияющий на величину параметра в], который уменьшается:

&(Ак+1) = (^(Ак) - Ы) + < ^(Ак),

а

I ац | + 7]

а параметр ^ остается без изменения. Знаменатель дроби —---— увеличи-

1 - вз

вается, следовательно дробь уменьшается: ^(Ад+1) < ^(Ад).

Утверждение 1 доказано. □

Замечание 1. Для обеспечения строгих неравенств ^г(Ад+1) > ),

(Ад+1) < (Ад) было использовано оговоренное ранее отсутствие нулевых элементов в матрице Ао и во всех последующих матрицах Ад.

Замечание 2. Доказательство утверждения можно построить иначе. Возьмем от выражения (3) производную по а:

,, N _ 7»(1 ~ аРг) ~ (\ац\ + а7»)(~А)2 _ Ъ + Рг\ац\2

Только одна из величин 7г и вг может быть нулем. Получаем неравенство ^'г(а) > 0, и ^г(а) есть строго монотонно возрастающая функция. Тогда при а > 1 имеем ^г(Ад+1) = ^г(а) > М®(1) = ^г(Ад). Аналогично показывается, что У] = г ^(а) есть строго монотонно убывающие функции (требуется отдельное рассмотрение случаев ] < г и ] > г). Данные вычисления также необходимы для утверждения 8 (см. ниже).

Утверждение 2. Решение задачи ^(Ak+1 ) —% min эквивалентно решению

задачи а = max aj, где aj — решения уравнений ^j (Ak+1 ) = ^ (Ak+1). j

Доказательство. Рассмотрим pj как функции от а. По замечанию 2 имеем, что pi (а) есть строго монотонно возрастающая функция, pj (а) при j = i — строго монотонно убывающие функции. На рис. 1 представлен пример для матрицы размерности 6.

Очевидно, что необходимое условие минимума р(а) = maxpj (а) есть условие

j

pj (а) = р^а) для одного из j, а именно того, при котором значение pj (а) будет наибольшим, что соответствует наибольшему aj. Доказано. □

Замечание 3. Может возникнуть иллюзия, что достаточно решить лишь одно

уравнение pp^) = pi (а), где p = argmax pj (Ak), т.е. искать пересечение самого

j

верхнего убывающего графика pp(а) с возрастающим графиком p{(а). Это не так, из неравенства pp(1) > pj (1), j = i, не следует pp(а) > pj (а), т.е. графики могут «меняться местами».

Замечание 4. Тем не менее, существенную часть уравнений можно не рассматривать. Пусть найден корень ар уравнения pp^) = pi^), где p = argmax pj (Ak).

j

Тогда уравнения pj (а) = pi^), для которых pj (1) < pi(аp), можно не рассматривать, так как их корни заведомо меньше ар. Численный эксперимент показывает, что на матрицах размерности 100 средний выигрыш из-за устранения лишних уравнений составляет 92%, а на матрицах размерности 1000 — 98%.

4. Итоговая вычислительная схема предлагаемого алгоритма.

1. Положим к = 0. Для всех i вычислим параметры Дг, pi начальной матрицы Ао. Время вычисления — O(n2).

2. Положим i = argmin pj (Ak), s = argmax pj (Ak).

jj

3. Найдем корень а8 уравнения ps(Ak+i) = pi(Ak+i). Рассмотрим два случая.

А. Если s < i, то в строке s меняется элемент asi, лежащий правее диагонали. Следовательно, меняется 7S : 7s(Afc+i) = A/s(Ak) — |aSj| + ^|aSj|. Уравнение ps(Ak+i) = pi(Ak+i) имеет вид

Ks| + Ъ ~ lasi| + ¿|asi| \аы \ + агц

1 - вз 1 - авг '

что дает квадратное уравнение относительно а: С2а2 + с\а + со = 0, где С2 = 7г(1 - вз) + Ьвг, С\ = |ай| - вз1ац1 - Ь + вМвг^ со = -|а8г|, Ь = |а88| + у - К»|. Из двух вещественных корней берем наибольший.

Б. Если в > г, то в строке в меняется элемент азг, лежащий левее диагонали. Следовательно, меняется /38: ¡Зв{Ак+1) = ¡в8(Ак) — |а8г| + Уравнение Рз(Ак+1) = рг(Ак+1) принимает вид

|азз 1 + 7з _ 1агг1 + а1г

1 - /3S + |asi| - ¿|asi| l-o-Pi

Соответствующее квадратное уравнение таково: С2а2 + с\а + со =0, где С2 = t7i + fti(\ass | + 7s), ci = 11 a jj | — 7i|asi|-|ass|-7s, со = — Ki||aii|, t = 1 — вs + Ki|-Из двух вещественных корней также берем наибольший. Сложность данных вычислений — O(1).

4. Для всех j = i,s выполняем следующее:

а) если fij < jj,i{as) = "_а ^ , полагаем aj = 1 и ничего не вычисляем;

б) в противном случае аj вычисляем как корень pj(Ak+i) = pi(Ak+i). Схема вычисления описана в пункте 3. Сложность вычисления всех векторов аj в «худшем» случае — O(n).

5. Положим а = max аj, где аj — найденные корни уравнений pj (Ak+i) =

j

pi(Ak+i).

6. Выполняем преобразование Ak+i = Dk AkD-i, Dk(i, а) = diag(1,... ,а,..., 1). Преобразование изменяет одну строку и столбец, сложность вычислений — O(n).

7. Произведем перерасчет Pj, Yj:

а) для j = i: Pj (Ak+i) = авj (Ak), 7j (Ak+i) = а7] (Ak);

б) для j < i: /Зз(Ак+1) = Pj(Ak), 7j(Ak+1) = 7j(Ak) - +

в) для j > i: Pj(Ak+1) = fij(Ak) - \a>ji \ + 7j(Ak+1) = 7j(Ak).

+Ъ(Ак+1.

А после проведем перерасчет всех pf. ¡j,j(Ak+1) = ■

Сложность данных вычислений — О(п).

8. Обновим В и к: В = ВкВ, к = к +1.

9. Проверяем критерий остановки: количество выполненных шагов или иной. Переходим к пункту 2.

На рис. 2 показан ход работы алгоритма для матрицы размерности 10. По оси абсцисс отложен номер шага алгоритма, по оси ординат — значения ^ на соответствующих шагах алгоритма (10 графиков). Данный пример показывает уменьшение

"у

О 5 10 15 20 25 30

Номер шага алгоритма Рис. 2. Значения на различных шагах алгоритма.

¡(А&) с величины примерно 1.3 (сходимость метода Зейделя не ясна) до величины примерно 0.9 (метод Зейделя гарантированно сходится).

5. Доказательство сходимости. Доказательство сходимости сформулируем в виде лемм 1-2 и утверждений 3-9. Подробные доказательства каждого из утверждений оставим предметом для последующей публикации.

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

Лемма 2 (об увеличении наибольшего собственного числа матрицы с неотрицательными элементами). Если А > 0, В > 0, то Атах(А + В) > Атах(А).

Утверждение 3 (необходимое условие оптимальной матрицы А). Если значение ¡¿(А) минимально (т. е. матрица А оптимальная), то Уг, ] ¡г(А) = ¡3(А).

Утверждение 4 (существование и единственность оптимальной матрицы А). Матрица А существует и является единственной с точностью до знаков элементов.

Утверждение 5 (необходимое и достаточное условие оптимальной матрицы А). Матрица А оптимальная тогда и только тогда, когда Уг, ] ¡г(А) = ¡з (А).

- ! 1 -Исходное ц

--По — По еле 200 шг еле 300 ше гов гов

SS /

Áy

/ «"V

у/ ^

*

¡¡Г***6"*

Дисперсия х103

Рис. 3. Зависимость д от дисперсии элементов матрицы.

Утверждение 6 (о двухсторонней оценке оптимального ß). Пусть

ß(A) = ß*, A — оптимальная матрица. Тогда справедлива оценка: minßi(A) <

i

ß* < maxßi(A), где A — начальная или любая другая матрица, получаемая алгоритмом .

Утверждение 7 (о двухсторонней оценке элементов матрицы Ak).

Пусть Ak получена на шаге k предложенным алгоритмом: Ak = DAD-1, D = diag(di,..., dn), A — начальная матрица, aps — ее элементы, aj — элементы Ak. Справедлива оценка:

Vi, j m2 < |akj| < 1, где m = min |aps|.

j p,s

Утверждение 8 (о двусторонней оценке a). Рассмотрим один шаг алгоритма: Ak+1 = DkAkD-1, Dk(i, a) = diag(1,..., a,..., 1), i = argminßs(Ak),

i s

a = maxaj, aj — решения уравнений ßj(Ak+1) = ßi(Ak+1). Справедлива оценка j

для a:

2

m2 1

— maxßs{Ak) - minßs{Ak) < a < —. 3 V s s J m2

Утверждение 9 (о линейной скорости сходимости алгоритма).

\ß(Ak+i) — < С ■ Jц(Ак) — ц*\, С = 1 — ß* = ц{А), А — оптимальная матрица.

6. Численные эксперименты. Численные эксперименты проведем на случайных матрицах размерности 100, элементами которых являются случайные величины с нормальным законом распределения с математическим ожиданием 0 и дисперсией, которую варьируем от 0.001 до 0.01. Для каждой матрицы отмечаем ее исходное

0.4-

0.35

0.3

0.25

Рис. 4• Зависимость ^ от размерности матрицы.

20 40 60 80 100 120 140 160 180 200 Размерность п

значение а также значения которые получены после 100, 200 и 300 шагов предложенного алгоритма. Результаты представлены на рис. 3.

Следующую серию экспериментов проведем на матрицах размерностей от 10 до 200. Элементами матриц являются случайные величины с нормальным законом распределения с математическим ожиданием 0 и дисперсией Выполнение предложенного алгоритма осуществляем за n, 2n и 3n шагов. Результаты представлены на рис. 4.

Из проведенных экспериментов видим, что предложенный алгоритм дает уменьшение скорости сходимости ^ в среднем на 40%, требуя на вычисления порядка O(n2) операций.

Последнюю серию экспериментов проведем на матрице фиксированной размерности 100 с фиксированной дисперсией элементов 0.005, демонстрируя скорость сходимости предложенного алгоритма. Результаты представлены на рис. 5.

По оси абсцисс отложен номер шага алгоритма. По оси ординат, представленной в логарифмическом масштабе, отображаются:

1) фактическая ошибка: ) — , где значение получено этим же алгоритмом за 2000 шагов;

2) апостериорная теоретическая ошибка: max ^(Ak) — min ^(Ak);

i i

3) априорная теоретическая ошибка: — ^^ ^тах/хДД)) — тт/хДД))^ .

Практически горизонтальный вид априорной ошибки обусловлен тем, что величина т5 весьма мала, величина 1 — близка к единице. В этом смысле пред-

ставленное доказательство линейной скорости сходимости имеет преимущественно теоретическую ценность доказательства самой сходимости, а не практический способ оценки ее скорости. Скачкообразный вид апостериорной теоретической ошибки обусловлен тем, что в ходе работы алгоритма обеспечивается уменьшение max /ii(Ak)

Рис. 5. Фактическая, априорная и апостериорная ошибки.

на каждом шаге алгоритма, однако увеличение min p,¿ (Ak) на каждом шаге не га-

i

рантировано.

7. Открытые вопросы. 1. В данной работе переход от системы x = Ax + f к системе X = AX + f осуществлялся посредством преобразования Al = DAD-1, где D — диагональная матрица. В работе [2] использовался подход, в котором матрицей D являлась матрица перестановок. Каждый из подходов, как видно из результатов численных экспериментов, дает существенное уменьшение р. Интересен вопрос об объединении этих подходов, а также использовании каких-либо иных видов матрицы D.

2. Продемонстрированные подходы дают уменьшение р на 20-50%. Интересен вопрос о том, насколько данное уменьшение можно считать значимым, насколько получаемое значение р близко к истинной скорости сходимости, определяемой спектральным радиусом матрицы B = (E — L)-1(D + R). Численные эксперименты показывают, что в зависимости от вида матрицы A возможны как почти точные совпадения «улучшенного» значения р и величины p(B), так и их существенные различия, что оставляет вопрос о достоверной оценке скорости сходимости метода Зейделя актуальным.

3. Предложенный алгоритм является итерационным, и, следовательно, интересен вопрос об оценке его вычислительной сложности, например, в сравнении со сложностью метода Зейделя, для которого он получает оценку. Очевидно, что алгоритм, получающий оценку скорости сходимости, должен иметь существенно меньшую сложность, чем алгоритм, для которого эта оценка получается (иначе в нем смысла). Для матриц размерности 100, описанных выше, эксперименты показывают, что вычислительная сложность составляет примерно 10% от сложности метода Зейделя.

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

Литература

1. Фаддеев Д. К., Фаддеева В.Н. Вычислительные методы линейной алгебры. СПб.: Изд-во Лань, 2002.

2. Борзых А. Н. Улучшение одной из оценок скорости сходимости метода Зейделя путем выбора оптимального порядка уравнений системы линейных алгебраических уравнений // Ж. вычисл. матем. и матем. физ. 2017. Т. 57, №1. С. 3—8.

Статья поступила в редакцию 8 сентября 2018 г.;

после доработки 10 декабря 2018 г.; рекомендована в печать 20 декабря 2018 г.

Контактная информация:

Борзых Алексей Николаевич — канд. физ.-мат. наук, доц.; [email protected]

Enhansing an estimate of the rate of the Seidel method convergence

A. N. Borzykh

St. Petersburg State University, Universitetskaya nab., 7—9, St. Petersburg, 199034, Russian Federation

For citation: Borzykh A. N. Enhansing an estimate of the rate of the Seidel method convergence. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2019, vol. 6(64), issue 2, pp. 185-195. https://doi.org/10.21638/11701/spbu01.2019.201 (In Russian)

The article discusses the Seidel method for solving a system of linear algebraic equations and an estimate of the rate of the Seidel method convergence. It is proposed to construct an equivalent system for which the Seidel method also converges, but the rate of convergence is better. An equivalent system is constructed by a separate iterative process, where each single step requires O(n) operations. Stability of this iterative process is proved. Results of numerical experiments are presented showing an improvement of the estimate of the rate of convergence.

Keywords: Seidel method, one-step cyclic process, SLAE, system of linear algebraic equations, iterative solution methods, Seidel method convergence, estimate of the rate of the Seidel method convergence.

References

1. Faddeev D.K., Faddeeva V. N., Computational methods of linear algebra (Lan' Publ., St. Petersburg, 2002). (In Russian)

2. Borzykh A. N., "Improving an estimate of the convergence rate of the seidel method by selecting the optimal order of equations in the system of linear algebraic equations", Comput. Math. Math. Phys. (57) (1), 1-6 (2017). https://doi.org/10.1134/S0965542517010055

Received: September 8, 2018 Revised: December 10, 2018 Accepted: December 20, 2018

A u t h o r's i n f o r m a t i o n: Alexey N. Borzykh — [email protected]