Быстродействующие алгоритмы диагонализации трехдиагональных симметричных матриц на основе элементарных плоских вращений Текст научной статьи по специальности «Математика»

УДК 519.614

Быстродействующие алгоритмы диагонализации трехдиагональных симметричных матриц на основе элементарных плоских вращений

В.И. Иордан

Алтайский государственный университет (Барнаул, Россия)

High-Performance Algorithms for Tridiagonal Symmetric Matrices Diagonalization Based on Elementary Planar Rotations

V.I. Jordan

Altai State University (Barnaul, Russia)

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

БОТ 10.14258Лгуа8и(2017)1-15

This paper discusses two algorithms for diagonal-ization of real symmetric tridiagonal matrices. The proposed algorithms preserve the invariant (unchanged) tridiagonal form and provide high performance comparable with the most high-performance QR-algorithms and similar techniques. Contrary to the transformations of reflections, the presented algorithms of matrices di-agonalization using elementary planar rotations have absolute stability of calculations (Wilkinson stability). In this paper, the absolutely stable convergence of the proposed algorithms is proved. Convergence criteria are defined inherently by "integral" properties in contrast with "differential" properties of neighborhood convergence of QR-algorithms (in practice, QR-algorithms are prone to instability and loss of accuracy for eigenvalues computation under certain conditions). A priori estimates of maximum values of diagonalization error for proposed algorithms are provided. Also, the paper discusses the results of numerical experiments for a wide class of symmetric matrices conducnted to elaborate matrices di-agonalization time dependencies for QR-techniques and proposed algorithms.

Key words: tridiagonal symmetric matrix, planar rotations, high-performance, algorithms for matrices diagonalization.

Введение. Как известно [1], наиболее эффективные методы Гивенса и Хаусхолдера приводят квадратную симметрическую матрицу размерности п к трехдиаго-нальной симметричной форме за конечное число шагов (объем вычислений пропорционален п3). Алгоритмы диагонализации, использующие ортогональные преоб-

разования и сохраняющие инвариантной (неизменной) трехдиагональную форму, характеризуются меньшим объемом порядка п2. К таким эффективным алгоритмам относятся алгоритмы рЯ-метода, различающиеся стратегией выбора так называемого «сдвига» (сдвиги Релея, Уилкинсона и др. [1-7]). Несмотря на очень высокое

быстродействие и теоретически обоснованную устойчивость вычислений в рЯ-методе, при использовании «отражений» на практике, хотя и редко, но все-таки происходят случаи потери точности или неустойчивости вычислений [1, 8]. В случае патологически близких собственных чисел матрицы значения «сдвигов» могут приближаться то к одному, то к другому собственному числу матрицы (возможна неустойчивость сходимости [8]). То есть в рЯ-методе сходимость определяется близостью значения сдвига к собственному числу (малой окрестностью), иначе говоря, сходимость зависит от дифференциальных свойств плотности распределения спектра собственных чисел матрицы.

В данной работе рассматриваются два алгоритма диагонализации трехдиагональных симметричных матриц на основе плоских вращений (обозначение первого - М83БТв90, второго - М83БТ01ас), сохраняю-

щие инвариантной трехдиагональную симметричную форму и обеспечивающие абсолютную сходимость и устойчивость, а также высокое быстродействие, сравнимое с быстродействием алгоритмов рЯ-метода.

1. Математическое обоснование алгоритмов М83БТ090 и М83БТО^е. В предлагаемых алгоритмах производится последовательность так называемых «больших» итераций (БИ), состоящих из (п-1) элементарных плоских вращений в плоскостях с центрами (/,/+1), где / = 1,2,....,(п-1). На (к+1)-ой «большой итерации (БИ)», где к = 0,1,2,......, в результате

1-го вращения в плоскости (1,2), осуществляемого двусторонним матричным преобразованием с помо щью матрицы и[к+1\ нулевые значения элементов аЦ и сг^' матрицы. 11к1 пересчитываются в ненулевые элементы («выступы» г на рисунке а, элементы 3 диагоналей матрицы 5-го порядка обозначены символом х).

(*1

а) б) в)

