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

В результате память, необходимая для хранения Я, снижается с пз машинных слов до 0(п 1об п). Итак, на выходе алгоритма получаем матрицу 1»' в «факторизованном» виде, представленном множителями Я' для всех узлов дерева. Такое представление Я вполне достаточно для наших целей, поскольку позволяет умножить ее на любой вектор за время 0(п 1ой" и), используя метод РММ.

5.3.4. Бисекция и обратная итерация В алгоритме бисекции теорема Сильвестра об инерции (теорема 5.3) используется для вычисления только к желаемых собственных значений (вместо всех п) с затратой 0(п/с) операций. Напомним, что 1пег11а(А) = (и, ~, я), где и, ~ и и дают, соответственно, число отрицательных, нулевых и положительных собственных значений матрицы А. Если Х вЂ” невырожденная матрица, то, согласно теореме Сильвестра, 1пег11а(А) = 1пегИа(ХтАХ). Предположим теперь, что посредством гауссова исключения найдена факторизация А — вЕ = йР1т, где Ь вЂ” невырожденная матрица, а Р диагональная. Тогда 1пегба(А — в1) = 1пег11а(Р), а инерция диагональной матрицы Р находится тривиально.

(Ниже мы используем обозначение «1~ Ни < О» для количества чисел 4,, меньших нуля, и другие обозначения того же типа.) 1пег11а(А — в1) =(З]Г «1п <0,4 би =0,1В «(н )О) = (количество отрицательных собственных значений матрицы А — в1, количество нулевых собственных значений матрицы А — в1, количество положительных собственных 5.3. Алгоритмы длн симметричной проблемы собственных значений 241 значений матрицы А — г1) = (количество собственных значений матрицыА, меньших, чем г; количество собственных значений, равных г; количество собственных значений, ббльших, чем г). Предположим, что гь < гг, и пусть найдены 1пегба(А — гь1) и 1пегМа(А— гг1). Тогда число собственных значений матрицы А, попадающих в полуинтервал [вы гг), определяется формулой: (число собственных значений, меньших, чем гг) — (число собственных значений, меньших, чем гь).

Чтобы превратить это наблюдение в алгоритм, определим величину ХейсоипФ (А, г) = количество собственных значений матрицы А, меньших, чем г. Алгоритм 5.4. Бисекц я: найти все собственные значения матрицы А, принадлежащие полуинтервалу [а, Ь) с заданным допуском на ошибку ьо1ь п, = Медсоип1(А, а) пь = Федсоипг(А, Ь) если п = пь, то выход из алгоритма, так как в [а, Ь) нет собственных значений матрицы А поместить [а,п, Ь,пь] в ИгогЫ1з! /» %огЫьзг — это список полуинтервалов [а, Ь), содержащих собственные значения с номерами от п — п до и — пь+ 1. Эти полуинтервалы в алгоритме делятся пополам до тех пор, пока их длина не станет меньше допуска со1.

»/ ьвйь1е КогЫьзг непусто до удалить [1оьр,пм,пр,п„р] из И'огЫьзг ь1(пр — 1он ( 1о1) йеп вььвести сообщение»в полуинтервале [!очг, пр) находятся п„р — тч собственных значений» е1зе пзЫ = (1он+ пр)/2 п га — — 1»" едсоип1(А, шЫ) ьз'пым > п! йеп ... в [!оы,шЫ) есть собственные значения поместить [1оьч,пььч,шЫ,п ы] в ИгогЫьз! епд ь1 ь1п„р > п, ы йеп... в [шЫ,ир) есть собственные значения поместить [шЫ,п м,пр,п„р] в И»огЫьзг епд ь1 епд ь1 епд ьлйь!е Пусть оь » . о„— собственные значения матрицы А. Изложенную идею можно использовать и для определения чисел од с номерами ь' = зь,зо + 1,... 1ь.

Действительно, мы знаем, что в полуинтервале [!орь, цр) находятся собственные значения с номерами от о„„,.„до сь„„„,ть. 242 Глава 5. Симметричная проблема собственных значений Если А — плотная матрица, то число Хейсопп1(А, з) можно определить с помощью симметричного гауссова исключения с выбором главного элемента, как это описано в равд. 2.7.2. Однако такой метод не эффективен, поскольку требует 0(нз) операций для каждого значения г. В то же время, для трехдиагональной А число г1ейсоппс(А, г) найти весьма просто, выполняя разложение без выбора главного элемента: ь1 а1 — г ь аг г 707т А — з7 = ь„ Ь„1 а„— г 1 11 Таким образом, а, — з = д1, д111 = Ь1 и далее 11 1д1 1+ д« = а; — г, д111 = Ь,. Подстановка выражения 1; = Ь;/д1 в предпоследнее равенство дает простую рекурсию ь2 д1 = (а; — г) — —. (5.17) д1-1 Поскольку выбора главного элемента не производится, этот процесс можно заподозрить в опасной неустойчивости, особенно если встретится малое число д1 1.

