Аналогично, матрица C=:UiB отличается от матрицы Втолько элементами k-n и 1-й строк:

cji = bji при j фк, I и 1 t п.

(44)

Следовательно, матрица C = UAU отличается от матрицы А лишь двумя строками и двумя столбцами. Формулы для вычисления элементов этих строк и столбцов написать нетрудно, но в этом нет необходимости; удобнее программировать на ЭВМ непосредственно формулы (43)-(44). Заметим, что если матрица А эрмитова, то матрица С также будет эрмитова; тогда в изменившихся столбцах -И строках достаточно вычислить только половину элементов и тем самым вдвое уменьшить объем расчетов.

Найдем такую последовательность элементарных вращений, которая приводит произвольную (неэрмитову) матрицу А к верхней почти треугольной форме.

Можно так подобрать угол поворота в матрице Uki, чтобы уничтожить элемент С;, расположенный непосредственно перед левым нижним углом подматрицы плоского поворота (рис. 33, а). Из формул (43) и (44) видно, что для этого надо положить

Рис. 33.

V I а*, к-1 ? +1 «/, к-1 f

ai, h-i Ok, k-i

(45a)

Сами углы вычислять нет необходимости, ибо в формулы для преобразования матричных элементов они не входят. Отметим, что для вещественных матриц величина р тоже будет вещественной; тогда формулы (45) удобнее записать следующим образом:

«*, к-1

а/, к-1

Vl k-i+l к~1 /4, А-1

(456)

Теперь будем аннулировать те элементы матрицы и в том порядке, как это указано цифрами на рис. 33, б. Первый элемент уничтожается при помощи матрицы U23, обозначенной на рисунке сплошным квадратом. Второй уничтожается вращением и, обозначенным пунктирным квадратом. При втором вращении в матрице А меняются элементы вторых и четвертых строк и столбцов. Значит, аннулированный элемент «1», лежащий в трегьеп строке, так и останется равным нулю.

Продолжая эти рассуждения, можно убедиться, что однажды уничтоженный элемент при такой последовательности исключения будет оставаться равным нулю. Поэтому после окончания всех исключений матрица станет верхней почти треугольной матрицей (ау = 0 при {>/ + !). Это справедливо для произвольной -(неэрмитовой) матрицы.

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

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

Для полученной трехдиагональной (или почти треугольной) матрицы можно вычислять собственные значения и собственные векторы способами, изложенными в § 1. Найденные собственные значения будут одновременно собственными значениями исходной матрицы. А собственные векторы л,- исходной матрицы связаны с собственными векторами трехдиагональной матрицы соотношением

Xi = U3U2i ... Un-i.nyi- (46)

Проще всего вычислять их, последовательно умножая требуемый вектор у слева на матрицы вращения. Структура матриц такова, что при умножении на и меняются только k-я и /-я компоненты вектора

•Vft = a( ;-p*(/;, xi = yii + ayi,

х/ = у,- при jk, I.

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

На приведение эрмитовой матрицы к трехдиагональной форме и нахождение всех собственных значений в методе вращений требуется около 2rfi-{-b0n арифметических действий и ячеек оперативной памяти. Для нахождения каждого собственного вектора надо затратить еще 3ifi действий. Собственные значения и собственные векторы в этом методе определяются устойчиво (если унитарность U/ii не нарушена ошибками округления).

3. Итерационный метод вращений. Несмотря на свою быстроту, описанные выше прямые методы не вполне удовлетворительны. Так, их алгоритм состоит из разнородных частей: преобразования исходной матрицы, вычисления корней многочлена, нахождения собственных векторов обратными итерациями. Кроме того, их формулы не упрощаются для некоторых употребительных специальных форм матриц (например, ленточных); тем самым они невыгодны для таких матриц. Поэтому разработан и используется ряд итерационцых методов, в общем случае более медленных, но обладающих какими-то частными преимуществами.

Для эрмитовых матриц наиболее известен итерационный метод вращений, предложенный Якоби в 1846 г.; но в численных расчетах он начал использоваться только после появления работы [52]. Метод основан на подборе такой бесконечной последовательности элементарных вращений, которая в пределе преобразует эрмитову матрицу А в диагональную. При этом используются преобразования вращения с матрицами (42) такого же типа, как и для прямого метода вращений, но последовательность поворотов и их углы подбираются совершенно иным способом.

Рассмотрим, как действует элементарное вращение на сферическую норму матрицы (точнее, квадрат этой нормы):

5 = И= 2 ауР. (48)

*) Это является частным случаем общего утверждения, которое мы не доказываем: сферическая норма любой матрицы не меняется при умножении с любой стороны на унитарную матрицу.

Для определенности рассмотрим сначала умножение справа, В = = AUki. Из формулы (43) видно, что при этом элементы -го и /-Г0 столбцов меняются так, что попарные суммы квадратов модулей сохраняются:

элементы остальных столбцов остаются неизменными. Отсюда следует, что ЦЛС/Ц = ЦЛЦ, т. е. сферическая норма матрицы А не меняется при умножении справа на матрицу вращения. Аналогичное утверждение легко доказать для умножения на матрицу U или U" слева*). Описываемый метод основан на сохранении сферической нормы при вращениях.

Разобьем сумму, входящую в сферическую норму (48), на диагональную и недиагональную части:

5i= Ё \ап\\ S,= 21 lifl- (4)

i =1 г, / = 1

При элементарном преобразовании вращения UkiAUi недиагональные элементы а,-, ац и aki, ац при 1фк, I меняются так, что попарные суммы квадратов их модулей сохраняются; это легко видеть из формул (43)-(44). Кроме этих элементов вне диагонали есть еще один меняющийся элемент -это а. Поэтому величина меняется при элементарном вращении настолько, насколько изменится 1а«. Будем подбирать вращения так, чтобы уменьшалась.

Чтобы максимально уменьшить за одно вращение, подберем угол поворота так, чтобы аннулировать элемент at. Для простоты ограничимся вещественными эрмитовыми (т. е. симметричными) матрицами. Тогда = а;* - вещественные числа, и матрицы вращений Ui тоже вещественны. Из формул (44) и (43) с учетом вещественности всех величин следует

Сы = bkia -f ЬггР = а {а - -f (а„ - а**) ар.

Полагая Cw = 0 и вспоминая условие нормировки (42), получим систему уравнений для определения параметров поворота

+ (50)

аы{а-) = {аии-аи)а.

Возводя второе уравнение (50) в квадрат и исключая из него р при помощи первого уравнения, получим биквадратное уравнение