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

Поэтому задача построения полностью надежного и быстрого варианта ЯГситерации остается открытой. Неявная г1-теорема Наша окончательная реализация цГь-итерации будет опираться на следующую теорему. Теорема 4.9 (неявная Я-теорема). Пусть ЯтАЯ = Н вЂ” неразложимая верхняя матрица Хгссенбгрга. Тогда первый столбец матрицы Я однозначно (с точносгпью до знаков) определяет ге г-й,..., и-й столбцы. Из этой теоремы следует: чтобы в цВ;алгоритме вычислить по А; матрицу Аг»-1 = Я, Аг ги достаточно 1, вычислить первый столбец матрицы Я, (~соторый параллелен первому столбцу матрицы А, — о,1 и может быть получен из него всего лишь нормировкой) и 2.

выбрать остальные столбцы матрицы ц; так, чтобы ц; была ортогональна, а Аыа была неразложимой хессенберговой матрицей. В этом случае неявная Я-теорема гарантирует, что мы правильно вычислим матрицу А;»и, поскольку Я; определена однозначно с точностью до знаков столбцов, а зти знаки не существенны. (Изменение знаков столбцов в ьЪ равносильно замене разложения А; — о;1 = Я;Нг на (Я;Яг)(Я,Н;), где Я; = дйа8(~1,..., х1). Но тогда А +1 — — (Я,Н,) ОЩЯ») + о 1 = Я,(Н,ьЪ + о 1) $о Это ортогональное подобие всего лишь изменяет знаки строк и столбцов в Аььь.) Доказательство неявной Я-теоремы.

Предположим, что ЯтАЯ = Н и Из АИ = С вЂ” неразложимые верхние матрицы Хессенберга, причем матрицы ц и И ортогональны, а их первые столбцы равны. Примем обозначение 1Х), для 1-го столбца матрицы Х. Мы хотим показать, что Я)г = х1$'); для всех г' > 1 или, иначе говоря, что Иг = г'~Я = 61а81х1,..., х1). Поскольку И' = РтЯ, имеем СИт = Сг'тЯ = 'г'тАЯ = 'г'тЯН ИгН. Из равенства СИ' = И'Н выводим С(Ит)г = 1СИ'); = (ИгН).

~,, и.;1И')гч откуда )г;» 1,,(И'),+1 — — С1И'); — ~ 1, Ь1;1И') .. Учитывая, что (Ит)1 — — 11,0,..., 0)т, а С имеет верхнюю форму Хессенберга, можем индукцией по г' показать, что в столбце (Иг), поддиагональные элементы равны нулю, т.е. матрица И' — верхняя треугольная. Так как И' в то же время ортогональна, то Иг — диагональная матрица: И' = 61а81х1,..., х1). П 4.4.

Алгоритмы для несимметричной проблемы собственных значений 181 Неявный ь)К-алгоритм с одинарным сдвигом Покажем с помощью примера 5-го порядка, как использовать ц-теорему для вычисления матрицы Аг по Ае = А. х х х х х х х х х х + х х х х о о х х х о о о х х гмт так, чтобы Аг = Я~АЯг —— Выбор сг и зг мы обсудим позже. Пока что эти параметры могут соответствовать произвольному вращению.

Символ + в позиции (3, 1) называется емспгупом (или горбом). Мы должны избавиться от него, чтобы восстановить форму Хессенберга. 2. Выбрать х х х х х сг зг — зг сг х х х о х х х х х х т 'ег так, чтобы Я~~Аг —— о о х х х о о о х х Аг = Яг АгЯг = Таким образом, выступ «вытеснен» из позиции (3, 1) в позицию (4, 2). 3. Выбрать х х х х х х х х о х х х х х х Ст так чтобы Яз Аг —— т сз зз — ез сз 1 о о х о о о х х х х Аз — = Яз АгЮз = Пример 4.10.

