3. Итерационный метод вращений.

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

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

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

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

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

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

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

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

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

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

Можно выбрать любой из четырех корней этого уравнения, тогда определится однозначно. Для определенности положим

Сами углы поворота находить не требуется.

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

Какой именно недиагональный элемент целесообразно аннулировать при очередном повороте? Конечно, если уничтожать максимальный по модулю внедиагональный элемент, то скорость убывания будет наибольшей. При ручных расчетах это наилучший способ. Но на ЭВМ перебор элементов матрицы для определения максимального элемента требует неприемлемо большого числа действий. А если Аннулировать элементы в заранее определенном порядке — циклом, то сходимость будет очень медленной. Наиболее выгодным оказалось уничтожение так называемого оптимального элемента.

Составим суммы квадратов модулей внедиагональных элементов строк:

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

Суммы (52) вычисляются тоже экономично, ибо при каждом вращении из них меняются только две и причем их можно вычислять по таким формулам:

штрихи относятся к значениям после вращений.

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

и тем самым стремится к нулю при Следовательно, процесс Якоби с выбором оптимального (и тем более максимального) элемента всегда сходится.

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

Итак, вдали от решения сходимость не хуже линейной, а вблизи решения - квадратичная, т. е. быстрая. Это позволяет получать все собственные значения с высокой точностью. Обычно процесс сходится за 6—8 циклов вращений, или за элементарных вращений. Интересно, что чем больше кратных собственных значений, тем быстрей сходится метод.

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

Но чтобы получить с той же точностью собственные векторы, надо или вычислять их по формулам

или делать еще один цикл вращений.

Хотя теоретически в методе Якоби могут накапливаться ошибки, фактически устойчивость и точность очень высоки. Кратные собственные значения получаются столь же точно, как и простые, а собственные векторы практически ортогональны друг другу.

Рис. 34.

Метод Якоби с выбором оптимального элемента требует обычно около арифметических действий и ячеек памяти для нахождения всех собственных значений. Для нахождения всех собственных векторов требуется еще около действий. Таким образом, этот метод раз в 10 медленнее метода отражений. Основное его достоинство — надежность и единообразность вычислений, что позволяет легко запрограммировать метод. Итерационный метод вращений применяется там, где важна точность, надежность и простота расчета и менее существен объем вычислений.

Замечание. Этот метод можно ускорить в 1,5 — 2 раза, не теряя его достоинств. У произвольных матриц недиагональная часть сферической нормы в среднем много больше диагональной, , а у трехдиагональных в среднем . Значит, трехдиагональная матрица является выгодным начальным приближением для итерационного метода Якоби. Поэтому целесообразно предварительно привести исходную эрмитову матрицу к трехдиагональной форме при помощи прямого метода вращений и затем первым ходом метода Якоби аннулировать все нечетные или все четные поддиагональные элемецты (рис. 34). После этого можно переходить на обычный вариант итерационного метода вращений с выбором оптимального элемента.