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

Матрица Твтв не дает никакой информации о левых сингулярных векторах матрицы В. Доказательство леммы вынесено в вопрос 5.19. Итак, к любой из трехдиагональных матриц леммы 5.5 теоретически можно применить ьгй.-итерацию, алгоритм «разделяй-и-властвуй» или бисекцию в комбинации с обратной итерацией, а затем извлечь из полученного спектрального разложения сингулярные числа и сингулярные векторы (возможно, только левые или только правые). Однако этот простой подход игнорирует особенности сингулярного разложения, что приводит к потерям в скорости и точности. Приведем два аргумента в поддержку этого утверждения.

Во-первых, применение к Тр, симметричной трехдиагональной с)В;итерации или алгоритма «разделяй-и-властвуй» было бы неэффективным, так как оба алгоритма вычисляют все собственные значения (и, если требуется, собственные векторы) матрицы Тр„тогда как, согласно лемме 5.5, нам нужны лишь неотрицательные собственные значения (и, возможно, соответствующие собственные векторы).

Имеются также трудности, связанные с точностью вычисления сингулярных векторов для малых сингулярных чисел. Во-вторых, явное формирование любой из матриц Тлвт или Твтв может привести к численной неустойчивости. В малых сингулярных числах матрицы В можно, в действительности, потерять половину верных разрядов. Положим, например, ц = е/2, тогда в арифметике с плавающей точкой число 1 + ц будет округляться до 1. Пусть В =; сингулярные числа „1, ~~.м в= ~, — 1+у округляется до точно вырожденной матрицы Тптн —— ~ .

