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

вопрос 5.13). Мы докажем для этого случая, что сходимость почти всегда является кубической. Теорема 5.9. НО.-итерации локально сходится кубически, т. е. если собственное значение простое, а ошибка достаточно мала, то число верных разрядов приближенного собственного значения утраивается на каждом шаге. Доказательство.

Мы утверждаем, что достаточно проанализировать случай диагональной матрицы А. Действительно, пусть ц' Ац = Л, где ц — ор- тогональная матрица, составленная из собственных векторов матрицы А, а Л вЂ” диагональная матрица собственных значений. Заменим переменные в В14- итерации по формулам х; = Я~хг и у; = Я~у,. Тогда т -т т - -т и Яу;+з — — 1А — р;1) гахн поэтому у;~.~ = Я~(А — р;1) ~Яхг = Я~ АЦ вЂ” р,1) ~хг = (Л вЂ” рг1) 'х;. Следовательно, проведение КЯ-итерации с А и хо равносильно ее проведению с Л и хо.

Поэтому, без ограничения общности, предположим, что А = Л— уже диагональная матрица; таким образом, собственными векторами являются столбцы е; единичной матрицы. Не ограничивая общности, будем считать, что векторы х; сходятся к еы Тогда можно написать х, = е1+ д„где 5о;!)з = е « 1.

