Модифицированный метод Гивенса для ускоренного приведения действительной матрицы к форме Хессенберга Текст научной статьи по специальности «Математика»

УДК 512.64

Модифицированный метод Гивенса

для ускоренного приведения действительной

матрицы к форме Хессенберга

В.И. Иордан

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

The Modified Givens Method for Accelerated Reduction of a Real Matrix to the Hessenberg Form

V.I. Jordan

Altai State University (Barnaul, Russia)

Рассматривается проблема, связанная с повышением быстродействия широко известного метода Гивенса, основанного на элементарных плоских вращениях и преобразующего исходную несимметричную матрицу к форме Хессенберга (почти треугольной форме), а симметричную матрицу — к трех-диагональной симметричной форме. Предложенная модификация метода Гивенса за счет использования свойств рекуррентности пересчета некоторых элементов матрицы позволила уменьшить число операций умножений и тем самым увеличить быстродействие стандартного алгоритма Гивенса приблизительно в 1,4 раза. В стандартном и модифицированном алгоритмах Гивенса гарантируется устойчивость и точность преобразований исходной матрицы к форме Хессенберга, а накопленные ошибки округлений оказываются одного порядка. Численные эксперименты показали, что быстродействие модификации метода Гивенса практически сравнимо с быстродействием аналогичного по назначению высокоскоростного метода Хаусхолдера, который не всегда обеспечивает гарантированную точность и устойчивость преобразований матриц больших порядков к форме Хессенберга. Ключевые слова: метод Гивенса, симметричная матрица, ортогональные преобразования, плоские вращения, быстродействие, модификация алгоритма, форма Хессенберга.

The paper considers the problem of speeding up the widely known Givens method based on the elementary plane rotations. The Givens method reduces an initial asymmetric matrix to the Hessenberg (almost triangular) form, and an initial symmetric matrix — to the tridiagonal symmetric form. The proposed modification of the Givens method utilizes features of recurrent recalculation of certain elements of a matrix, and, thereby, the standard Givens algorithm speeds up to 1.4 times approximately. The standard and modified Givens algorithms provide guaranteed stability and accuracy of an initial matrix reduction to the Hessenberg form with cumulative round-off errors of the same order. Numerical experiments demonstrate the modified Givens method performance to be almost equal to the high-speed Householder method. However, the Householder method does not provide guaranteed stability and accuracy in all cases of large matrices reduction to the Hessenberg form.

Keywords: Givens method, symmetric matrix, orthogonal transformations, plane rotations, performance, modified algorithm, Hessenberg form.

DOI 10.14258/izvasu(2015)1.2-20

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

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

Модифицированный метод Гивенса...

ходной матрицей А имеет одинаковый спектр собственных значений [1].

В настоящей работе рассматривается модификация метода Гивенса [2, 3], ускоряющая «стандартный» алгоритм Гивенса [1] за счет использования свойств рекуррентности пересчета некоторых элементов матрицы, что позволяет уменьшить число операций умножений в 1,4 раза.

1. Математическое обоснование модификации метода Гивенса. Перед обоснованием модификации метода Гивенса в кратком виде представим основные этапы и особенности преобразований в стандартном алгоритме Гивенса [1, 2, 3].

Пусть Ат—1 обозначает матрицу, полученную из А0 = А после приведения столбцов с индексами 1,2,..., т — 1 к виду, изображенному на рисунке 1. В случае симметричной А0 кроме столбцов в симметричной матрице Ат—1 будут приведены и строки с теми же индексами (рис. 1б). На т-ом шаге преобразования матрицы Ат—1 к матрице Ат последовательно производятся вращения в плоскостях (т +1,т +1 + к), где к = 1,2,...,п — т — 1. В результате этих вращений аннулируются элементы в т-ом столбце под кодиагональю (для симметричной матрицы дополнительно аннулируются и соответствующие элементы в т-ой строке). То есть, как показано на рисунке 1, аннулироваться будут подчеркнутые в ниже изображенных структурах элементы (например, порядок матрицы п = 6 и номер шага т = 3).

Ат— =

А , =

X X | X X X X

X X | X X X X

0 X | X X X X

0 0| X X X X

0 0| X X X X

0 0| X X X X

а

X X | 0 0 0 0

X X | X — — 0

0 X | X X X X

0 0| X X X X

0 0| X X X X

0 0| X X X X

—

б

Рис. 1. Структуры матриц, полученные на (т - 1)-ом