1. Выбрать сг ег — зг сг 1 1 1 Выступ вытеснен из позиции (4, 2) в позицию (5, 3). х х х х х х о х х о + х о о о х х х х х х о х х о о х о о + х х х х х х х х х х х х х х х х х х х х 182 Глава 4. Несимметричная проблема собственных значений 4. Выбрать 1 1 1 х х х х х х о х х х х х х так, чтобы 414 Аз = т х х х х о о х о о о С4 З4 — З4 С4 х х х х х х х х х х х х 14 = Ю4 Аз'44 т о х х о о х о о о х х х х х х Таким образом, форма Хессенберга восстановлена.

В целом, имеем: для матрицы С1 Х Х В1 Х Х зз х зз х х х х х х х х В4 Х 14 — а 1 а 2 а 3 ьг 4— матрица Яз АЯ вЂ” верхняя хессенбергова. По неявной Я-теореме, первый столбец [с1, з1,0,..., 0]~ матрицы Я однозначно (с точностью до х1) определяет остальные столбцы. Выберем первый столбец для Я так, чтобы он был пропорционален первому столбцу матрицы А — о Х, т. е. вектору [аы — о, аз1, О,..., 0]т. Тогда Я совпадет с ортогональной матрицей в ЯГс-разложении А — о1, что и требовалось. О Неявный 44К-алгоритм с двойным сдвигом В этом разделе будет показано, как, выполняя сдвиги и и о совместно, сохра- нить возможность вычислений в вещественной арифметике. Это существенно для эффективности практической реализации алгоритма,но не для понима- ния его математики.

Поэтому при первом чтении книги данный раздел можно опустить. Результаты последовательного использования сдвигов и и о можно описать формулами Ао — о1 = 4г1Н1, А, = НЯ1+ о1, так что А1 —— Я1~АсЯ1, А1 — 1тХ = ЯзЛз, Аз = НзЯз+ оХ, так что Аз — — Яз~А19з = Яз~А1 АсЯ1Яз. Лемма 4.5. Матрицы Щ и Яз мозсно вь1брать так, чтобы: (1) произведение сглаз бьто вещественным и, следовательно, (2) матрица Аз была вегцественной а (3) первый столбец матрицы Я1Яз легко вычислллсл.