Структурное отображение преобразований матрицы А(к) на (к+1)-ой БИ

г)

На фрагментах рисунка показаны результаты последующих вращений Гивенса в плоскостях с центрами (/, /+1), где / = 2,3,4. В результате двустороннего преобразования вращениями Гивенса [1] в конце (к+1)-ой БИ выступы г «сгоняются» за пределы матрицы и восстанавливается трехдиагональный вид:

А(к+1) = Т1 ■ А(к) • Т , (1)

к+1 к+1' 4 '

где Тк+1 — матрица ортогонального преобразования (Тк'+1 — транспонированная матрица) определяется факторизацией матриц элементарных вращений П'(к+1):

Т, +, = и(к+1) ■ и(к+1) ■■■ и(к+1). (2)

к+1 12 п-1 4 '

В матрице и.Ш) отличны от нуля только лишь элементы [1, с. 127]: а) в позициях (/,/) и (/ +1, / +1) определено значение е(к+1) — косинуса угла вращения;

б) оставшиеся диагональные элементы равны 1;

в) в позициях (/,/ +1) и (/ +1,/) определены, со-

(к+1) / (к+1К

ответственно, значения синуса и (-«/ ). Определенное число БИ производит в верхней части матрицы отделение двух собственных чисел для алгоритма М83БТв90, а для алгоритма М83БТв1ас —

трех собственных чисел (происходит «исчерпывание матрицы сверху» [1, 8]). То есть в верхней части матрицы «удлиняется» диагональная структура собственных чисел, а в нижней части сокращается трех-диагональная структура (рис., г). В конце всех этапов исчерпывания матрица имеет диагональный вид.

На каждой (к + 1)-ой БИ для первого алгоритма М83БТв90 начальное вращение в плоскости (1,2) осуществляется на + 90 градусов [8]

с1(к+1' = 0, sl(k+1' = + 1, (3)

Во втором алгоритме М83БТв1ас начальным является вращение Якоби [8], т.е. тангенс двойного угла начального вращения определяется выражением

{8(2-срГ)) = 2-Ь^/(а^-а^), (4)

где через Ь^к) обозначены элементы двух симметрично совпадающих кодиагоналей матрицы А(к), а через а(к) — элементы диагонали матрицы А(к). Значение угла вращения в плоскости (1,2) согласно (4) необходимо выбирать наименьшим по модулю. Отметим, что вращение (4) обнуляет элемент Ь(к),

и поэтому следующее вращение Гивенса для «сгона выступа вниз» (рис., а) по своей алгоритмической сути совпадает с начальным вращением алгоритма М83БТв90, так как с2к+1) = 0 и 5(к+1) = <к+1)), т.е. 52*+1) = 1. Последующие враще12ия для алго ритма М83БТв1ас в плоскостях (', /+1), где ' = 3,4,..., п-1, структурно и по алгоритмической сути совпа-

дают с последовательностью вращений алгоритма М83БТв90 с той лишь разницей, что эти вращения в алгоритме М83БТв1ас совершаются, начиная с плоскости (2,3). Ниже для этих алгоритмов приведены общие выражения (согласно рис.) для вращений в плоскостях (', /+1), где '=2, 3,..., п-1:

С целью уменьшения числа мультипликативных операций (умножений, делений) запишем алгоритм

в более рациональной форме [8]:

(5)

В (к + 1)-ой БИ диагональные элементы а'к), кроме первого и последнего, пересчитываются два раза вращениями в плоскостях (' -1, ') и (/,/ +1), а симметричные элементы Ь'к) двух кодиагоналей, кроме первого и последнего, пересчитываются три раза вращениями в плоскостях (' -1,'), (',' +1) и (' +1,' + 2). Элемент а'к) после вращения в плоскости (/',/' -1) обозначен как ~(к+1), а его результирующее значение а'к+1) в (к + 1)-ой БИ определяется вращением в плоскости (',' +1). Элемент Ь.к) вращением в плоскости (' -1,') пересчитывается как Ь.к) • с,1—+1), после вращения в пло-