шаге преобразований: а — для несимметричной исходной матрицы А0 = А ; б — для симметричной исходной матрицы А0 = А

Косинусы и синусы вращений в плоскостях (т +1, т +1 + к), соответственно, определяются следующим образом:

(0,т—1)

ск = Ь—/ ьк, $к = а(0+т—!' / Ьк,

к к—1 к ' к т+1 + к ,т к'

где

К = К2— +

(а(0,т—1)

^ т+1 + к,т

)Т

к =

(0,т—1)

(1)

(2)

В этих вращениях элементы т -го столбца а(т+1+кт, начиная с номера строки (т + 2) и до последней п -ой, обнуляются в количестве (п — т — 1). Последнее вращение в плоскости (т +1,п) определяет результирующее значение элемента а(0+™)т, равного значению Ъп—т—, = [(Ъ02 + £ (а%т—1))2]1/2. Если не-

}=т+2

(0,т—1)

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

Строка и столбец с индексами (т +1) пересчи-тываются при каждом вращении в плоскостях (т +1, т +1 + к), где к = 1,2,..., п — т — 1. Аналогично, строка и столбец с индексами (т +1 + к) пере-считываются при вращении в плоскостях (т +1, т +1 + к). Кроме того, определенные элементы строки и столбца с индексами (т +1 + к) пересчи-тываются предыдущими вращениями в плоскостях (т +1,т +1 + /), где } = 1,2,...,к — 1. Ниже приведен пересчет лишь (т +1) -ой и (т +1 + к) -ой строк в соответствии со стандартным алгоритмом Гивенса (столбцы пересчитываются аналогично, [2, 3]).

с, + а

а<) ,т—1) = а< 1—1,т—1)

т+1 + к,т+1 т + 1+к ,т + 1

() ,т—1) = а(0,т—1)

т+1 + к,т+1 +1 т+1 + к,т+1 + }

а(к,т—1) = (к—1т—1)

т+1,т+1 +1 т+1,т+1 + }

,(0,т)

(0,т—1) _ ^

т+1 + к,т+1 + ] ] '

Ск +

а— = а(к—1-т—1)

т+1 + к,т+1 +1 т+1 + к,т+1 + )

к = 2,3,., п — т — 1; ] = 1,2,., к — 1;

(1—1,т—1) Г — П • Ъ

1 т + 1 + к ,т+1 -3;'

(к—1,т—1) _

т+1 + к ,т+1 +1 ' к >

ск— а (кг,1' т-,1+. • $к,

к т+1,т+1 +1 к 3

(3)

где первый верхний индекс у пересчитываемых элементов связан либо с индексной переменной ] внутреннего цикла, либо с индексной переменной к внешнего цикла. Пока не завершен полностью шаг преобразования с номером т, второй индекс у пересчитываемых элементов остается равным ( т — 1). Элемент а^к т+1+. в формуле (3) в последующих итерациях цикла по к больше не пересчитывается, поэтому он представляет результат т-го шага преобразования, и его первый индекс обнуляется (для следующего (т+1)-го шага преобразования он считается как исходное данное).

Переходя к модификации алгоритма Гивенса, с учетом определения косинусов и синусов вращений

согласно (1), умножив на Ь. первую формулу в групповой формуле (3), получим

() ,т—1)

ат-"-; +1 • Ь. = а^+Г1' • ь. 1

т+1 + к ,т+1 . т + 1+к ,т + 1 1

, а(о,ш—1) • а(о,т—1)

т + 1+к, т + 1+. т + 1+., т '

(4)

В формуле (4) выполняется рекуррентное свойство для пересчитываемого элемента. Аналогично, умножив на Ьк третью формулу в групповой формуле (3), получим

ь = а(к—1-т—1) Ь

т+1, т+1 +) к ит + 1,т + 1+. ик—1

(к ,т—1)

+а(к—1,т—1) ^ а(0,т—1)

т+1 + к,т+1+. т + 1+к,т'

Введем обозначения

Я(.+п1+1к) +, = а (.+"1Т1к) +, • Ь.,

т+1 + к ,т+1 т + 1 + к ,т + 1 .

~хк = (а

(0, т—1)

т+1 + к, т+1 + к /~(к—1, т—1) т+1 + к ,т+1

— я

(к—1,т—1)) а<°.т—1) т+1,т+1 ' т + 1+к ,т

-Я

(к—1,т—1) \ _

а(0,т) = а(0,т—1)

т + 1+к ,т + 1+к т+1 + к ,т+1 + к

