Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 54
Текст из файла (страница 54)
по этому поводу равд. 1.5 и вопрос 1.18. К сожалению, эти эффективные алгоритмы моделирования не пригодны для компьютеров Сгау 2, УМР и С90, где округления производятся недостаточно точно. В конце концов, была найдена другая формула, устраняющая необходимость в моделировании арифметики повышенной разрядности. Она основана на следующей теореме Левнера [129, 179]: Теорема 5.10 (Левнер). Пусть для диагональной матриц»1 Р=Жак(д1,..., д„) справедливы неравенства д„«. ° д1, и пусть заданы числа оп « п1, удовлетворяющие условиям перемежаемости дь < сгь « ' ' ' ч1+1 < сг««-1 < дг < 111 « ' ' ч1 < с»1 ° Тогда найдется вектор й, такой, что сг, суть точные собственнь1е значения матрицы Р = Р + йй .
Компоненты вектора й выражаются формулами 5.3. Алгоритмы для симметричной проблемы собственных значений 237 Доказательство. Характеристический многочлен матрицы Р можно запи- сать, с одной стороны, как 11еь(Р— М) = П,(сг — Л), а с другой стороны, как пес(Р— Лу) = П(ау — Л) 1=1 ш у и П и2 = П(д -') Е„', 2=1 у — 1 1 УФ1 П (дб — Л) 2=1 УФ1 . й2 Здесь использованы формулы (5.13) и (5.14). Полагая Л = дэ и приравнивая оба выражения для 11ес(Р— ЛГ~, имеем П( ° -') =< П(д,-дэ) 1=1 1=1 утн Сформулируем теперь устойчивый алгоритм вычисления собственных значений и собственных векторов (чтобы упростить изложение, положим р = 1).
Алгоритм 3.3. Вычисление собственных значений и собственных векторов матрицы .Р + ии~ . Решая вековое уравнение 1+ 2,'," "'„= О, найти (приближенные) собственные значения сг, матрицы Р + ииг. Используя теорему Левнеро, вычислить вектор й, такой, что ои являются «точными» собсп1венными значениями матрицы Р+ иит. С помощью леммьг 5.3 вьэчислить собственньге векторьэ матрицы Р +им~. Приведем краткое обоснование устойчивости этого алгоритма. Анализируя критерий астапова в программе, решающей вековое уравнение, можно показать, что ~уии — ййт ~Я2 < 0(е)ЯР32 + Оиит'52)'. Это означает, что Р + ии~ и Р+ ййт настолько близки, что собственные значения и собственные векторы второй матрицы суть устойчивые приближения собственных значений и собственных векторов первой.
Заметим теперь, что в формулу для й; в теореме Левнера входят только разности Н вЂ” с(, и сг — дэ чисел с плавающей точкой, произведения и частные этих разностей и квадратный корень. Предположим, что л«ашинн я арифметика с плавающей «почкой достаточно точна в том смысле, что Я(а О Ь) = (а О Ь)(1+ 6) для всех операций О Е (+,—, х,7') и 1 Говоря более подробно, для справедливости этой оценки необходимо, чтобы программа определяла не о„а наименьшую из разностей о, — д, и д,»1 — о,. П,= ( — д1) П,=1, у~*(113 — д1) Используя свойство перемежаемости, можно показать, что дробь в правой части положительна, поэтому извлечение квадратного корня возможно и дает требуемое выражение для й,.
П 238 Глава 5. Симметричная проблема собственных значений вс)гс(а) = з/а (1+ б), где )б! = О(г). Тогда вычисления по формуле Левнера можно провести с высокой относительной точностью. В частности, легко показать, что, в отсутствие переполнений и машинных нулей, где ~б~ = 0(г). С высокой относительной точностью можно реализовать и формулу из леммы 5.2, поэтому собственные векторы матрицы Р + ййт могут быть вычислены с высокой относительной точностью. В частности, с большой точностью выполняются соотношения ортогональности.
Резюмируем сказанное таким образом: если машинная арифметика с плавающей точкой достаточно точна, то алгоритм 5.3 с очень высокой точностью определяет собственные значения и собственные векторы матрицы Т) + ййт, слабо отличающейся от исходной матрицы В + им~. Это свойство означает численную устойчивость алгоритма.
Читатель должен отметить, что именно нарушение этого требования достаточной точности машинной арифметики привело к тому, что на некоторых компьютерах Сгау не работает моделирование учетверенной точности, предложенное в [234, 204]. Пока нам не удалось разработать алгоритма, который бы надежно работал на этих компьютерах.
Необходим еще один трюк. Дело в том, что недостаточно точными операциями этих компьютеров являются только сложение и вычитание, причем это связано с отсутствием так называемого вспомогательного разряда в аппаратной реализации арифметики с плавающей точкой. Это означает, что при сложении или вычитании самый правый бит операнда может рассматриваться как О, даже если, в действительности, он равен 1. Бсли болыпинство старших разрядов взаимно уничтожаются, то этот «потерянный битв становится значимым. Например, при вычитании 1 из ближайшего к ней, но меньшего числа с плавающей точкой (это вычитание приводит к взаимному уничтожению всех разрядов, кроме последнего) на компьютере Сгау С90 получается результат, вдвое больший правильного, а на Сгау 2 получается О. Однако если самый правый бит действительно равен О, то никакой ошибки не происходит.
Это подсказывает такой трюк: перед применением теоремы Левнера или леммы 5.2 в алгоритме 5.3 намеренно заменить самые правые биты во всех числах д«на О. Алгоритм при этом сохраняет устойчивость, так как числа д» и ол претерпели лишь небольшие относительные изменения . 1 Более подробное описание данного алгоритма можно найти в (129, 1311. ЬАРАСК-программа и1ае63 является его реализацией. г Чтобы заменить самый правый бит числа с плавающей точкой б на О на машине Сгау, нь обходимо лишь выполнить присэаиэание б:= (б+ б) — б. Это несложное вычисление вовсе не меняет б на машине с «праэильной» двоичной арифметикой (исключая переполнения, которых легко избежать). По на машине Сгау оно действительно заменяет самый правый бит нулем.
Приглашаем читателя, знакомого с арифметикой этих машин, доказать данное утверждение. Едиистэеиная остающаяся трудность — помешать оптимизирующему транслятору удалить это присэаиэание из программы; некоторые чересчур усердные оптимизаторы, действительно, могут это сделать. Для нынешнего поколения транслятороэ указанную трудность можно преодолеть, эычнсляя б + й посредстэом обращения к функции, хранящейся э файле, отличном от основной программы.
Мы надеемся, что к тому времени, когда трансляторы будут усовершенствованы до такой степени, что смогут оптимизировать и эту ситуацию, арифметика машин Сгау отомрет окончательно. 5.3. Алгоритмы для симметричной проблемы собственных значений 239 Ускорение алгоритма «разделяй-и-властвуй» посредством метода РММ Метод РММ был первоначально предложен в [124] для совершенно другой задачи вычисления электрических сил взаимодействия в ансамбле из п заряженных частиц или гравитационных сил в системе из и масс.
Мы ограничимся указанием связи этих задач с вычислением собственных значений и собственных векторов, адресуя читателя за деталями к [131]. Пусть п1,..., и'„— трехмерные радиус-векторы и частиц с зарядами з, иь а а1,..., а„— радиус-векторы других п частиц с единичными положительными зарядами. Согласно закону обратных квадратов, сила, действующая на частицу в позиции а, со стороны частиц в позициях с11, ..., и'„, пропорциональна величине з1и1[111 — а.) — [[;-а1[[з Если электростатика моделируется в двумерном пространстве вместо трехмерного, то силы подчиняются закону обратной пропорциональности' з,и,[111 — а,) [[111 — 1[Я Поскольку п1 и а; — это векторы в Ра~, их можно интерпретировать и как комплексные числа. В таком случае п ь=Е„'"", 1=1 1 1 где п1 и а — числа, сопряженные соответственно с п1 и а .
Если все п1 и ау окажутся вещественными числами, то произойдет еще большее упрощение: Л = ~,"".1 Рассмотрим теперь матрично-векторное умножение ут = хгч', где Я' есть матрица собственных векторов для Р + иит. Согласно лемме 5.2, ф, = и1вт/[111 — а ), где масштабирующий множитель в выбран так, чтобы учй столбец был нормирован. Поэтому у-й элемент вектора ут = з~Я' равен А';, = ау~У„„* ' 1=1 1=1 т.е. [с точностью до множителя в ) той же сумме, что и электростатическая сила. Итак, самая трудоемкая часть алгоритма «разделяй-и-властвуй», а именно матричное умножение в последней строке алгоритма 5.2, эквивалентна вычислению электростатических сил.
Вычисление этой суммы для у = 1,..., и, по видимости, требует 0[па) флопов. С помощью РММ и других родственных методов [124, 23] ее можно вычислить [приближенно, но с очень хорошей точностью) за время 0(п 1окп) Технически вто означает, что потенциал удовлетворяет уравнению Пуассона в двумерном пространстве, а не трехмерном. 240 Глава 5. Симметричная проблема собственных значений или даже 0(п).
(См. лекции «Быстрые иерархические методы в задаче А' тел» на РАНА1Д Е1 НОМЕРАСЕ. В них даны все необходимые подробности.) Однако этой идеи самой по себе недостаточно для того, чтобы снизить цену алгоритма «разделяй-и-властвуй» до 0(п 1ой" п). В конце концов, матрица Я на выходе алгоритма состоит из и элементов, что, по-видимому, должно означать, что цена алгоритма никак не может быть меньше, чем пз. Поэтому мы должны научиться представлять я, используя менее чем пз независимых чисел. Это оказывается возможным, так как трехдиагональная п х п-матрица имеет лишь 2и — 1 «степеней свободы» (например, диагональные и наддиагональные элементы), в качестве которых можно принять п собственных значений, оставив еще п — 1 для ортогональной матрицы Я. Иначе говоря, не всякая ортогональная матрица может быть матрицей собственных векторов для симметричной трехдиагональной матрицы Т.
На эту роль годится лишь (и — 1)-мерное подмногообразие всего многообразия ортогональньгх матриц, размерность которого равна п(п — 1)/2. Матрица Я представляется с помощью дерева, вычисляемого алгоритмом 5.2. Вместо того чтобы производить накопление по формуле [ ].Я, мы сохраняем все матрицы (~, по одной на каждый узелдереа 0 ! ва. Кроме того, вместо хранения 1,1' в явном виде, мы храним соответствующие величины Р, р, и и собственные значения си матрицы Р + рви~. Мы можем действовать таким образом, потому что в методе РММ требуется лишь уметь умножать Я' на векторы.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