Для п хи-матрицы стоимость одного шага неявной ЯН-итерации составляет бпз + 0(п) флопов. 4.4. Алгоритмы для несимметричной проблемы собственных значений 183 аг„+ а»гам — 2(Яо)аы + )о(2 агг(аы + агг — 2(3«о)) амагг О Остальные столбцы матрицы Я»Я2 вычисляются неявно в соответствии с неявной СХ-теоремой. Процесс по-прежнему называется «вытеснением выступа», только выступ теперь размера 2 х 2, а не 1 х 1.

Пример 4.11. Проиллюстрируем вытеснение выступа с помощью примера б- го порядка. т Хет О г 1. Положим Ятг = ' ~, где первый столбец матрицы Щ определяется О приведенной выше формулой. Тогда х х х х х х х х х х х х х х х х х х + + о о х х х х + х х х о о х х о о о х х х х х х х х + х х о о х х х х х и А,=Ят»АЯ дтА— х х х х х х о о о о х х о о о х х Заметим, что появился выступ размера 2 х 2, указываемый плюсами. 2. Выберем отражение Я~т, которое, воздействуя только на 2-ю, 3-ю и 4-ю строки матрицы Аы аннулирует ее элементы (3, 1) и (4, 1) (это означает, что вне подматрицы, образованной строками и столбцами 2, 3 и 4, Ягт совпадает с единичной матрицей): х х х х х х о х х о + х + + о о о х х х х х х х х х х х х х х х о х х х х х х х х х х х х х х х х + х х х о о х х и Аг = Д2~А»Ц2 —— ««'г Аг= т о о о о х В результате 2 х 2-выступ «вытеснен» на один столбец вправо.

Док эательсгпео. Поскольку Ягйг = Аг — й1 = Лэсли + (о — й)Х, имеем »ег»ег«»2«»1»«1(Я1»»'1 + (о й)1)«»1 = Я»В»Я»В» + (и — й)Я»В» = (Ао — о.1)(Ао — о1) + (а — й)(Ае — оХ) = Аго 2(Яо)Ао + )о~21 М Таким образом, произведение (Ч»Я2)(В2В») есть ЯВ;разложение вещественной матрицы М. Поэтому матрицы чгцг и Вгйг могут быть выбраны вещественными. Но тогда и матрица Аг = ЯЩ2)г А(4ХЯг) вещественна. Первый столбец в Я»Я2 пропорционален первому столбцу матрицы Агс— 2%аАо + (о~21, который равен 184 Глава 4.

Несимметричная проблема собственных значений 3. Выберем отражение Я~т, которое, воздействуя только на З-ю, 4-ю и 5-ю строки матрицы А2, аннулирует ее элементы (4, 2) и (5, 2) (это означает, что вне подматрицы, образованной строками и столбцами 3, 4 и 5, Яз совпадает т с единичной матрицей): х х х х х х х х х х х х х х х х х х х х х х х х х х х х х х о х х х о о х х о о + х о о о о о х х х о о х х о о + х о о + + х х и Аз=Ма А2Яз = т Оз А2— т х х х х х х 4.

Выберем отражение (~~е, которое, воздействуя только на 4-ю, 5-ю и 6-ю строки матрицы Аз, аннулирует ее элементы (5, 3) и (6, 3) (это означает, что вне подматрицы, образованной строками и столбцами 4, 5 и 6, Я4 совпадает с единичной матрицей): х х х х х х х х х х х х о х' о о о о о о х х х х х х х х о х х х о + х х А4 = Я4АзА4 = 5. Выберем матрицу 1 1 1 1 с я — ее х х х х х х х х х х х х х х о х х х так, чтобы Аз = Яз А4Яв = т о о х х х х о о о х х х о о о о х х Выбор сдвига в ЯК-алгоритме Чтобы полностью охарактеризовать шаг хессенберговой (4Гс-итерации с одинарным или двойным сдвигом, мы должны указать способ выбора сдвига а (и а).

Вспомним (см. конец раздела 4.4.4), что разумным выбором значения одинарного сдвига, асимптотически обеспечивающим квадратичную сходимость к вещественному собственному значению, является правило а = а„„, т. е. значение нижнего углового элемента матрицы Ао Обобщением этого правила на случай двойного сдвига является правило 42рянсиса, Оно означает, что и и и суть собственные значения нижней угловой 2 х 2-подматрицы ~ а„,„1 аи,и матрицы А,. Это правило позволяет нам получить сходимость как к двум вещественным собственным значениям в нижнем 2 х 2-углу матрицы, так и к 2 х 2- блоку с комплексно-сопряженными собственными значениями.

Приближение к моменту, когда сходнмость наступила, проявляется в том, что элемент а„г „2 (а возможно, и а„„1) становится мал; в этом случае собственные значения 185 4.5. Другие типы иасииматричиых спектральных задач угловой 2 х 2-подматрицы являются хорошими приближениями к собственным значениям самой матрицы А. Действительно, можно показать, что правило Фрэнсиса асимптотически обеспечивает квадратичную сходимость.

Это значит, что если элемент а„» „г (а возможно, и а„,„1) уже достаточно мал, то с каждым следующим шагом его абсолютная величина возводится в квадрат и быстро стремится к нулю. На практике этот прием работает столь хорошо, что для почти всех матриц требуется в среднем два «гВ;шага для сходимости одного собственного значения. Это оправдывает трактовку ОВ;итерации как «прямого» метода.

Существуют примеры, для которых «4В:итерация со сдвигом Фрэнсиса не сходится в определенном выше смысле (например, матрица остается в этом алгоритме неизменной). Поэтому практические реализации алгоритма, используемые в течение десятилетий, включают в себя так называемый «особый сдвиг»; он применяется в том случае, если предыдущие 10 обычных сдвигов не привели к сходимости. Тем не менее недавно были обнаружены малые множества матриц, для которых не сходятся и эти реализации [25, 65].