а(к,т—1) = а(к—1,т—1) + я т + 1,т + 1 т+1,т+1 чк

т+1,т+1 + к) Хк — вк • Як > Я, ,

• Ьк,

д(к,т—1) = я т + 1+к ,т + 1 к

Я(к,т—1) = Я т + 1,т + 1+к Лк

к = 1,2,.....,п — т — 1.

Ск Ят+1,т+1 + к / Ск '

С _ я(к—1,т—1) / С

Ск — Ят+1 + к,т+1 ' Ск

я(',т—1) _ я('—1,т—1)

"г ,т+1 "г ,т + 1

а(0'т) = а(0'т—+ •

г ,т+1+1 г ,т+1+1

(0,т—1) ^ а(0,т—1) г,т + 1+1 т+1+1 ,т ' •5(1—1,т—1)

г = 1,2,...,т; 1 = 1,2, ...,п — т — 1;

а(0,т) = й("—т—1,т—1) / Ь

г ,т+1 г ,т+1 п

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

г ^ т +1; а(0+7 = Ь

' т+1,т п

(10)

(11)

а^""' = ^(п—т—1,т—1)

т+1,т+1+1 т+1,т+1+1

^(0,т)

т+1,т (0,т)

/ Ь

,= 0; г = 1,2,., п — т — 1.

Я(.—1,т—1) = а(1,т—1) • ь т+1 + к ,т+1 т+1 + к ,т+1 .—1

~(к,т—1)

(к,т—1)

Ят+1,т+1 +. ат+1,т + 1+. ' Ьк :

^(к—1,т—1) = а(к—1,т—1) ^ Ь

т+1,т+1 +. т+1,т+1 +. к—1 '

д. = / Ь.—1, . = 1,2,...,п — т — 1.

(5)

(6) (7)

В начале т-го шага в ячейках памяти необходимо вместо величин а^—1 и а(°^—+¥1+1 сформировать величины

~(0,т—1) = т—1) "г ,т+1 "г ,т + 1

~(0,т—1) = а(0,т—1) т+1,т+1 + к Ит + 1,т + 1 + к '

Ь0 (г = 1,2,...п; г ^ т + 1), Ь0, (к = 1,2,....,п — т — 1).

Введенные обозначения позволяют записать модифицированный вариант алгоритма Гивенса

дО,т—1) = 1,т—1)

т + 1 + к ,т + 1 т+1 + к,т+1

_ а(0,т—1) ^ а(0,т—1)

т+1 + к ,т+1 + . т + 1+.,т'

а(.,т—1) = а(0,т—1)

т + 1 + к ,т + 1+. т+1 + к ,т+1 +. ^

~(1,т—1)

С- — я т+,'+к +, • в-,

т + 1+к ,т + 1 Л. '

Я(к ,т—1) = я(к—1,т—1) т + 1,т+1 + . т+1,т+1 + .

а(к—1,т—1) ^ а(0,т—1) т + 1+к ,т + 1+. т+1 + к,т'

(8)

а(0,т) = а(к—1,т—1)

т + 1 + к,т + 1+. т+1 + к,т+1 +.

(к—1,т—1)

Ск — я!,,,ml,)l. • вк,

к ^т+1,т+1 +.

к = 2,3, ...,п — т — 1; . = 1,2,., к — 1;

Я (;,т—1) = я( 1,т—1) т + 1,т+1 + к т+1,т+1 + к

_ а(0,т—1) ^ а(0,т—1)

т+1 + .,т+1 + к т + 1+.,т'

а(;,т—1) = а(0,т—1)

т + 1+. ,т + 1+к т+1 +. ,т + 1+к ^

Я (к ,т—1) = я(к—1,т—1)

т + 1+. ,т + 1 т+1 + . ,т + 1

а(0,т) = (к—1,т—1)

т + 1+.,т + 1+к т+1 + .,т + 1+к ' ^к

~(1,т—1)

с - — я 0+. +,+, • в.,

т + 1,т + 1+к 1}1 (к—1,т—1) ^ а(0,т—1)

т + 1+.,т+1 + к т+1 + к,т'

(9)

(к—1,т—1)

Ск — Я! , 1 • вк,

т+1 +. ,т + 1

к = 2,3,...,п — т — 1; . = 1,2,., к — 1;