Для доказательства кубической сходимости нужно показать, что хьь1 — — ед + 4-ьы где ()дььз'вз = О~ з) Заметим прежде всего, что т ( +д)т1 +яг) т +2 т1 + 1т1 1+2д + з поэтому дп — — — ез/2. Следовательно, р; = хТЛх; = 1ез + дг) Л(ег+ дз) = е~~Лез + 2е~~Лд;+~йЛдг = а1 — у, где г1 = — 2езЛдг — о," Лдз = сг1е~ — атЛдг. Отсюда ф < 1сгз!е + 5Лйзйа15з < 25Лйзе . 228 Глава 5. Симметричная проблема собственных значений Таким обРазом, число Р, = аг — 21 = аг + 0(ег) ЯвлЯетсЯ очень хоРошим пРиближением к собственному значению аы Теперь можно написать (Л вЂ” р 1) ~хг ! т хп Х22 хс« 1 поскольку (Л вЂ” р,1) ~ = бба8 1,«»1 р' / ! 1+ 41 42 Ае поскольку х, = ег + бг ! 1 — ег/2 022 б;„ 21 аг ог+7/ оп йг+гг поскольку р; = аг — 9 и дгг = — е /2 г 1 — ег/2 [ 41221 4 21 [ ' (1 — ег/2)(аг — аг + 9) ' ' (1 — е~/2)(ое — ог + 21) 1 — ег/2 (е, + 22„.,).

21 Чтобы оЦенить [[4~2 $$2, использУем оЦенкУ ДлЯ знаменателей [оу — аг + 21$ > 8ар(1, Л) — $21$ и неравенство (5.12). В результате получим $$4$$2$9$ < 2$[Л [[2«2 (1 — ег/2)(8ар(1, Л) — $21$) (1 — ег/2)(бар(1, Л) — 2$[Л[[ег) ' или $$4»2$$2 — — 0(ег). НаконеЦ, посколькУ х»ю — — ег + с(»»г — — (ег + «[2~-г)/$[ег + 4~-2$$2~ то и $$4~-2$$2 = 0(е ). П 5.3.3. «Разделяй-и-властвуй» Это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матрицы, начиная с порядка п, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна.

В самом деле, хотя впервые метод был предложен еще в 1981 г. [59$, «правильный» способ его реализации был найден лишь в 1992 г. [127, 131]. Этот способ воплощен 1 АРАСК-программами веуечб (для плотных матриц) и весечб (для трехдиагональных матриц). В них стратегия «разделяй-и-властвуй» используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к Я11-итерации. 5.3. Алгоритмы длл симметричной проблемы собственных значений 229 Обсудим вначале структуру алгоритма в целом, а уже затем численные детали. Пусть Предположим, что нам известны спектральные разложения матриц Т, и Тас Т; = Я;ЛЯ~.

В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матри- 230 Глава 5. Симметричная проблема собственных значений цы Т и собственными значениями матриц Т~ и Тг. Имеем г=- [ т, 01 т О та+5 2 ,л, г ]+ """ о Ог л, +'""" о дт где С Х~ 0 ] [последний столбец матрицы Ят О СХг ] [ первый столбец матрицы Я так как и = [О,...,О, 1, 1,0,...,0]т.

Следовательно, Т имеет те же собственные значения, что и подобная ей матрица Р + рии, где Р = [ О Л ]— л, о Лг диагональная матрица, р = 5 — число, а и — вектор. Будем предполагать, не ограничивая общности, что диагональные элементы с(ы ..., Ы„матрицы Р упорядочены по убыванию: д„< ° . < 4.

Чтобы найти собственные значения матрицы Р + риит, вычислим ее характеристический многочлен, считая пока матрипу Р— Л1 невырожденной. Тогда де1(Р+ рииг — Л1) = де1((Р— Л1)(1+ р(Р— ЛХ) гиит)). (5.13) Поскольку Р— ЛХ невырожденна, бег(1+ р(Р— ЛХ) 'иит) = 0 тогда и толь- ко тогда, когда Л вЂ” собственное значение.

Заметим, что матрица 1+ р(Р— Л1) 'ии~ получается из единичной добавлением матрицы ранга 1. Определи- тель такой матрицы легко вычислить. Доказательство леммы составляет содержание вопроса 5.14. Таким образом, де1(1+ р(Р— Л1) 'ии ) = 1+ ри (Р— Л1) 'и =1+р~ ' =- Х(Л), г=1 (5.14) т.е. собственные значения матрицы Т суть корни так называемого векового уравнения 1(Л) = О. Если все числа д, различны и все и; ф 0 (случай общего положения), то 1(Л) имеет график типа показанного на рис.

5.2 (где и = 4 ир>0). Можно видеть, что прямая у = 1 является горизонтальной асимптотой для этого графика, а прямые Л = д; суть вертикальные асимптоты. Поскольку г 1'(Л) = Р 2,"., ~-"-',)т > О, фУнкциЯ возРастает всюдУ, кРоме точек Л = дь Поэтому корни функции разделяются числами д, и еще один корень находится справа от точки дг (на рис. 5.2 дг = 4). (При р < 0 функция 1(Л) всюду убывает Лемма 5.1. Справедливо равенство бег(1 + хут) = 1 + утх, где х и у— векторы. 5.3. Алгоритмы для симметричной проблемы собственных значений 231 2 3 Рис. 5.2. График функции 1(Л) = 1+ 1"'"х + з 'х + з 'х + 4'х.

и соответствующий корень находится слева от точки с(„.) Для функции Х(Л), монотонной и гладкой на каждом из интервалов (4+~,4), можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в (4~.ы 4). Мы обсудим детали позднее в данном разделе. Пока нам достаточно знать,что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление Х(Л) и Х'(Л) стоит 0(п) флопов, для вычисления одного собственного значения достаточно 0(п) флопов, а для вычисления всех п собственных значений матрицы Р + рии~ требуется 0(пз) флопов.

Для собственных векторов матрицы Р + риит мы легко можем получить явные выражения. Лемма 5.2. Если а — собственное значение матрицы Р+рии, то соответствующий собственный вектор равен (Р— а1) 1и. ХХоскольку матрица Р— а1 диагональная, для вычисления такого вектора достаточно 0(п) флопов. Доказательство. (Р+ рий)((Р— а1) ~и) = (Р— аХ+ а1+ рии )(Р— аХ) ги = и+ а(Р— а1) 1и+ и~рит(Р— а1) 1и) = и + а(Р— аХ) 1и — и поскольку ри (Р— аХ) 'и+1 = Х(а) = О = а[(Р— аХ) 'и), что и требовалось. П 232 Глава 5. Симметричная проблема собственных значений Для вычисления по этой простой формуле всех п собственных векторов требуется 0(п') флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений «г, может давать неортогональные приближенные собственные векторы и;.

Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе. Алгоритм в целом является рекурсивным. Алгоритм 5.2. Вычисление собственных значений и собственных векторов симметричной трехдиагональной матрицы посредством стратегии «разделяй-и-властвуй»н ртос дс егд(Т, 1з, Л) ... по входной матрице Т вычисляются выходные матрицы 9 и Л, такие, что Т = 1еЛ1е« Ясли Т вЂ” матрица размера 1 х 1 присвоить выходным параметрам значения 0 = 1, Л = Т иначе представить Т в виде Т = ~ О Т ~ + Ь ио г са11 дс ег д(Тг, б»1 ы Л г ) сай Йс е1д(Тг, Юг Лг) построить Р+ рии~ по Лы Лг, 1»1ы ссг найти матрицу собственных значений Л и матрицу собственных векторов Я' для матрицы Р + риит построить матрицу собственных векторов 1е для матрицы Т: 1д, присвоить выходны.м параметрам значения 1е и Л епдг/ Проанализируем сложность алгоритма 5.2.

Пусть 1(п) — число флопов при обработке матрицы размера и х и процедурой дс е15. Тогда 1(и) = 21(п/2) два рекурсивных обращения к «1с е!5(То Ян Л,) +0(пг) вычисление собственньгх значений матрицы .Р + рии +0(пг) вычисление собственных векторов матрицы Р+ рии +с. пг вычисление матрицы 1д = (д' ('„1 О О Чг Если Цы 1дг и ц' рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа с в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоемкую часть алгоритма в целом. Игнорируя члены порядка иг, получаем г(и) = 2с(и/2) + сиз. Решая это рззностное уравнение, находим с(п) с1пз (см. вопрос 5.15).

На практике константа с обычно гораздо меньше 1, потому что матрица 1'„1' весьма разрежена вследствие явления, называемого дефляцией. Мы обсудим дефляцию в следующем разделе, а затем перейдем к деталям, связанным с решениен. векового уравнения и устойчивым вычислением 5.3. Алгоритмы для симметричной проблемы собственных значений 233 собственных векторов. Наконец, мы поговорим о том, как ускорить метод с помощью техники РММ, используемой при моделировании электростатики заряженных частиц [124]. При первом чтении книги эти разделы можно опустить.