скости (',' +1) он обозначен как Ь11 +), а после вращения в плоскости (' +1,' + 2) определяется его результирующее значение Ь(к+1) в (к +1) — ой БИ. Значения

„(к+1)

и 5

(к+1)

в алгоритмах М83БТв90 и М83БТ01ас определяются, соответственно, формулами (3) и (4); остальные же вращения - группой формул (5).

1.1. Сходимость алгоритма MS3DTG90. В каждой (к+1)-ой БИ вращение в плоскости (1,2) определяет: с(к+1) = 0, 5(к+1) = +1. Смена знака 5(к+1) результат (к+1)-ой БИ не изменяет (положим 51(к+1) = 1). Вращения в плоскостях (1,2) и (2,3) дают результаты [8]

= а\

х = {а\к) -ф)■ «>)г, Вш = [(а<*>-а^)2 ■ (с™)1 + (Ь^)2]"2.

Ь^НФ^У+^У)1'2.

(6)

(7)

(8)

Теорема 1. Алгоритм М83БТв90, определяемый начальным вращением (3) и соотношениями (5), абсолютно сходится.

Как известно [1, 8], евклидова норма матрицы и ее квадрат (сумма квадратов элементов матрицы) являются инвариантами матрицы (сохраняются в про-

Доказательство. При А —»со из (6) следует цессе унитарных преобразований). Следовательно,

.

элемент Ь1(™) ограничен некоторой положительной

константой В, то есть й1(™) < В. Поэтому знакоположительный ряд X (Ь2(к')2 является абсолютно сходящимся, а по необходимому условию сходимости ряда следует условие ¿2к) ^ 0 при к ^ да. Отсюда с учетом (6), (7), (8) следует, что ^) ^ 0 при к ^ да. Начиная с некоторых к, значением Ь2к) можно пренебречь, поэтому в верхней части матрицы отделяется блок 2-го порядка. Вращение Якоби в этом блоке определяет два собственных числа (исчерпывание матрицы сверху). Таким образом, для алгоритма М83БТв90 выше приведено доказательство абсолютной сходимости, гарантированной независимо от свойств «обусловленности» матрицы.

1.2. Сходимость алгоритма М83БТО^е. В каждой (к + 1)-ой БИ в плоскости (1,2) производится вращение Якоби [8]:

(9)

в результате чего элемент Ь1(к+1) = 0. Следовательно, вращение в плоскости (2,3) аналогично начальному

вращению алгоритма MS3DTG90, так как с!,к+1) = 0

,(*+1) -,

и а2 = sign(s( )). Выбор определенного значения угла Ф(к1) из возможных связан с обеспечением сходимости алгоритма М83БТв1ас. В (к +1) -ой БИ начальное вращение в плоскости (1,2) определяет элемент йГ» как

«¡<+,) = 4' ■ {ct]) f + а[к) ■ (5Г11)2 +2 ■ Й;

] 'Ч

