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

Эти явления ложной сходимости объясняют важность вычислимой оценки погрешности из утверждения 2 теоремы 7.2 для контроля сходимости собственных значений [198]. Если оценка погрешности мала, то вычисленное собственное значение действительно является хорошим приближением к некоторому точному собственному значению, даже если какое-то другое оказалось «пропущенным». О 7.4. Алгоритм Ланцоша а арифметике с плавающей точкой 391 Имеется еще одна оценка погрешности, принадлежащая Каниэлю и Сааду, которая проливает свет на причины ложной сходимости.

Эта оценка .: етавлена в терминах углов между начальным вектором уг и нужными собственными векторами, а также чисел Ритца и нужных собственных значений. Другими словами, в ней участвуют величины, неизвестные в ходе вычислений, поэтому с практической точки зрения она бесполезна. Однако эта оценка показывает, что если у«почти ортогонален к разыскиваемому собственному вектору или же разыскиваемое собственное значение почти кратно, то следует ожидать медленной сходимости. По поводу деталей см. [197, разд. 12-4]. 7.4.

Алгоритм Лаицоша в арифметике с плавающей точкой Пример, разобранный в предыдущем разделе, описывал поведение «идеального» алгоритма Ланцоша, когда округления по существу отсутствуют. Мы назвали тщательную, но дорогостоящую реализацию алгоритма 6.10 л«етодол«Ланцоша с полной переортогонализацией, противопоставляя ее исходной необременительной реализации, называемой методом Ланцоша без переортогонализации (см. НОМЕРАСЕ/Мас)аЬ/1 апсзоеХоиеогФЬой.ш). Приведем запись, объединяющую оба этих метода. Алгоритм 7.2.