В формулах (8) и (9) пересчитываются, соответственно, строки и столбцы с номерами (т +1) и (т +1 + к), за исключением «опорных» элементов из плоскости вращения (т +1, т +1 + к), которые пе-ресчитываются в формуле (10). Опорными элементами являются а(к+т—1\ , а(0'т)+, +,+, , а также элементы

т + 1т+1 ' т+1+к ,т+1+к

Я(k+"1_!) +. и Я(к+"—1+1+к, которые «восстанавливаются»

т + 1+к ,т + 1 т + 1,т + 1+к

в формуле (11).

Число операций умножений в модификации алгоритма Гивенса составляет

2 6 • (п — т) • (п — т —1)/ 2 +

т

3 • (п — т) •т=(5/ 2) • п3,

т

т. е. в 4/3 раза меньше, чем в стандартном алгоритме Гивенса [1], в котором число умножений — (10/3) • п3. Число же операций сложений и вычитаний не изменилось. Выражения для векторов матрицы преобразования Ят записываются следующим образом:

~(0,т—1)

(0, т—1)

Ьп,

як,т— 1) = я<к—1, т—1) + ^ а(0,т—1)

1,т+1 (0, )

,т+1 + к т + 1 + к,т 5 (к—1,т—1)

• вк >

(12)

(0, т—1)

^ = г • С — Я

1, т+1 + к 1, т+1 + к к г ,т+1

к = 1,2,..., п — т —1,

г(0,т) = я(п—т—1,т—1) / Ь

1, т+1 1, т+1 п—т—1'

г = 1,2,п.

Таким образом, алгоритм, определенный соотношениями (8)-(12), соответствует результирующему матричному преобразованию

Ат = ^т • А0 • ^т = Нт ' Ат—1 • Нт ,

Модифицированный метод Гивенса.

где Я'т — результат транспонирования матрицы Ят = Ят • Нт; Н — результат транспонирования матрицы Нт — матрицы преобразования (перехода) т -го шага.

Алгоритм (8)-(12) реализован программно на различных алгоритмических языках в режиме «двойной точности». Программы позволяют привести также и симметрическую матрицу А0 = А к трехдиагональ-ной форме. Однако используя симметрию таких матриц, алгоритм (8)-(12) можно упростить путем преобразования лишь верхнего треугольника элементов матрицы. Тогда число умножений в модифицированном алгоритме Гивенса за счет симметрии матрицы сокращается до ^ 6 • (п — т) • (п — т — 1) / 2 = п3,

т

а число операций сложений и вычитаний равно

^4 • (п — т) • (п — т — 1) / 2 = 2 • п3 / 3. Для симметрич-

т

ного варианта стандартный алгоритм Гивенса использует (4/ 3)- п3 умножений [1], таким образом, модификация алгоритма Гивенса и для симметричного варианта дает уменьшение числа умножений в 4/3 раза. Число же операций сложений и вычитаний совпадает с числом аналогичных операций стандартного алгоритма Гивенса. Модификация алгоритма Гивенса для симметричных матриц также программно реализована в режиме «двойной точности».

2. Анализ результатов численных экспериментов. Численная устойчивость стандартного и модифицированного алгоритмов Гивенса реализуется надежным вычислением величин ск, $к, qk, которые, в свою очередь, обеспечивают свойства ортогональности каждого элементарного вращения с гарантированной точностью. Сопоставление погрешностей вычислений, возникающих в стандартном алгоритме Гивенса и его модификации, главным образом определяется возникновением ошибок округлений

(}—1,т—1)

при вычислении пары произведений я„+1т+1+к', ] • ^ (для стандартного алгоритма) и пары

произведений е^+к • qj и ят^т+й; • qk (для ег°

модификации). В [3] проведен подробный анализ точности этих алгоритмов и показано, что погрешности соответствующих друг другу пар произведений одного порядка, а именно, они равны величине (5е1 + 7е2), где е1 — относительная погрешность, переносимая с предыдущей итерации на последующую; е2 — относительная погрешность при выполнении арифметической операции или квадратного корня [3].

В [1] с помощью достаточно громоздких выкладок учитывается монотонный переход от т-го шага преобразования к следующему, что в конечном итоге для метода Гивенса дает возможность оценить норму матрицы «невязки» М0 в преобразовании исходной матрицы А0 к форме Хессенберга. Оценка нормы «невязки» соответствует неравенству в виде

||м0|| < k1 • n312 • s • ||А0 , для метода Хаусхолдера норма «невязки» выше: ||М01| < k2 • n2 • s • ||A01|.

