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

В этом разделе будет показано, как с помощью метода Якоби найти сингулярное разложение плотной матрицы С; именно, алгоритм 5.8 из равд. 5.3.5 нужно неявно применить к симметричной матрице А = Ст'С. Отсюда следует, что сходимость этого метода вычисления 8УР примерно такая же, как алгоритма 5.8; в частности, и для задачи вычисления БУР метод Якоби является самым медленным из имеющихся. Тем не менее, интерес к методу Якоби сохраняется, потому что для некоторых типов матриц С он способен вычислять сингулярные числа и сингуляр- 262 Глава 5.

Снмметрнчнгя проблема собственных значений ные векторы намного точнее, чем другие обсуждавшиеся нами методы. Для этих матриц С метод Якоби обеспечивает высокую относительную точность (в смысле ргзд. 5.2.1) вычисленных приближений к сингулярным числам и сингулярным векторам. После описания неявного метода Якоби мы покажем, что он вычисляет БУР матрицы С с высокой относительной точностью, если С может быть представлена в виде С .= РХ, где Р— диагональная матрица, а Х хорошо обусловлена. (Наличие такого представления означает, что С плохо обусловлена в том и только в том случае, когда матрица Р имеет и большие, и малые диагональные элементы.) Более общо, мы получаем выигрыш всякий раз, когда Х обусловлена значительно лучше, чем С.

Мы проиллюстрируем это положение с помощью матрицы, для которой всякий алгоритм, включающий в себя приведение к двухдиагональной форме, необходимо теряет все верные разряды во всех сингулярных числах, кроме старшего; в то же время, метод Якоби определяет все сингулярные числа этой матрицы с полной машинной точностью. Затем мы дадим обзор других классов матриц С с тем же свойством, т.е. метод Якоби для них значительно точнее, чем методы, опирающиеся на двухдиагонализацию. Заметим, что в разд.

5.4.2 было показано: БУР двухдиагональной матрицы С можно вычислить с высокой относительной точностью, применяя бисекцию или алгоритм дск)г (см. разд. 5.4.1). Проблема состоит в том, что преобразование плотной матрицы в двухдиагонгльную форму может сопровождаться ошибками, достаточными для потери высокой относительной точности; именно это покажет наш пример. Поскольку в методе Якоби заданная матрица обрабатывается без предварительного приведения к двухдиагональному виду, метод имеет гораздо больше шансов достичь высокой относительной точности.

Неявный метод Якоби математически эквивалентен применению алгоритма 5.8 к матрице А = СтС. Иначе говоря, на каждом шаге вычисляется вращение Якоби д, с помощью которого матрица СтС неявно пересчитывается в .1тСтСХ; вращение выбрано так, чтобы пара внедиагональных элементов из СтС обратилась в нули в матрице дтСтСд. Однако ни СтС, ни .7тСтСд не вычисляются в явном виде; вместо них вычисляется матрица Сд. По этой причине мы называем алгоритм методом односторонних вращений. Алгоритм 5.12.

Вычислить вращение Якоби в координатной плоскости т', й и применить его (односторонне) к матрице С: ргос Опе-Вгдед-ЮасоЬь'-Во1айоп (С, ), й) В,„и и„, а (СтС) а (СтС) „и а (СтС)„„ г/~аде~ не слишком мал т = (а,, — агг)/(2 а г) 1 = г1яп(т)/(/т/ + Я + тг) с = 1/,/1 + гг г=с-$ С = С В(т',й,д) ...где с = созд и г = гшд г/ нулсны правые сингулярные векторы д = д Л(), й, В) епй г/ епд г/ 5.4.

Алгоритмы вычисления сингулярного разложения Заметим, что в данной процедуре вначале определяются элементы (~,1), (,1, к) и (й, к) матрицы А, после чего вращение Якоби вычисляется таким же образом, как в процедуре ЛасоЪ1-КоФа11оп (алгоритм 5.5). Алгоритм 5г13. Метод односторонних вращений: пусть С вЂ” матрица размера и х п. Алгоритм вычисляет сингуллрнтяе числа и;, матрицу У левых сингуллрных векторов и матрицу И правых сингулярных векторов. Если положить Е = Йад(от), то С = УЕ\гз'. тереаФ ~от1=1 Фоп — 1 (от й = т + 1 1о п обратиться к процедуре Опе-Етдед-Лавой~'-йо1а1топ(С,1', й) епд )от епд )от пока СтС не стпанета достаточно близка к диагональной матрице Положить о, = ((С(:, т)Йг (2-норма т'-го столбца в С) Положить У = (иы..., и„], где и, = С(: т)/от Положить И = У, где .У вЂ” накопленное произведение вращений Якоби В вопросе 5.22 предлагается доказать, что матрицы Е, У и )т, вычисляемые методом односторонних вращений, действительно дают 51тР матрицы С.