- («<*>+ а[ку)! 2 + [(о®- о«*')/ 2] - cos<2<p'1<+l)) + ЬШ-8ш(2ф{*+1))

в силу тригонометрических тождеств

(cft+,))2 =[1+COS(2(pi*+1))]/2 и (5]С*+1))2 =[1-соб(2ф^'1,)]/2.

Обозначим через signl = sign(a{k] + а[к)),

sign2 = ) И = Y+[(«;*'

Определив cos^tp^1») = signl ■ (<> -af)/(2 ■ SQ), sintZipi'^") = signl • h\k) I SQ, юторые, очевидным образом, удовлетворяют (9), можно записать

= {а\К>+а[к))!2 + signl SQ, = (а'Г + а\К>)/2 -signl ■ SQ,

(t+i).

„ш

b(tt,) -

= 0,

"i

= ■s][\+cos(2(pf*1))] /2, = signl ■ sign2 ■ Л/[1-со5(2ф,|<+1>

)]/2

(10)

Используя формулы (5), имеем результаты вращений в плоскостях (2,3), (3,4):

b(k+1) = Ь('

»г2

(ii)

(12)

Теорема 2. Алгоритм М83БТв1ас, опреде- sign2 = sign(b(k)), абсолютно сходится. ляемый начальным вращением (9), соотношени- Доказательство. С учетом (10) получаем нераями (5) и условиями = венство

анализируя которое, можно записать

(a((k+1))2 > (a((k))2 + (b(k))2, если |a((k^ > Ia2k)|; (a((k+1))2 > (a2k))2 + (b(k))2, если Ia2k^ > \a(k)|.

либо

То есть в любом случае (a((k+())2 > (a((k))2 + (b((k))2 и, устремив k ^ да, получим

(aD2 >(a(0))2 + £ (b((k))2.

(13)

s

Величина (а^)2 ограничена некоторой константой В, так как евклидова норма матрицы сохраняется в унитарных преобразованиях. Следовательно, знакоположительный ряд 2 (Ь, ) абсолютно схо-

Т-Г ^ ^ к=0 1

дится. По необходимому условию сходимости ряда ^ 0 при к ^ да. Начиная с некоторых к, значением Ь^к) можно пренебречь, и в верхней части матрицы отделяется собственное число, равное значению элемента а[к). Однако одновременно с условием ¿,(к) ^ 0 при к ^ да в соответствии с (11) и (12) и величина st¡k) ^ 0. Следовательно,

С<к) ^ 1 и начиная с некоторого к = р, в «асимптотике» (12) можно вы-

1 (к+1) //! (кК 2 ,, (кК 2

разить в виде Ь2 ' (Ь2 ) + (Ь3 ) , которое в силу свойства рекуррентности позволяет записать выражение (Ь2да))2 = (¿2р))2 + 2 (Ь3(к))2.

к=р 3

В силу сохранения евклидовой нормы матрицы в унитарных преобразованиях знакоположительный

ряд 2 (Ь3(к))2 абсолютно сходится. По необходимо-

к=р

му условию сходимости ряда Ь3(к) ^ 0 при к ^ да. Начиная с некоторых к значением Ь((к ) можно пренебречь, и тем самым отделяется блок 2-го порядка с центрами в плоскости (2,3), вращение Якоби в которой определяет два собственных числа. В итоге производится исчерпывание матрицы сверху в количестве трех собственных чисел. Таким образом, для алгоритма М83БТ01ас гарантируется абсолютная сходимость независимо от свойств «обусловленности» матрицы.

2. Оценки погрешностей и результаты численных экспериментов. Скорость сходимости алгоритмов (процесс исчерпывания) зависит от спектральных свойств исходной матрицы. Для алгоритма М83БТв90 с учетом формул (6)-(8) можно записать

^ = Ь2к)/Ь(к+1) = (Вк / Ь™). |*2к)| = Ош • |*2к)| , где

°к+1 = ви /Ь1(к+1). Так как Ь(к+1' растет, а ^5 и Ь{к) стремятся к нулю при к ^ да, тогда, начиная с некоторых к, величина Ок+1 становится «сжимающим оператором» (Ок+1 < 1). Чем меньше к, начиная с которого Ок+1 < 1, и чем Ок ближе к нулю, тем быстрее происходит отделение клетки второго порядка сверху, а значит и быстрее происходит полная диагонализа-ция матрицы.

Для алгоритма М83БТ01ас с учетом (11) выражение (9) можно переписать в виде

^(2<рГ0) = С,+1 где вм =2-\ьГ)\/ФГ -<>)• Если на первых итерациях разность а[к) - а(к) близка

к нулю, тогда из (9) следует |ср|к+1) = к / 4 и '" | = 42 / 2,

а из (10) следует |а<*+1)-а?+,)\=2-5£>=2$к).

В выражении (1() показано монотонное возрастание а(к) при к ^ да. Следовательно, при к ^ да

для величины \Ок+2| = 2 • |Ь2к)|/|а(к+1) - а2к+ч| в асимптотическом режиме сходимости будет выполняться условие \5Ш\< 1, а в соответствии с «замечательными пределами», а также с вышеполученным

выражением ■ ср|А+п) = С(+| ■ ', можно записать:

= \вк+1 • ^)|/2 или < ^>|/2.

Из последнего неравенства следует, что не более 24 итераций требуется для достижения «одинарной» относительной точности порядка 10-7 (при этом отделяется три собственных числа). То есть не более восьми итераций достаточно на отделение одного собственного числа, что подтверждается численными экспериментами [8]. Анализ вычислительных погрешностей алгоритмов М83БТв90 и М83БТ01ас довольно громоздкий. В [8] подробно изложен вывод априорных оценок погрешностей для этих алгоритмов, имеющих вид:

а) для алгоритма М83БТв90:

