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

Таким образом, единственность разложения Холесского позволяет нам заключить, что Я = ВгВо. Тем самым найдена связь между двумя шагами ВК-итерации и одним шагом ЯН-итерации. Мы используем ее следующим образом: из То — — СзЯ следует, что Т' = ЯЯ = Ш~(ЯЯ ') = ЯЯЯ)Я ' = ЯТоЯ ' поскольку То — — (~Я = (В1Во)(Во Во)(В~Во) " поскольку Я = В1Во и То = ВтВо : В1ВоВо ВоВо Вг = Вг(ВоВо )В1 = В1(В1 В1)В 1 поскольку ВоВо = Т1 — — В~В1 = В1ВТ = Тг, что и требовалось. П Ни алгоритм 5.9, ни лемма 5.6 не требуют, чтобы матрица То была трехдиагональной; достаточно, чтобы она была симметричной и положительно определенной. Используя связь между 1 В- и ьЗВ;итерациями, установленную в лемме 5.6, можно показать, что значительная часть анализа сходимости, проведенного для с)В;итерации, переносится на ЬВ-итерацию; однако здесь мы не будем этого делать.

Алгоритм с)цдз, являющийся нашей целью, математически эквивалентен ВН-итерации. Однако он реализуется не так, как описано в алгоритме 5.9, по- сколькУ последний пРедполагает фоРмиРование матРицы Тою — — В,В, + тг1 в явном вцде. В равд.

5.4 было показано, что такая операция может быть численно неустойчивой. Поэтому матрица В;+1 будет вычисляться непосредственно из В,, а промежуточная матрица Т;+1 не будет формироваться вовсе. Чтобы упростить обозначения, положим, что В; имеет диагональные элементы аз, ..., а„и надциагональные элементы Ьы ..., Ь„м а Ваьг — диагональные элементы ам ..., а„и налдиагонзльные элементы Ьы..., Ь„г, Условимся, что Ьо = Ьо = Ь„= Ь„= О. Запишем соотношение, связывающее В, и В,+1.

В~~~Визг + та+11 = Тььг = В,ВТ+ тг1. (5.20) Приравнивая элементы (~,1), г' ( п, в левой и правой частях (5.20), имеем 5.4. Алгоритмы вычисления сингулярного разложения 257 где б = т~ь, — тэ. Будем считать, что последовательность чисел 7.э сходится к младшему собственному значению матрицы Т, монотонно возрастая (тогда все матрицы Т; будут положительно определены и алгоритм является корректным); следовательно, б ) О. Приравнивая теперь квадраты элементов (у, у'+ 1) в левой и правой частях (5.20), находим (5.22) Объединяя равенства (5.21) и (5.22), получаем следующий пока еще не окончательный алгоритм: Гогу = 1 Фо и — 1 аэ = аз+ Ьй — Ьэ — б 1 1 52 52, (а2 /52) уег епг1 Гог "2 э 52 Эта весьма экономичная версия алгоритма содержит в своем внутреннем цикле всего пять операций с плавающей точкой. Квадраты элементов матрицы В, отображаются непосредственно в квадраты элементов матрицы В;~г.

Извлекать квадратные корни не имеет смысла до самого конца алгоритма. В самом деле, на современных компьютерах извлечение квадратного корня, как и деление, может потребовать в 10-30 раз больше времени, чем сложение, вычитание или умножение, поэтому этих дорогостоящих операций нужно по возможности избегать. Чтобы подчеркнуть, что вычисляются именно квадраты элементов, заменим переменные по формулам 91 = а~~ и е, = Ьэ. Это приводит к нашему предпоследнему алгоритму цпэ (название снова связано с историей алгоритма, на которую мы не будем отвлекаться ~209]). Алгоритм 5.10.