Мы покажем сейчас, что, несмотря на округления, метод односторонних вращений способен вычислять ЗЧР с высокой относительной точностью, если матрица С может быть представлена в виде С = РХ, где Р— диагональная матрица, а Х хорошо обусловлена. Теорема 5.15. Пусть С = РХ вЂ” матрица порядка и, причем Р и Х невырожденны и Р— диагональная матрица. Пусть С вЂ” матрица, полученнал из С в результате т-кратного обращения к процедуре Опе-ЕтдедУасоЫ-Яо1а1топ(С,1', и) в арифметике с плавающей тпочкой. Обозначи.м через от » о„и о, » ° о„сингулярные числа матриц С и С. Тогда ( 0(тв)к(Х), (5.25) ттт где к(Х) = ()ХЙ 'ИХ т)~ естпь число обугловленностпи матприцы Х. Другими словами, относительная погрешность вычисленных сингулярных чисел мала, если мало число обусловленности матрицы Х.

Доказательство. Рассмотрим вначале случай т = 1, т. е. результат применения одного вращения Якоби, а затем обобщим наш анализ на случай произвольного т. Изучение процедуры Опе-ЗЫет1-ЛасоЪ|-гсосас1оп(С, т', тт) показывает, что С = й(С гт), где  — вращение Гивенса, вычисленное в арифметике с плавающей точкой. Но построению, А отличается от некоторого точного вращения В матрицей с нормой 0(е). (Неважно (и даже не обязательно верно), чтобы Л отличалась на 0(е) от . подлинногот вращения Якоби, которое процедура ОпеЗвдет1-дасоЪ|-Нотастоп(С, т', й) вычислила бы в точной арифметике.