В численных экспериментах было использовано большое число матриц с различной структурой (ленточные и заполненные матрицы). Сравнение по точности и быстродействию стандартного и модифицированного алгоритмов Гивенса производилось не только между собой, но и с аналогичным по назначению высокоскоростным методом Хаусхолдера, который по быстродействию не менее чем в 1,5 раза превосходит стандартный метод Гивенса [1, 2, 3]. Однако в практических численных экспериментах оказалось немало матриц, для которых стандартные программы метода Хаусхолдера (например, из библиотеки программ на Фортране) давали результаты преобразований с большой потерей точности. Например, для некоторых матриц евклидова норма преобразованной матрицы в форме Хессенберга, которая должна быть инвариантной в процессе преобразований, содержала неверные цифры в первом знаке после запятой. Значительная потеря точности, как отмечают Пауэл и Рид [4], может возникать в случае, когда на m-ом шаге все элементы аннулируемого m-го столбца, с помощью которого строится матрица преобразования Хаусхолдера, оказываются превосходящими заданный порог малости eps и с большой относительной погрешностью.

Сравнение по точности и быстродействию модификации метода Гивенса со стандартным методом Гивенса и методом Хаусхолдера производилось с помощью соответствующих программ, реализованных на языках Фортран, PASCAL и C/C++. Для программной реализации стандартного метода Гивенса и метода Хаусхолдера были взяты за основу стандартные программы, соответственно, MSTG и MSTU из библиотеки языка Фортран, затем они были переведены на языки PASCAL и C/C++.

Численные эксперименты были проведены на IBM-совместимом компьютере с процессором Pentium 1,6 ГГц для ленточных и заполненных матриц со случайным распределением значений элементов.

Например, для ленточной структуры симметричных матриц (все элементы ленты равны 1, а вне ленты — 0; ширина ленты составляла 9 диагоналей, т. е. кроме главной диагонали имелось рядом по 4 ко-диагонали с каждой стороны), порядки которых поочередно выбирались равными 150, 200 и 250, в преобразованиях Хаусхолдера не сохранялся инвариантным «квадрат евклидовой нормы матриц» Б2. Для ленточной структуры с порядком матрицы N = 150 вместо Б2 = 0.1330-104 получилось Б2 = 0.3641769-104. Для N = 200 вместо Б2 = 0.1780 104 получилась величина Б2 = 0.18722 104 Аналогичная ситуация встречалась и для заполненных матриц порядка 250-300 со случайным распределением значений элементов матриц.

Для приведенных примеров ленточных матриц стандартный метод Гивенса и предлагаемая здесь модификация метода Гивенса обеспечивали корректные результаты (погрешность для Б2 возникала в 13-м и 14-м знаках после запятой).

Поэтому можно сделать вывод, что метод Гивенса и предлагаемая модификация метода Гивенса в аспекте обеспечения устойчивости и точности преобразований на практике предпочтительнее метода Хаусхолдера.

Обобщая результаты численных экспериментов по быстродействию упомянутых выше методов, можно заключить следующее:

1. Подтвердилась оценка повышения быстродействия модификации метода Гивенса по отношению к его стандартному варианту: быстродействие модификации метода Гивенса приблизительно в 1,4 раза выше быстродействия стандартного метода Гивенса.

2. Для ленточных (разреженных) матриц быстродействие метода Хаусхолдера в 1,0 ^ 1,08 раза

выше быстродействия модифицированного метода Гивенса.

3. Для заполненных матриц быстродействие метода Хаусхолдера в 1,1^ 1,15 раза выше быстродействия модифицированного метода Гивенса.

4. Для некоторых «сильно разреженных» матриц модифицированный метод Гивенса оказывается эффективнее по быстродействию метода Хаусхолдера в 1,1 раза.

Заключение. Таким образом, по результатам анализа точности и быстродействия можно сделать следующий вывод. В задачах диагонализации матриц на этапе приведения исходной матрицы А0 к компактной форме (форме Хессенберга или трехдиагональ-ной форме) использование модифицированного метода Гивенса может оказаться предпочтительнее в силу того, что в отличие от метода Хаусхолдера модификация метода Гивенса обеспечивает гарантированную точность и устойчивость преобразований и имеет с ним практически одинаковое быстродействие.

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

1. Уилкинсон Дж.Х. Алгебраическая проблема собственных значений. — М., 1970.

2. Иордан В.И. Модифицированный метод Гивенса // Рукопись представлена АлтПИ. Деп. в ВИНИТИ 27.06.88, № 5118-В88.

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

4. Парлетт Б. Симметричная проблема собственных значений. Численные методы. — М., 1983.