Таким образом, округление числа 1+ ц до 1 изменяет найденное младшее сингулярное число с первоначального значения вблизи ~/ц72 = ~/е/2 до нуля. Напротив, обратно устойчивый алгоритм не должен изменять сингулярные числа более 254 Глава 5. Симметричная проблема собственных значений чем на 0(е)[[В[[а = О(е). Для 1ЕЕЕ-арифметики двойной точности е 10 'в и ~lе/2 10 в, поэтому ошибка, вызванная формированием матрицы В~В, в 10з раз превьппает ошибку единичного округления, что намного больше уровня ошибки в обратно устойчивом алгоритме. Подобная же потеря точности может происходить при явном формировании матрицы Тввт. Вследствие неустойчивости, вызываемой вычислением Тввт или Тяги, хорошие БЧП-алгоритмы работают непосредственно с В или, может быть, с Трм Резюмируя, опишем алгоритмы, используемые для практического вычисления ВЧП.

1. ОБ-итерация и ее вариантьь При правильной реализации (см. [104))— это наиболее быстрый метод отыскания всех сингулярных чисел двухдиагональной матрицы. Более того, все сингулярные числа определяются с высокой относительной точностью в смысле равд. 5.2.1. Это означает, что у всех найденных сингулярных чисел, даже у самых малых, все разряды будут верны. Сопоставим этот факт с тем, что у малых собственных значений, вычисленных трехдиагональной с4й-итерацией, верных разрядов может не быть вообще. Если нужны и сингулярные векторы, то применяется другой вариант ОЕ-итерации [80]: для вычисления сингулярных векторов, отвечающих младшим сингулярным числам, цп-итерация используется с нулевым сдвигом; в результате сингулярные векторы вычисляются так точно (а соответствующие сингулярные числа почти столь же точно), как это описано в равд.

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

Для большинства приложений такой точности достаточно. 3. Бисекция и обратная итерация. Применяя бисекцию в комбинации с обратной итерацией к матрице Т„, из первой части леммы 5.5, можно находить только сингулярные числа из заданного интервала. Этот алгоритм гарантирует высокую относительную точность сингулярных чисел, хотя, как указано в разд. 5.3.4, вычисленные сингулярные векторы могут терять ортогональность. 4.

Метод Якоби. Сингулярное разложение плотной матрицы О можно найти, применяя метод Якоби из равд. 5.3.5 к матрице Ос«т (или ОтО) неявно, т.е. так, что ни одна из названных матриц в явном виде не формируется, чем избегается потеря численной устойчивости, связанная с этим формированием. Для некоторых классов матриц С, а именно тех, к которым применима теория относительных возмущений из равд. 5.2.1, можно показать, что сингулярные числа и сингулярные векторы находятся методом Якоби с высокой относительной точностью (в смысле указанного раздела).

255 5.4. Алгоритмы вычисления сингулярного разложения В последующих разделах мы даем более подробное описание некоторых из перечисленных алгоритмов. Прежде всего, мы обсуждаем с411-итерацию и ее вариант, называемый с1цс)з (разд. 5.4.1), затем показываем высокую точность как с1цс1з, так и бисекции (разд. 5.4.2), наконец, рассматриваем метод Якоби (разд. 5.4.3). Изложение алгоритма «разделяй-и-властвуй» опущено, поскольку в целом он аналогичен алгоритму, обсуждавшемуся в разд. 5.3.3; относящиеся к нему детали можно найти в [130]. 5.4.1.

ь4й.-итерация и ее модификации для двухдиагональной ВЧР-задачи 14ВгитеРациЯ в пРименении к вычислению БЧП имеет долгУю истоРию модификаций, направленных на достижение возможно большей эффективности и точности. Хороший обзор этого вопроса дан в [200]. Алгоритм, реализованный ЬАРАСК-программой в1»йвцг, первоначально основывался на [80]. Позднее он был изменен в части, касающейся вычисления только сингулярных чисел. Для этого случая стал использоваться алгоритм из [104], называемый (по причинам, связанным с его историей) с1цс1з!. Мы начинаем с описания этого элегантного, быстрого и точного алгоритма. Выводу формул с1цс1з предпошлем изложение алгоритма, который хронологически предшествовал !3Вгитерации. Этот алгоритм, называемый ЬВг итерацией, будет описан для случая симметричных положительно определенных матриц.

Алгоритм 5.9. ЬВ;итерацият пусть То — произвольная симметричная положительно определенная матрица. Стпроится последовательность симметричных полоокительно определенных матриц Ти подобных матрице То. »=0 гереа! Выбрать сдвиг т,', меньший младшего собственного значения матрицы Т;. Вычислить разложение Холесского Т! — т21 = ВтВ; (Вс — верхняя треугольная матрица с положительной диагональю.) Т;< ! — — ВсВ2 + тз1 «=«+1 пока не будетп достигнута сходимость По своей структуре, ЬВгитерация весьма схожа с !4Вгитерацией: сначала вычисляется разложение текущей матрицы, а затем сомножители перемножаются в обратном порядке, что дает следуюшую матрицу Тс«!. Легко видеть, что Тс ь! и Т подобны.

Тс«.! — — В.В2 +т21 = В ВтВсВт+тгВ Вт — В Т»Вт В действительности, можно показать, что при тз = 0 два шага ЬВгитерации порождают ту же самую матрицу Тз, что и один шаг !.4Вгитерации. г с10ав — это сокращение от «гййегепт!а1 Чпог!епг-ссгяегепсе а1йогнйгп н!СЬ вмяв» (209) (дифференциальный алгоритм частных и разностей со сдвигами). 256 Глава 5. Симметричная проблема собственных значений Лемма 5.6.

Пусть Тг — матрица, порождаемая двумя шагами алгоршпма 5.9 при гг = О, а Т' — матрица, получаемая из То одним шагом ЯЯ- итерации (т. е. ьзЯ = То, Т' = Ж„)). Тогда Тг — — Т'. Доказательство. Используя симметрию матрицы То, разложим Тот в произведение следующими двумя способами. Во-первых, Тг = То То — — (сзЯ)г сзЯ = ЯтЯ.

Без ограничения общности, можно считать, что все элементы Яи положительны. Итак, То представлена как произведение нижнетреугольной матрицы Я~ и (транспонированной) матрицы Я; это представление должно совпадать с представлением Холесского в силу единственности последнего. Матрицу То можно факторизовать иначе, а именно То = Во ВоВо Во. Согласно ало т т горитму 5.9, Т1 = ВоВо т= ВтВм поэтому Тог ВотВоВотВо = Вот(В1тВ1)Во = (В1Во)т(В~Во). Снова получено представление Тг произведением нижнетреугольной матрицы (В1Во) и транспонированной матрицы; следовательно, оно также должно совпадать с разложением Холесского.