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

В этом случае комбинация бисекции с «исправленной» обратной итерацией была бы наиболее предпочтительным алгоритмом во всех случаях, независимо от того, сколько собственных значений и собственных векторов требуется вычислить. Мы надеемся описать такой алгоритм в новом издании книги. Обратим внимание на удивительную внутреннюю параллельность, присущую бисекции и обратной итерации — ведь каждое собственное значение, а затем и собственный вектор можно вычислить независимо от других.

(Такое утверждение предполагает, что поправленная обратная итерация устраняет необходимость в переортогонализации больших групп приближенных собственных векторов.) Это обстоятельство делает названные алгоритмы очень привлекательными для параллельных компьютеров ~76]. 5.3.5. Метод Якоби В отличие от уже рассмотренных алгоритмов, метод Якоби не предполагает первоначального приведения к трехдиагональной форме, а работает с исходной плотной матрицей. Обычно он работает много медленнее, чем конкурирующие алгоритмы, и продолжает представлять интерес только по той причине, что иногда способен вычислять малые собственные значения и отвечающие им собственные векторы с гораздо большей точностью, чем ранее описанные методы. Кроме того, он легко распараллеливается.

Здесь мы изложим лишь основы теории метода, а обсуждение вопросов, связанных с высокой точностью, отложим до равд. 5.4.3. По заданной симметричной матрице А = Ае в методе Якоби строится последовательность ортогонально подобных матриц Аы Аз,....

В конечном счете, она сходится к диагональной матрице, на диагонали которой стоят собственные значения. Матрица Аг ы получается из А; по формуле А»»« —— ,Ут А;,У;, где У; — ортогональная матрица, называемая вращением Якоби. Таким образом, Ав« = У|» 1 Аюв-1 Ув«-1 т т = У„1У~, зА -зу -зЮ -1= ". Лт " ' Лет Ае Ло " .У -1 = гтАу 5.3. Алгоритмы для снмметрнчной проблемы собственных значений 245 сов 0 — в)п 0 дг = Ло,)с,В) = япВ сов 0 1 1 выбирается из условия, чтобы элементы (у,й) и (к, )) в А;+1 стали нулями. Чтобы найти 0 (фактнчески определяются сов В и вш О), запишем Оы) 1 т~ О) Вд а ь 1 ~ сов — япВ ) ) а,, а ь ~ ) сов — в1пВ 1 а~'ч') ~ ~ вшВ совВ ~ ~ аб) а~') ~ ~ вшВ совВ зь азу ~аз =~а — собственные значения подматрицы с Оы) ьу где Лг и Лг Вычислить сов В и яп О, для которых мы введем сокращенные обозначения с вв совВ и в в— в япВ, несложно.

Опуская, для простоты, индекс (1) в приведенном выше матричном равенстве и выполняя перемножение матриц, получаем, с учетом симметрии, ! Л1 О 1 ~ ауусз+ азьвг+ 2зса а вс(азз — а ) +а з(сг — зг) О Лг 1 '( зс(азз — а ) + ауз(сг — зг) а .вг + аьасг — 2зса ь Приравнивая внедиагональный элемент нулю, имеем а - — аьа сг — зг сов 20 — с152В = т.

2ауь 2вс яп 20 Положим 1 = -, '= сяВ и заметим, что 1г + 2т1 — 1 = О. Решая квадратное уравнение, находим 1 = ' '-, с = и в = 1 с. Суммируем этот вывод )г)+ '1+ ~' 1+гт в виде следующего алгоритма: Алгоритм 5.5. Вычислить вращение Якоби е координатной плоскости у, к и применить его к матрице А: При подходящем выборе вращений Уг матрица А для больших т будет близка к диагональной матрице Л. Следовательно, можно написать Л д Ад илн т дЛзт А. Поэтому столбцы матрицы 1 являются приближенными собственными векторами для А.

Матрица,Уг выбирается так, чтобы сделать нулями пару внедиагональных элементов матрицы Аг.ы — — дтАе,Уо Именно, угол 0 вращения у' )с 24б Глава 5. Симметричная проблема собственных значений ртос дасоог-Ко|а|гоп (Ад',й) г/ (а г~ не слишком мал т = (а„— агг)/(2.

а г г = згура(т)/(~т~ + |/1 + тг) О/0 з=с 1 А = Втаб,й,В) А Л(з,й,В) ... где с= согВ из = з|пВ г/ пуз|сны собственные векторы .тыд ЛЦ,й,В) -д/ епд г/ Стоимость применения вращения ВЦ, й, В) к матрице А (или 1) составляет лишь 0(п) флопов, так как в А изменяются только строки и столбцы с номерами |' и й (а в 1 только столбцы у и й). Алгоритм Якоби в целом можно описать Алгоритм 5.6. Метод Якоби для вычисления собственных значений с метричной матрицы: тереа| выбрать пару индексов у', й обратиться к процедуре дасоМ-Ногае|оп(Ад',й/ пока А не станет достаточно близка к диагональной матрице Остается решить, каким образом выбирать пары у',й. Имеется несколько возможностей.

Для их описания и в качестве меры продвижения к пределу определим величину оК(А) |— и ~~| агг. 1<д<й<п Таким образом, оК(А) есть норма наддиагональной части матрицы, поэтому А тогда и только тогда будет диагональной матрицей, когда оК(А) = О. Мы хотели бы, чтобы величина оК(А) уменьшалась быстро. Следующая лемма показывает, что оК(А) монотонно убывает с каждым вращением Якоби.

Лемма 5.4. Пусть з' ф й и А' обозначает матрицу, полученную в результате обращенил к процедуре дасоЬь'-Ко~пегов(Ад',й). Тогда оКз(А') = оКг(А) — аг . Доказательство. Заметим, что А и А' отличаются только строками и столбцами с номерами у и й. Имеем оКг(.4) = г г яг г ача +аз=Я +аь. |<р<г'<ь р~з или М~г Точно так же, оКг(А') = Я' + а',г — — Я', посколькУ а'дг = О в РезУльтате применения процедуры дасоЪ|-Вогасюп(А,у,й/. Используя соотношения !Щг = 1~ЯХ~~г и ~~Х~)г = 1!ХЯ()г, верные для произвольной матрицы Х и всякой ортогональной матрицы с/, можно показать, что о'~ = з' .