Алгоритм Ланцоша (соответственно с полной переортогонализацией и без переортогонализации) для вычисления собственных значений и собственных векторов матрицы А = Ат « у« = («/][(«[]г, 17о = О, уо = О /агу =1 сок г = Ау. о« = у. г « < г = г — ~ ««(г~уе)у;, г = г — Я~~ «(г~уг)уг полная переортогонализация г««' у«' /««'-«у«' — « переортогонализация отсутствует А = ][г[[г если 3« = О, то прекратить выполнение лгоритма у«~- = /«1« Вычислить собственные значения и собственные векторы матрицы Т и оценки погрешностей в них епа /от Полная переортогонализация соответствует повгпорному проведению процесса ортогонализации Грама — Шмидта «г = г — ~ «,(г~уе)уг» для того, чтобы почти гарантировать, что г будет ортогонален векторам уы ..., у« (См.

алгоритм 3.1, а также [197, разд. 6-9] и [171, гл. 7] по поводу обсуждений вопроса о том, когда «двух ортогонализаций достаточно».) В равд. 6.6.1 мы показали, что в точной арифметике г ортогонален уы ..., 91 «без какой-либо переортогонализации. К сожалению, это свойство ортогональности, на котором до сих пор основывался весь наш анализ, теряется из-за округлений. 392 Глава 7. Итерационные методы для задач иа собственные значения Потеря ортогональности не заставляет алгоритм вести себя совершенно непредсказуемым образом. В самом деле, мы увидим, что цена, которую мы платим за потерю, — это приобретение повтиоряыт копий сошедшихся чисел Ритца.

Иными словами, для больших 7г матрица Тк может иметь не одно собственное значение, очень близкое к А,(А), но много таких собственных значений. Это не является катастрофой, если нас не интересуют кратности собственных значений и если замедление сходимости внутренних собственных значений, вызванное появлением копий крайних собственных чисел, несущественно для нас. Детальное описание реализации алгоритма Ланцоша, оперирующей именно таким образом, можно найти в [57], а сами программы — в ХЕТ11В/1апсгоз. (В этих программах используются некоторые эвристики для оценки кратностей собственных значений.) Однако если знание точных кратностей является важным, то нужно поддерживать (почти) ортогональность векторов Ланцоша.

В таком случае можно прибегнуть к алгоритму Ланцоша с полной переортогонализацией, как было сделано в предыдущем разделе. Легко проверить, правда, что 7с шагов этого алгоритма потребуют 0(йзп) флопов вместо 0(йп) и 0(7сп) слов памяти вместо 0(п). Такая цена может оказаться слишком высокой. К счастью, имеется вариант, промежуточный между полной переортогонализэцией и отсутствием всякой переортогонализации; он соединяет в себе почти в полном объеме лучшие качества этих двух крайних возможностей. Оказывается, что векторы 9ь теряют ортогональность весьма систематическим образом, а именно, приобретая большие компоненты в направлениях уже сошедшихся векторов Ритца.

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

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

Пример 7.2. Для матрицы примера 7.1 были проведены 146 шагов Ланцоша. Картина сходимости метода представлена на рис. 7.7. Диаграммы правого ряда соответствуют методу с полной переортогонализацией, а диаграммы левого ряда — методу без переортогонализацин. Эти диаграммы аналогичны диаграммам рис. 7.3 с тем отличием, что удалены графики глобальной ошибки, загромождающие средние диаграммы. На рис. 7.6 наименьшее сингулярное число о ыЯг) изображено как функция от номера шага й. В точной арифметике матрица Щ имела бы ортонормированные столбцы, а потому было бы о ыЯь) = 1. В арифметике с округлениями Яь теряет ортогональность столбцов, начиная примерно с 7с = 70. К 80-му шагу значение и маг) падает до уровня 0.01, и именно с этого момента две верхние диаграммы рис.

7.7 начинают заметно различаться. 7.4. Алгоритм Ланцоша в ариг(>метине с плавающей точкой 393 Наименьшее сингулярное число матрицы, составленной иа первых к векторов Ланцоша ю' га ! гв" О ю ю' М и ~ ге о о ге" 8В ге ~ш гс" го" а ю м ю ю гю кю !40 Номер шага (без переоргогоншгнзапнн) Рис. 7.6.

Алгоритм Ланцоша без переортогоналнзацни в применении к матрице А. Для й от 1 до 149 показано наименьшее сингулярное число !г !„(Яа) матрицы (Г)а), составленной из векторов Ланцоша. В отсутствие округлений столбцы матрицы Щ ортонормированы, и все сингулярные числа равны единице. При наличии округлений ГЕ» теряет свойство полноты ранга. В частности, начиная с й = 80, второе по величине (красное) собственное значение Ла(Та), которое уже сошлось к Ла(А) — 2.7 почти в 16 разрядах, вдруг в несколько шагов поднимается до уровня Лг(А) ш 2.81 (см, верхнюю левую диаграмму рис. 7.7).

В результате получается «вторая копия» числа Лг(А); первую копию дает Лг(Та), которому соответствует черный цвет. (На рисунке эти черные плюсы трудно разглядеть, потому что прямо на них нанесены красные.) Это «превращение» числа Ла(Та) сопровождается прыжком оценки ошибки (пунктирная красная линия на левой средней диаграмме). Кроме того, этому превращению предшествует возрастание компоненты по направлению ег на левой нижней диаграмме: на шаге 50 черная кривая начинает расти вместо того, чтобы, продолжая убывать, снизиться до уровня машинного эпсилон, как это происходит при полной переортогонализации на правой нижней диаграмме.

И то, и другое свидетельствует о том, что алгоритм отклоняется от своего точного пути (и что, следовательно, необходима некоторая выборочная ортогонализация). Вслед за завершением сходимости второй копии числа Лг (А) начинает снова убывать компонента векторов Ланцоша в направлении ег (этот момент наступает вскоре после )с = 80).

Аналогично, начиная примерно с )с = 95 появляется вторая копия числа Лз(А): синяя кривая (соответствующая Л«(Та)) на левой верхней диаграмме перемещается с уровня Лз(А) 2.6 на уровень Ла(А) ке 2.7. В этот момент мы имеем две копии для Лг(А) 2.81 и две копии для Лз(А). Это не так легко увидеть из диаграмм, поскольку плюсы одного цвета затеняют плюсы другого (красный цвет поверх черного и синий поверх зеленого). Превращение сопровождается резким прыжком оценки ошибки для Л«(Та) (пунктирная 304 Глава 7.