Дж. Деммель - Вычислительная линейная алгебра, страница 41

3. О сходимости судят по «достаточной малости» подциагональных элементов матрицы А,. Чтобы найти практичный критерий «малости», воспользуемся понятием обратной устойчивости. Поскольку А; получается из А ортогональным подобием, в любом случае следует ожидать, что А; включает в себя ошибки округлений порядка 0(г11А11). Таким образом, всякий поддиагональный элемент в А;, по абсолютной величине меньший, чем 0(г~)А~!), вполне мог бы быть нулем, поэтому мы и заменяем его на нуль'. Если А— верхняя хессенбергова матрица, то замена элемента ар+цр нулем превращает А в верхнюю блочно-треугольную матрицу А = ~ 4, где ~ Аы Агг 22 Ап — размера р х р и обе подматрицы Аы и А22 — верхние хессенберговы. После этого собственные значения матрицы А могут быть получены путем независимого вычисления собственных значений для Ам и А22. Когда все диагональные блоки станут размера 1 х 1 или 2 х 2, алгоритм закончит работу.

4.4.6. Приведение к форме Хессенберга Для заданной вещественной матрицы А нужно найти такую ортогональную матрицу Я, чтобы ЯАцг была верхней хессенберговой матрицей. Соответ- г Нв практике используется несколько более жесткий критерий, в котором вместо ~бА~2 берется норма соответствующеМ подматрицы в А. Он учитывает наличие твк называемых «гргдуироввнных» матриц с большим разбросом в величинах элементов. Ноддизгонвльный элемент может быть заменен нулем и в том случае, когда достаточно мало произведение ар»» рог»2, »г двух последовательных подднзгонвльных элементов. Детали можно найти в ЬАРАСК-прогрвмме я!вьцг.

4.4. Алгоритмы для несимметричной проблемы собственных значений 177 ствующий алгоритм является простой вариацией идеи, использованной для (1Е-разложения. 1. Выбрать Щ так, чтобы х х х х х х х х х х х х х х х о х х х х о х х х х о х х х х х х х о х х х х х х и А =а Аат= 1„1А = х х х х о х х о х х Умножение на ц1 не изменяет первую строку в А, а умножение на Ц1 не т изменяет первый столбец в 1,1А, сохраняя полученные в нем нули. 2. Выбрат/ цз так, чтобы х х х х х х х х х х х х х х х х х х х х о х х х х х х И А2 = Я2А1Ч2 = Ю2'11— о х х о о х х х о о х о о х х х х х о о х х х Умножение на Яз изменяет лишь три нижние строки в А1, а умножение на Чзт не изменяет первые два столбца в ЩА1, сохраняя имеющиеся в них нули. 3.

Выбрать 1гз так, чтобы х х х т, х х х х х х х х х х х о х х х х х х х х х о х х о о х и Аз = С/зА2С/з —— сгзА2 = х х о о х о о о х х х х х х о о о х х МатРиЦа Аз — веРхнЯЯ хессенбеРгова. В Целом, имеем Аз = (ЯзЯ291) . АЯ2112Я1)т = ЯАЯт сг Общий алгоритм для приведения к форме Хессенберга выглядит так: Алгоритм 4.6. Приведение к верхней форме Хессенбергаг если нужна матрица Я, то положить Я = Х /ог 1' = 1: п — 2 и, = Нопзе(А(2 + 1: п,з)) Рз = 7 — 2и;и1 /' Я1 = с)1ак(Рх' Р.) '/ А(2+ 1: п,з': п) = Рз А(1+ 1: п,з; п) А(1: п,з+ 1: п) = А(1: п,з+ 1: п) Р, если нужна матприца с/ Я(2'+1:п,з:п)=Р, Я(1'+1:п,з':п) /*Я=Щ Я*/ епа 1/ епа /ог Пример 4.8.

Общую последовательность действий в приведении к форме Хес- сенберга проиллюстрируем с помощью примера 5-го порядка. В нем Я, озна- чает отражение порядка 5, выбранное так, чтобы аннулировать в 1-м столбце матрицы элементы 1 + 2,..., п, не меняя элементов 1,..., 2. 178 Глава 4. Несимметричиав проблема собственных значений Как и в случае ЯК-разложения, матрицу Р; не формируют в явном виде: умножение на 1 — 2иги~ производится посредством матрично-векторных операций. Векторы и; могут храниться в подциагональной части матрицы, снова по аналогии с 11В;разложением. Отражения могут применяться с использованием ВТАБ-операций уровня 3, как это описано в вопросе 3.17.

Алгоритм приведения к форме Хессенберга реализован в Ма11аЬ'е командой Ьевв, а в 1 АРАСК'е — программой вйейтб. Легко подсчитывается, что число флопов в этом алгоритме равно з и + 0(пг), а если вычисляется и произведение О = ь1„1...Яы то ~~п + 0(п~). Преимущество формы Хессенберга для Яй-итерации заключается в том, что один шаг для нее стоит только бпз + 0(п) флопов вместо 0(пз) в общем случае. Кроме того, эта форма сохраняется алгоритмом, так что матрица остается верхней хессенберговой. Предложение 4.5.

Форма Хессенберга сохраняется ЯН-итерацией. Доказательство. Легко убедиться в том, что при 14й.-разложении верхней хессенберговой матрицы (например, матрицы А, — о.1) получается верхняя хессенбергова матрица Я (потому что у-й столбец в Я есть линейная комбинация первых з столбцов матрицы А, — о1). После этого легко проверить, что матрица Ргц остается верхней хессенберговой. Добавление о1 не изменяет формы матрицы. П Определение 4.5. Вер яв хессенбергова матрица Н называется неразложимой, если все злементьг на ее первой поддиагоноли отличны от нуля. Легко видеть, что если 6, ь13 — — О, т.

е. матрица Н разложима, то ее собственными значениями являются собственные значения ведущей главной подматрицы порядка г и дополнительной подматрицы порядка и — 1; обе они хессенберговы. Итак, во всех случаях можно ограничиться рассмотрением неразложимых матриц. 4.4.7. Приведение к трехднвгоналъной н двухднвгональной формам Если матрица А симметрична, то каждый полный шаг приведения к форме Хессенберга дает симметричную матрицу, поэтому нули создаются в симметричных позициях.

Это означает, что нужно обрабатывать лить половину матрицы; в результате число операций снижается до 4~пз + 0(пз), а если формируется матрица Я„1... Яы то до впз + 0(пз). Такой алгоритм называется приведением к трехдиагональной форме. Мы используем его в гл. 5.

Алгоритм реализован ЬАРАСК-программой ввустб. Несколько предвосхищая обсуждение вопроса о вычислении БУР (см. равд. 5.4), напомним (равд. 3.2.3), что собственные значения симметричной матрицы АтА суть квадраты сингулярных чисел матрицы А. Наш окончательный ЯУй-алгоритм будет опираться на этот факт, поэтому хотелось бы найти для А форму, которая обеспечив гла бы трехдиагонвльность матрицы АтА. В качестве такой формы возьмем верхнюю двухдиагональную матрипу, т.е.

матрицу, ненулевые элементы которой могут располагаться только на главной диагонали и первой нвддиагонали. Итак, требуется вычислить ортогональные матрицы Ц 4.4. Алгоритмы для несимметричной проблемы собственных значений 179 и г' таким образом, чтобы матрица цА г' была верхней двухдиагональной. Соответствующий алгоритм называется приведением к двухдиагональной форме и очень похож на приведение к форме Хессенберга и трехдиагональной форме. Пример 4.9.

Проиллюстрируем общую последовательность действий в при- ведении к двухдиагональной форме с помощью примера 4-го порядка: 1. Выбрать Ог так, чтобы х х о о о х х х о х х х о х х х :] и $'г так, чтобы Аг = ЯгАЪ'г —— Ц,А= [ Матрицы Яг и Ъг суть отражения, причем гг не меняет первый столбец матрицы Яг А. 2. Выбрать Яг так, чтобы о о х х о о ] х х и Ъг так, чтобы Аг х х х х х х о х ЮгАг= о о = ЯгАг)'г= Здесь с7г — отражение, не изменяющее первую строку в Аы а Иг — отраже- ние, не меняющее первые два столбца матрицы ЯгАы 3. Выбрать Яз так, чтобы х х о о о х х о о о х х о о о х ЯзАг = и положить Ъз — — Х, так что Аз = ЯзАг. Отражение цз не меняет первые две строки матрицы Аг. 4.4.8.