Необходимо лишь, чтобы Л отличалась на О(е) от некоторого вращения. Это требование сводится к соотношению сг + зг = 1 + 0(е), которое легко проверить.) 264 Глава 5. Симметричная проблема собственных значений Ниже мы покажем, что д = ОВ(1+ Е) для некоторой матрицы Е с малой нормой: ЦЕЦг = 0(е)к(Х). Если бы Е была нулевой матрицей, то д и СВ имели бы одни и те же сингулярные числа, поскольку А — (точно) ортогональная матрица. Если норма (ненулевой) матрицы Е меньше единицы, то относительные разности сингулярных чисел можно оценить с помощью следствия 5.2: ( Ц(1'+ Е)(1 + Е)~ 1Ц ЦЕ+ Е~ +~ЕЕ~ ( 3ЦЕЦ пг = 0(е)к(Х), (5.26) что и требовалось.

Построим требуемую матрицу Е. Поскольку Л является для О правым множителем, каждая строка в д зависит лишь от соответствующей строки в О; в обозначениях Маг!аЬ'а, это можно записать так: д(г,:) = Я(О(г,:) Л). Положим Е = д — ОЕ. Используя лемму 3.1 и соотношение С = РХ, имеем ЦЕ(г,:)Цг = Цд(г,:) — О(«,:)ВЦг = 0(е)ЦО(«,:)Цг = 0(е)ЦЕХ(г,:)Цг, откуда Цб,,~Е(г,:)Цг —— 0(е)ЦХ(«,:)Цг или ЦР ~ЕЦг = 0(е)ЦХЦг. Учитывая соотношения Е ' = т«т и О ' = (РХ) г = Х 'Р г, получаем д = ОЕ+ Е = ОВ(1 + ЕтО 'г) = ОЕ(1 + Г«~Х 'Р 'г ) ив з ОЕ(1+ Е) где ЦЕЦг ( ЦВ ЦгЦХ ЦгЦР гЕЦг = 0(е)ЦХЦгЦХ Цг = 0(е)к(Х).

Именно это и требовалось. Остается распространить этот результат на нг вращений, где т > 1. Заметим, что в точной арифметике мы имели бы д = ОЕ = РХЕ = РХ, где к(Х) = к(Х). В этом случае к каждому из т шагов была бы применима оценка (5.26), что влекло бы за собой (5.25). Вследствие округлений, к(Х) может приобретать на каждом шаге множитель, не превосходящий к(1+ Е) ( 1 + 0(е)к(Х). Эти множители, очень близкие к 1, учитываются символом 0(пге) в оценке (5.25).

П Чтобы завершить описание алгоритма, нужно указать критерий астапова, т.е. прояснить смысл условия «если ~агн~ не слишком малг в алгоритме 5.12 (т. е. в процедуре Опе-ВЫес1-ЛасоЬЬКогаг1оп). Подходящим критерием является условие )арг~ > е~/атгаеь.

Дальнейшее его обсуждение дается в вопросе 5.24. Пример 5.9. Рассмотрим, в известном смысле, экстремальный пример матрицы 0 = РХ, все сингулярные числа которой вычисляются методом Якоби с полной машинной точностью. В то же время, любой метод, опирающийся на приведение к двухдиагональному виду, вычисляет с полной машинной точностью лишь наибольшее сингулярное число АЗ; во всех остальных приближенных сингулярных числах верных разрядов нет вовсе (хотя их абсолютные погрешности суть величины вида хО(е) ~/3, как и следует ожидать от обратно устойчивого алгоритма).

Для нашего примера е = 2 ег 10 ге (двойная 265 5.4. Алгоритмы вычисления сингулярного разложения точность 1ЕЕЕ-арифметики), а тт = 10 зо (в действительности, годится почти любое значение тт < е). Положим =В.Х С точностью до, по меньшей мере, 16 разрядов, сингулярными числами матрицы С являются 3, з/3. тт, тт и тт. Чтобы убедиться в том, что при приведении С к двухдиагональной форме происходит потеря точности, достаточно рассмотреть первый шаг алгоритма из равд. 4.4.7. Этот шаг состоит в левом умножении на отражение с целью аннулировать элементы 2, 3 и 4 столбца 1. В точной арифметике С имела бы вид Однако, вследствие малости числа тт, округления приводят к матрице — 2тт —.5 —.5 —.5 0 —.5 —.5 —.5 0 —.5 —.5 —.5 0 —.5 —.5 —.5 Обратим внимание на то, что в трех последних столбцах матрицы Ст впотеряна» вся информация об тт.

Поскольку эти столбцы совпадают, Ст (точно) вырождена; в действительности, ее ранг равен 2. Таким образом, два младших сингулярных числа изменили свои значения с тт на О, т. е., в относительном смысле, точность в них полностью потеряна. При отсутствии новых округлений мы привели бы Ст к двухдиагональной матрице — 2п К755 1.5 0 0 0 0 имеющей сингулярные числа з/3, з/3 и, 0 и 0; первые два из них суть точные сингулярные числа матрицы С. Так как, однако, при приведении Ст к двух- диагональной форме округления происходят, в нулевые элементы матрицы В будут внесены ненулевые возмущения порядка 0(е); в результате неточными станут все три младших сингулярных числа. Два наименьших ненулевых сингулярных числа, найденных алгоритмом, суть величины порядка е и полностью состоят из ошибок округлений.