Шаг алгоритма дде: /ог 1 = 1 го и — 1 ф=дд+ед — е, г — б е, = ед (щ.г/91) епд /ог д„ = д„ вЂ” е„ г — б Наш окончательный алгоритм дсГдэ выполняет примерно такое же количество работы, что и цаэ, однако значительно более точен, что будет показано в разд. 5.4.2. Выделим из первой строки алгоритма 5.10 выражение д — е, г — б и перепишем его следующим образом: д=д — е г — б = д — ~„~ — б согласно (5.22) 91 — г 91 г — еу г о, — б согласно (5.21) 91-г 1 — ° ду г — б ГВ 91-г 258 Глава 5. Симметричная проблема собственных значений Это позволяет нам записать внутренний цикл алгоритма 5.10 так: ф =д.+еу ед — — ед (яд+с(с)1) д,+ =дэ (уд Нд) — 5 Заметим, наконец, что с(,.вс может быть записано на место д, и отношение 1 = ддз.сЩ достаточно вычислить один раз.

Это приводит нас к окончательному алгоритму с1с1с1в. Алгоритм 5.11. Шаг алгоритма дс1дэс д= дс — б (ог1'=1 гоп — 1 ф = с(+е г = (одэ., /с),) е,=е, 1 д=д $ — б епд 1ог дп = а Во внутреннем цикле алгоритма с1цс1в то же самое число операций с плавающей точкой, что и в алгоритме с1с)в, только одно вычитание заменено умножением. В следующем разделе мы увидим, что это изменение полностью окупает себя, гарантируя высокую относительную точность результатов.

Мы не коснулись двух важных вопросов: выбора сдвига Ю = тг+, — тг и распознавания сходимости. Их подробное обсуждение можно найти в [104]. 5.4.2. Вычисление ЯЧП двухдиагональной матрицы с высокой относительной точностью Этот раздел, опирающийся на содержание равд. 5.2.1, может быть опущен при первом чтении книги. Возможность вычисления БЧР двухдиагонвльной матрицы В с высокой относительной точностью (в смысле равд. 5.2.1) основывается на приводимой ниже теореме 5.13. Согласно этой теореме, малые относительные возмущения элементов матрицы В влекут за собой лишь малые относительные возмущения ее сингулярных чисел. Лемма 5.7. Пусть  — двухдиагональнол матрица с диагональными элеменспами аы ..., а„и наддиагональными элементами бы ..., Ь„ы Пусть дана есце одна деухдиагональнал матрица В с диагональными элементами а; = акис и наддиагональными элементами ус = б,с,с.

Тогда В = РсВРг, где Р 1 ( ХгХс ХвХгХг Хь ''Хс 66 ' 'ч -с" 6эс ' Рг =с11а8 1,—,—,..., ч'с ~гй 1п-с чг "с Хс ХгХс Х -с'''Хс) Лемма проверяется несложным вычислением (см. вопрос 5.20). Применяя теперь следствие 5.2, приходим к следующему утверждению.

5.4. Алгоритмы вычисления сингулярного разложения 259 Теорема 5.13. Пусть В и  — матрицы, определенные в лемме Б.7. Предположим, что найдется числот > 1, такое, что т ' < зч < т и т ' < ~; < т. Иначе говорл, число е г— и т — 1 есть граница длл оглносительного различия произвольного элемента матрицы В и соответствующего элемента матрицы В.

Пусть оп « ог — сингулярные числа матрицы В, а д « дг — сингулярные числа матрицы В. Тогда ~о, — о;~ < и;(т4п г — 1). Если и, ф 0 и т — 1 = е «1, то можно написать )дг — о;) <,4п г 1 (4п 2) +О(г) ог Таким образом, относительные возмущения сингулярных чисел, т.е. величины ~ог — о,~/он ограничены произведением относительного возмущения е матричных элементов и числа 4п — 2. При небольших дополнительных усилиях множитель 4п — 2 может быть уменьшен до 2п — 1 (см.

вопрос 5.21). Можно показать, что и сингулярные векторы определяются весьма точно; погрешность в них будет обратно пропорциональна относительной отделенности (характеристике, введенной в рззд. 5.2.1). Покажем, что сингулярные числа двухдиагональной матрицы можно найти с высокой относительной точностью, пользуясь как бисекцией (алгоритм 5.4 в применении к матрице Тр, из леммы 5.5), так и методом дцдз (алгоритм 5.11). Начнем с бисекцин.

Вспомним, что собственными значениями симметричной трехдиагонзльной матрицы Тг, являются сингулярные числа матрицы В и числа, им противоположные. Из леммы 5.3 следует, что инерция матрицы Тр, — Л1, вычисленная с помощью формулы (5.17), есть точная инерция матрицы того же вида, построенной по некоторой матрице В, такой, что относительные разности одноименных элементов в В и В не превосходят 2.5е.

Поэтому, согласно теореме 5.13, относительные разности точных н вычисленных сингулярных чисел (т.е. сингулярных чисел матрицы В) не превысят величины, примерно равной (10п — 5)е. Обратимся теперь к алгоритму 5.11. Пользуясь теоремой 5.13, покажем, что сингулярные числа матрицы В (т. е. матрицы на входе алгоритма 5.11) и сингулярные числа матрицы В (матрицы на выходе этого алгоритма) согласуются с высокой относительной точностью.

Отсюда будет следовать, что, когда после многих шагов алгоритма дс1дз матрица В станет почти диагональной и за ее (приближенные) сингулярные числа можно будет принять диагональные элементы, эти элементы будут приближениями высокой относительной точности к сингулярным числам первоначальной матрицы на входе алгоритма. Проще всего разобраться в ситуации с нулевым сдвигом д. В этом случае в ЙуЬ выполняются только сложения положительных чисел, умножения н деления; взаимного уничтожения верных разрядов нигде не происходит. Упрощая, можно сказать, что всякая последовательность выражений, построенных из указанных основньгх операций, гарантированно дает высокую относительную точность в каждом вычисленном результате.

Итак, В будет вычислена с высокой относительной точностью, а тогда, согласно теореме 5.13, с высокой относительной точностью будут согласованы и сингулярные числа матриц В и В. Общий случай, когда Б > О, труднее для анализа (104]. 260 Глава 5. Симметричная проблема собственных значений Теорема 5.14. В арифметике с плавающей точкой шаг алгоритма 5.11, при- менлемого к матрице В и вычисляющего матрицу В, эквивалентен следую- щей последовательности операций: 1. Произвести малые относительные возмущения (не превышающие по абсолюптой величине 1.5е) в каждом элементе матрицы В. Пусть В— полученная матрица.

2. Применить к В шаг алгоритма 5.11 в точной арифметике. Пусть В— полученная матрица. 3. Произвести малые относительные возмущения (не превышающие по абсолютной величине е) в каждом элементе матрицы В, получив в результате матрицу В. По теореме 5.13, шаги 1 и Я влекут за собой лишь малые относительные изменения сингулярных чисел двухдиагональной матрицы. Поэтому сингулярные числа матриц В и В должны быть согласованы с высокой относительной точностью. уд = (Нд + е ) (1 + е; «) 11 = (Яз'.«г(уд)(1+ ед 1) е, = ез 11(1+ едлг) дз«г = (Нд 11(1+ ед.г) — б)(1+ ед,— ) Подставляя выражение для уд из первой строки во вторую, находим уд+1 1 + ед / Н +е, 1+е1+ Подставляя это выражение для 8 в последнюю строку алгоритма и деля все члены на 1+ е, имеем Й1Я1«г (1 + е11)(1 + ед «г) б.

(5.23) 1+с, Нд+е 1+с е Отсюда видно, как следует определить матрицу В. Положим дз-«г дз.«г = 1+ ед е. е+г = 1+ед г (1 + ед /)(1 + ед *г) Уд+г = чз.«г 1 + ел« (5.24) Тогда (5.23) принимает вид г. дд а+1 б 1+1— ф +е. Доказательство. Перепишем внутренний цикл алгоритма 5.11, снабдив пере- менные Н и 1 индексами с тем, чтобы можно было различать их значения на разных итерациях, а также учтя округления посредством членов вида (1+ е), тоже снабженных индексами: 5А, Алгоритмы вычисления сингулярного разложения Формулы (5.24) показывают, что наибольшее относительное изменение в элементах В по сравнению с В не превышает 1.5е (поскольку в выражении для дт+г — — 53 „,,+, присутствуют три множителя вида 1+ е). Определим теперь величины 51 и е, для матрицы В формулами Дз=4 +е Ь = Иу+1/чз) е =ез ° 1~ 41+, =41 13 — 5.

Они описывают шаг алгоритма бцбэ, применяемого к матрице В в точной арифметике; результатом является матрица В. Остается показать, что В отличается от В относительным изменением каждого элемента, не превосходящим по абсолютной величине е. Это следует из представлений 51 = с~а+ в. + 1+ее 1 1+аз г 1 = (с~у + е ) (1 + е гь ) . 1 Ц 1 (1+ етгь)(1 + е1 э ) еу=еэ 11 е вв1+ 1+аз г 41 е ' 11(1 + ез ~2)(1 + ез ц — ) 1 + еу 1+9, г = е гз(1+ еэ.г) 1+аунг 1+ ет,,а 5.4.3. Метод Якоби дли вычислении Б 1г0 При обсуждении метода Якоби в равд. 5.3.5 было отмечено, что этот метод является самым медленным из имеющихся алгоритмов вычисления собственных значений и собственных векторов плотной симметричной матрицы А.