Следова,г тельно, оКг(А') = оКг(А) — аг„, что и требовалось. П 5.3. Алгоритмы для симметричной проблемы собственных значений 247 Опишем исходную версию алгоритма (указанную самим Якоби в 1846 г.). Она сопровождается привлекательным анализом, но, к сожалению, является слишком медленной для практического применения. Алгоритм 5.7. Классический алгоритм Якоби: юЫе оК(А) > ьо1 (здесь Со1 — параметр осгпанова, задаваемый пользователем) выбрать ) и к так, чтобы аде был наибольшим по абсолютной величине внедиагональным элементом обратиться к процедуре Хасай-Яогасгоп (А,у,й) епд юЫе Теорема 5.11. Пусть Ю = -"-(-" —:11 — число наддиагональных элементов матрицы А.

Тогда длл вращения Якоби в классическом алгоритме Якоби вьтолняетсл соотношение оК(А') < ф — фоК(А). После к вращений величина оК( ) не превышает (1 — Д) "7гоК(А). Доказательство. По лемме 5.4, оКг(А') = оК~(А) — агг. Поскольку а г — внедиагональный элемент с наибольшим модулем, то оКг(А) < -"-~-"1)агин или иначе аг„> -„-~-„-" — ~~ггоК~(А). Поэтому оКг(А) — агг < (1 — ф)оК~(А), что и требовалось. П Таким образом, если говорить о величине оК(А), то классический алгоритм Якоби сходится, по меньшей мере, со скоростью геометрической прогрессии, знаменатель которой не превосходит (1 — —. На самом деле, асимптотически метод сходится квадратично. Теорема 5.12.

Метод Якоби локально сходится с квадратичной скоростью в том смысле, что длл достаточно больших 1 справедливо соотношение оК(А;чзч) = 0(оКг(А;)). Тем самым, сравниваются состолния метода через Ж шагов (что достаточно длл того, чтобы выбрать по одному разу каэюдый элемент а и). Классический алгоритм Якоби не используется в практических вычислениях, потому что поиск наибольшего элемента слишком замедляет процесс: нужно провести поиск среди — "" элементов прежде, чем можно будет выполнить вращение, которое обойдется всего лишь в 0(п) флопов.

Это означает, что для больших ц время поиска преобладает в общем времени алгоритма. Поэтому вместо классического правила выбора пар з, к используется следующий простой метод: Алгоритм 5.8. Строчный циклический метод Якоби: построчньгй циклический обход внедиагон льнгях элементов матрицы А. гереа$ (агу=1 акоп — 1 уог к =) +1 го и обратиться к процедуре Лавой-Яо1аСгоп(Ад',Ц епд )ог епд )ог пока А не станет достаточно близка к диагональной матрице 248 Глава 5, Сиыметрнчная проблема собственных значений Если в течение всего внутреннего цикла процедура Ласоо1-Во1а«1оп дает лишь значения с = 1 и е = О, то матрица А уже более не изменяется.

Как и классический метод Якоби, циклический алгоритм асимптотически сходится квадратично [262, с. 270]. Стоимость одного «цикла» метода Якоби (в котором по одному разу выбирается каждая пара 1, Й) составляет примерно половину стоимости приведения к трехдиагональной форме и вычисления собственных значений и собственных векторов посредством ОН-итерации. Она превышает стоимость алгоритма «разделяй-и-властвуй». Поскольку для сходимости метода Якоби нередко требуется от 5 до 10 циклов, ясно, что он сильно проигрывает конкурирующим алгоритмам.

5.3.6. Сравнение производительности В этом разделе анализируется производительность трех самых быстрых алгоритмов для симметричной проблемы собственных значений: ОВ.-итерации, бисекции с обратной итерацией и алгоритма «разделяй-и-властвуй». Более подробные сведения можно найти в [10, гл. 3] или 1«ЕТЫВ/1арас)г/1пб/1арас)с)п3.

1»1ш1. Начнем обсуждение с наискорейшего метода, а потом сравним с ним другие алгоритмы. Мы пользовались ЬАРАСК-программой ввуечд. Если нужны только собственные значения, то программа реализует приведение к трех- диагональному виду, сопровождаемое ОВ-итерацией; соответствующее число операций равно «пз + 0(пз). Если же требуются и собственные значения, и собственные векторы, то вслед за приведением к трехдиагональному виду применяется алгоритм «разделяй-и-властвуй».

Измерялось время работы программы авуечб на компьютере 1ВМ ВБ6000/590, представляющем собой рабочую станцию с пиковой производительностью 266 мегафлопов. Оптимизированное матричное умножение выполняется на этой машине со скоростью 233 мегафлопа для матриц порядка 100 н 256 мегафлопов для матриц порядка 1000. В таблице приведены данные для программы ввуечб.