Однако, в действительности, можно показать, что вычисления по формуле (5.17) очень устойчивы благодаря тому, что матрица А — г7 трехдиагональная [73, 74, 156). Лемма 5.3. Величины а1, вычисленные в арифметике с плаваю»цей точкой по формуле (5.17), имеют те же знаки (и, следовательно, определяют ту же инерцию), что и величины д1, точно вс«числеяные по матрице А, очень близкой к А: (А)а = а, = а; и (А),,+1 = Ь, = Ь,(1+ е;), где )е,( < 2.5е+ 0(ег).

Доказательство. Пусть д1 — величины, вычисленные с учетом ошибок окру- гления по формуле (5.17): д« = ~(а, — г)(1+ с,»г) — ' ' . (1+ е14) (1+ в,гд), (5.18) ( Ь» 1(1 + е«д) д«-1 где все «эпсилоны» ограничены по абсолютной величине машинной точностью е, а их индексы указывают, какая операпия их породила (например, е г, происходит от второгс вычитания при вычислении д1).

Определим новые 5.3. Алгоритмы дли симметричной проблемы собственных значений 243 переменные 4 (1+ е 1 ~)(1+ е дй) з у (1+ е,;)(1+ е7,) Ь, 1 = 5; г ' " = Ь; 1(1 + е;). (1 + е — Ай)(1 + е д; г)(1 + е д г — 1)1 (5,19) Полный анализ должен был бы принимать во внимание возможность переполнений и машинных нулей. В действительности, опираясь на правила 1ЕЕЕ- арифметики в отношении обработки особых ситуаций, вычисления можно вести без опаски даже в том случае, когда некоторое число г1, г является точным нулем. В самом деле, тогда имеем 4 = — оо, 4~г = аг ы — г, и далее вычисления идут как обычно [73, 8Ц. Цена одного обращения к процедуре Мейсоцпс составляет для трехдиагональной матрицы А, самое большое, 4п флопов.

Поэтому общая стоимость вычисления й собственных значений равна 0(лп). Этот подход реализован ЬАРАСК-программой ввсебя. Отметим, что метод бисекции сходится линейно: каждое деление интервала пополам прибавляет один верный разряд. Существует много способов ускорения сходимости, опирающихся на использование метода Ньютона или родственных методов для вычисления нулей характеристического многочлена (значенне которого в точке х равно произведению всех чисел г(;) [173, 174, 175, 176, 178, 2б9].

После того как вычислены (некоторые) собственные значения, вычислить соответствующие собственные векторы можно посредством обратной итерации (алгоритм 4.2); она реализована ЬАРАСК-программой ввсе1п. Поскольку в качестве сдвигов можно брать очень хорошие приближения к собственным значениям, для сходимости обычно хватает одной нли двух итераций. Так как шаг обратной итерации состоит в решении трехдиагональной системы уравнений (см. равд, 2.7.3), то для вычисления одного собственного вектора достаточно 0(п) флопов.

Если вычисленные собственные значения ггы..., гх близки друг к другу, то вычисленные для них собственные векторы ф, ..., 91 могут не быть ортогональны. В этом случае алгоритм предусматривает переортогоиализацию собственных векторов: вычисляется ЯН-разложение [д,,..., 91] = Ятг и каждый вектор дь заменяется й-м столбцом матрицы сЗ. Это гарантирует, что новые векторы дь будут ортонормированными. ЯВ;разложение обычно вычисляется посредством модифицированного процесса Грама-Шмидта (алгоритм 3.1), т. е. из каждого вектора явным образом вычитаются компоненты в направлениях ранее вычисленных векторов.

Если размер кластера и мал, то невелика и стоимость переортогонализации 0(лтв), поэтому, в принципе, все собственные значения и собственные векторы можно вычислить всего лишь за 0(пз) флопов, пользуясь вначале бисекцией, а потом обратной итерацией. Заметим, что е[, и 4 имеют один и тот же знак, а [е;[ ( 2.5е + 0(ез). Подставляя (5.19) в (5.18), получаем -"г= 1 д =а.— х — —.' г[г-г что завершает доказательство.

П Глава 5. Симметричная проблема собственных значений Это гораздо меньше, чем 0(п~) флопов, необходимых для ЯВ;итерации или (в наихудшем случае) алгоритма «разделяй-и-властвуй». Однако такое ускорение удается получить не всегда: если размер кластера велик, т. е. составляет ощутимую долю от п, то общая стоимость процесса возрастает до тех же 0(пз) флопов. Е це хуже то, что нет никакой гарантии, что вычисленные собственные векторы достаточно точны и приблизительно ортогональны. (Проблема заключается в том, что взаимное уничтожение разрядов, происходящее при переортогонализации почти линейно зависимой системы векторов дь, может привести к новым векторам, почти целиком состоящим из ошибок округления.) В преодолении этих проблем в последние годы наметился известный прогресс ~105, 83, 201, 203], и теперь кажется возможным, что обратную итерацию удастся «отремонтировать» таким образом, чтобы всегда получать достаточно точные и ортогональные приближенные собственные векторы, не тратя более чем 0(п) флопов на каждый.