ЯК-итерация с неявными сдвигами В этом разделе будет показано, как эффективно реализовать ЯК-итерацию для верхней хессенберговой матрицы Н. Эта реализация будет неявной в том смысле, что ЯК-разложение матрицы Н не вычисляется в явном виде. Матрица Ч строится неявно как произведение вращений и других простых ортогональных матриц. Излагаемая ниже неявная Я-пгеорема показывает, что эта неявно построенная матрица Ч' и есть та ортогональная матрица, которая вам нужна. Далее мы указываем, как внедрить в нашу процедуру сдвиг и, В общем случае и х п-матрицы А будут получены ортогональные матрицы м' = ~,г -г...Яг и И = 1'г... 1~„-г такие, что матрица ЯАЪ' = А' — верхняя двухдиагональная.

Заметим, что А'"А' = ИгАг ЯтЯАЪ' = Иг Аг АИ, т. е. А'г А' имеет те же собственные значения, что и АтА, а потому А' имеет те же сингулярные числа, что и А. Стоимость приведения к двухдиагональной форме составляет зпа + 0(пг) флопов, а если вычисляются матрицы Я и И, то 4пз + 0(пг) флопов. Этот алгоритм реализован ЬАРАСК-программой вбебгд. 180 Глава 4. Несимметричная проблема собственных значений необходимый для ускорения сходимости.

Затем мы обсуждаем прием двойного сдвига, используемый для того, чтобы проводить вычисления в вещественной арифметике да.ке в присутствии комплексных собственных значений. Он состоит в комбинировании двух последовательных ЯГь-итераций с комплексно- сопряженными сдвигами о и о; матрица после этого двойного сдвига снова становится вещественной. Наконец, мы обсудим стратегии выбора сдвигов о и о, надежно обеспечивающие квадратичную сходимость. Недавно были обнаружены некоторые (редкие) ситуации, когда алгоритм не сходится ~25, 65].