\\(АА)м-АА\\*61-к-п™-гг\\А\\, (14)

б) для алгоритма М83БТ01ас:

¡(Л,)«45-*-И- (15)

где л ^ л; = 0,6 • 10"7 — в режиме одинарной точности, либо £\ = 0,22 • 10~15 — в режиме двойной точности вычислений; ЛА и (ЛА )м — диагональные матрицы, соответственно, точных и вычисленных алгоритмами собственных чисел трехдиагональных симметрических матриц. Оценки (14) и (15) показывают максимально возможные значения погрешностей при определении собственных чисел матриц, соответственно, алгоритмами М83БТв90 и М83БТ01ас.

Вычислительные эксперименты для широкого класса трехдиагональных симметричных матриц показали следующие зависимости временных затрат [8]:

а) для СЖ.-алгоритмов: I = а п2 : б) для алгоритма М830ТС90: / = |38 • и14; в) для алгоритма М830ТС1ас: / = 0,8 • |38 • п 1. Параметр Ре оказался на два порядка больше параметра ае, но для алгоритмов М83ЭТС90 и М83БТ01ас в асимптотике роста п (например, уже при п >50) становится «регулярным» исчерпывание матрицы снизу [8], что и объясняет понижение пока-

зателя степени до 1,4. Поэтому выравнивание затрат на диагонализацию матриц алгоритмами М83БТв90, М83БТв1ас и ОЯ-методом происходит при п > 700.

Заключение. В задачах диагонализации трехди-агональных симметричных матриц могут быть использованы алгоритмы М83БТв90 М83БТв1ас, со-

храняющие инвариантной трехдиагональную форму и обеспечивающие абсолютную сходимость и устойчивость вычислений, а также высокое быстродействие, сравнимое с быстродействием ОЯ-метода при п > 700.

Библиографический список

1. Уилкинсон Дж.Х., Райнш К. Справочник алгоритмов на языке Алгол. Линейная алгебра. — М., 1976.

2. Парлетт Б. Симметрическая проблема собственных значений. Численные методы / пер. с англ. Х.Д. Икрамова и Ю.А. Кузнецова. — М., 1983.

3. Watkins D.S. The Matrix Eigenvalue Problem: GR and Krylov Methods // D.S. Watkins. — SIAM. — 2007.

4. Prodi G. Eigenvalues of non-linear problems // G. Prodi (ed.). — Berlin Heidelberg, 2010.

5. Новиков М.А. Одновременная диагонализация трех вещественных симметричных матриц // Известия ВУЗов. Математика. — 2014. — № 12.

6. Кочура А.Е., Подкользина Л.В., Ивакин Я.А., Нид-зиев И.И. Сингулярные матричные пучки в обобщенной симметричной проблеме собственных значений // Труды СПИИРАН. — 2013. — Вып. 3(26).

7. Кузнецов Ю.И. Проблема собственных значений симметричной теплициевой матрицы // Сибирский журнал вычислительной математики. — 2009. — Т. 12, № 4.

8. Иордан В.И. Эффективные методы определения энергетического спектра матриц большой размерности в задачах экспериментальной физики : дисс. ... канд. физ.-мат. наук: 01.04.01. — Барнаул, 2003.