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

Кроме того, мы будем вычислять лишь столько первых столбцов в Я, сколько необходимо для получения достаточно хорошего приближения к решению (системы Ах = 5 или задачи Ах = Лх). На практике, обычно бывает достаточно очень небольшого (в сравнении с порядком п матрицы) числа столбцов.

Используем ЯВ;разложение матрицы К, чтобы записать К = ЯН. Тогда откуда Так как Н и Л "— верхние треугольные матрицы, а С вЂ” верхняя хессенбергова, то, как легко проверить, матрица Н = НСН ~ также является верхней 316 Глава 6. Итерационные методы для линейных систем хессенберговой (см. вопрос 6.11). Иначе говоря, мы привели А к верхней хессенберговой форме посредством ортогонального преобразования с матрицей Я. (Это первая часть алгоритма вычисления собственньгх значений несимметричных матриц, обсуждавшегося в разд.

4.4.6.) Заметим, что при симметричной матрице А матрица ЯтАЯ = Н также симметрична; симметричная же матрица, имеющая верхнюю хессенбергову форму, должна быть и нижней хессенберговой, т. е. трехдиагональной. В этом случае мы будем писать ЯтАЯ = Т. Нам еще остается показать, как вычислять столбцы матрицы Ц «по одному», а не все сразу. Положим Я = [дм...,д„]. Поскольку из ЯгАЯ = Н следует АЯ = ЯН, можно приравнять столбцы у' в обеих частях матричного равенства, что дает дъ1 Аду = ~~ Ь,ддь а=1 Учитывая, что д; — ортонормированные векторы, и умножая обе части этого соотношения иа д~, получаем тв1 дтАд.=~,й дгд =Ь .для1(т(~' в=1 и 1ц+1ддд«.1. = Ад, — ~~ Ь;ддь 1=1 В результате приходим к следующему алгоритму.

Алгоритм 6.9. Алгоритм Арнолъди для (частичного) приведения к форме Хессенберга: д, = ЬДЬ|~г /" й — это число определяемых столбцов в Я и Н '/ )ог у = 1 1о Й г = Ад )ог « = 1 »о,у г = г — 6;уд1 епд 1ог Ьз»16 — — !!42 если 1» +1 . = 0 то выход дхщ = г/)» .»11 епд ~ог Векторы ди вычисляемые в алгоритме Арнольди, часто называют векторами Арнолъди. Цикл по «, в котором пересчитывается вектор г, можно интерпретировать как применение к этому вектору модифицированного алгоритма Грома — Шмидта (алгоритм 3.1): из г вычитаются компоненты в направлениях от д» до д ч в результате чего г становится ортогонален этим направлениям.

Для вычисления векторов ды ..., дь нужны Й матрично-векторных умножений с участием матрицы А и другие операции, стоимость которых составляет 0(1«~п). Если алгоритм остановить на этом месте, то что можно узнать о 317 6.6. Методы крыловского подаространстаа матрице Ау Положим Я = [Яг, сг'„], где (вг = [В, °, дг] и сге = [чг+и, Яв]. Заметим, что мы вычислили лишь Яг и ог+1, прочие столбцы матрицы Я„не известны. Имеем Н= Я АЯ = [Юг,Юв] А[Як,Чв] = ~ фАд ~фАд" т т [ ДтАЯг ЯтАЯ„ й и — й — ь [в. й. (6.30) о1 А Т= дв-1 А~-1 ое Приравнивая столбцы г в обеих частях равенства АЯ = ЯТ, получаем Аод = Я яд 1 + од ду + Цгй.+и Умножая обе части этого соотношения на 6т и учитывая, что столбцы в мат- рице Я ортонормальны, находим дтАа = о .

Сказанное обосновывает следу- ющую версию алгоритма Арнольди, называемую алгоритмом Ланцоша: Алгоритм 6.10. Алгоритм Ланцоша для [частичного) приведения к симме- тричной трехдиагональной форме. 6~ = Ь)'[]Ь[)г, ~3о = О, Чо = 0 )ог 1 = 1 1о й в=Ад. о=де г = г — одЦ1 — Ру — 1Ц-1 61 = И[г если Д = 0 то выход Чг~- = l 61 епд )ог Отметим, что Нг — верхняя хессенбергова матрица„поскольку вся матрица Н обладает этим свойством. По той же причине, матрица Нг„имеет единственный (возможно) ненулевой элемент, стоящий в ее правом верхнем углу, а именно элемент йг+г г. Итак, Н„и Нег не известны, мы знаем лишь блоки Нг и Нг„. Если А — симметричная матрица, то Н = Т симметрична и трехдиагональна. Алгоритм Арнольди значительно упрощается, поскольку большинство чисел йб являются нулями.

Положим 318 Глава б. Итерационные методы для линейных систем Векторы вн вычисляемые в алгоритме Ланцоша, часто называют векторами Ланцоша. Вот что мы знаем о матрице А после й шагов этого алгоритма; Т вЂ” дтАд — Р„д )тАДь д ] ятА(~ ятА(~ ~ АЯь Я Ац„ й й и — й и й Т Т (6.31) Поскольку матрица А симметрична, нам известны блоки Ть и Ть„= Тты но не блок Т„. Матрица Тьа имеет единственный (возможно) ненулевой элемент, стояший в ее правом верхнем углу, а именно элемент Ды Отметим, что этот элемент неотрицателен, поскольку вычисляется как норма вектора х.

Введем некоторые стандартные обозначения, связанные с частичным разложением матрицы А, вычисляемым в алгоритмах Арнольди и Ланцоша. Определение 6.16. Подпространство Крылова Кь(А, Ь) — это линейная обо- лочка векторов Ь, АЬ, АаЬ, ..., А" 1Ь. Мы будем писать Ке вместо Кь (А, Ь), если А и Ь ясны из контекста. Если не происходит досрочного выхода из алгоритма вследствие и = О, то векторы 9ы вычисляемые в алгоритмах Арнольди и Ланцоша, образуют ортонормированный базис подпространства Кы (Можно показать, что Кь имеет размерность й в том и только том случае, если алгоритм Арнольди или Ланцоша может вычислить вектор дь без досрочного выхода; см.

вопрос 6.12.) Будем называть Нь (или Ть) проекцией матрицы А на подпространство Крылова Кь. Наша цель — построить алгоритмы решения системы Ах = Ь, использующие лишь информацию, вычисляемую на й шагах алгоритма Арнольди или Ланцоша. Мы рассчитываем на то, что й может быть много меньше, чем и, что сделает алгоритмы эффективными. (В гл. 7 та же информация будет использована для вычисления собственных значений матрицы А. Мы можем уже сейчас вкратце описать, как это будет сделано. Заметим, что если элемент йет1 ь окажется нулем, то матрица Н (или Т) будет блочно верхнетреугольной, а потому все собственные значения блока Нь будут в то же время собственными значениями для Н, следовательно, и для А, поскольку Н и А подобны. Правые собственные векторы блока Нь определяют собственные векторы матрицы Н; умножив их на Яь, мы получим собственные векторы матрицы А. Если элемент й~+~ ь не равен нулю, но мал, то можно ожидать, что собственные значения и собственные векторы блока Не будут хорошими приближениями к собственным значениям и собственным векторам матрицы А.) Закончим это введение указанием на то, что при наличии округлений многие из обсуясдаемых алгоритмов ведут себя совершенно иначе, чем они делали бы это в точной арифметике.

В частности, векторы оы вычисляемые в алгоритме Ланцоша, могут быстро потерять ортогональность и, в действительности, 319 б.б, Методы крыловского оодоространства часто становятся линейно зависимыми. Эта кажущаяся катастрофической численная неустойчивость была причиной того, что методы данной группы были иа ряд лет преданы забвению после их обнаружения. Однако со временем исследователи научились стабилизировать эти алгоритмы или убедились, что сходимость в них имеет место вопреки неустойчивости! Мы вернемся к этим вопросам в разд. 6.6.4, где анализируется сходимость метода сопряженных градиентов для решения системы Ах = Ь [который «неустойчив», но тем не менее сходится), и в гл.

7, особенно в рэзд. 7.4 и 7.5, где показано, как вычислить собственные значения (и основной алгоритм модифицирован так, чтобы обеспечить устойчивость). 6.6.2. Решение системы Ах = Ь с помощью крыловского надпространства Как же решать систему Ах = Ь, имея лишь информацию, доступную после к шагов алгоритма Арнольди или Ланцоша? Поскольку единственными известными нам векторами являются столбцы матрицы Я», единственное место, где следует «искать» приближенное решение, — это надпространство Крылова Кы натянутое на векторы йв Иначе говоря, мы будем искать «наилучшее» приближенное решение вида где к = [кы.,.,к»] т Теперь нужно определить смысл слона «наилучшее».

Имеется несколько естественных, но различных определений, ведущих к разным алгоритмам. Будем обозначать точное решение через х = А гЬ, а невязку Ь вЂ” Ах» через г». 1. «Наилучшее» х» минимизирует []хь — х[]ю К сожалению, у нас нет достаточной информации в крыловском подпространстве для вычисления такого вектора хы 2. «Наилучшее» хь минимизирует ]]г»[[ю Это определение вполне реализуемо и соответствующий алгоритм называется М1НВЕЯ [от ппп1тит ге»Ыиа!— минимальная невязка) в случае симметричной матрицы А [194] и СМВЕЯ (от уепега!«»ей т?п»тит ге»Ыиа! — обобщенный алгоритм минимальных невязок), если матрица А несимметрична [215]. 3. «Наилучшее» хь обеспечивает условие г».!.Кы или Я» гь = О.

Оно иногда называется условием ортогон линой невлзки, или условием Галеркина по аналогии с родственным условием в теории конечных элементов. Если А симметрична, то соответствующий алгоритм называется ЯУММЩ [194]. Для несимметричной матрицы А имеется вариант метода СМКЕБ [211]. 4. Посредством симметричной положительно определенной матрицы А можно задать норму ][с[]л- — — (г~А 'г)'7» (см. лемму 1.3). Будем считать, что «наилучшее» хь минимизирует [[г»[[л- .

Эта норма есть то же самое, что ][х» — х[]л. Соответствующий метод называется алгоритмом сопряженных градиентов [145]. Для симметричной положительно определенной матрицы А два последних определения слова «наилучшее» оказываются эквивалентными. 320 Глава 6. Итерационные методы дли линейных систем Теорема 6.8. Пусть А — симметричнол матрица, Тг = ЯгтА1гг и тг = 6— Ах«, где х б Кг. Если матрица Тг невырожденна и хг = ЯгТг ег!!Ь|!«, где е~" = [1,0,...,0]т, то Я~~тг = О. Если матрица А егце и положительно определена, то Тг обязана быть невырожденной матрицей и указанный выбор для хг минимизирует !!тг!!и-1 по всем хг е Кг.

Кроме того, тг = ~!!тг!!«угь1. Доказательство. Чтобы упростить обозначения, опустим индекс (г. Предположим, что матрица Т = Я~АД невырожденна, н пусть х = ЯТ 'е1!!Ь!!г и т = 6 — Ах. Проверим, что Яхт = 0 следующим вычислением: дтт дт(6 Ах) дтЬ дтАх = е1!!Ь!!г — С«~АЯТ ~еЯ!Ь!!г) поскольку первым столбцом матрицы Я является вектор Ь/!!Ь!!г, а прочие столбцы ортогональны к Ь = е1!!Ь!!г — (9~АД)Т 1е1!!Ь|!г = е1!!ЬЦ!г — (Т)Т 1е1!!Ь!!г поскольку Я~АД = Т = О. Предположим теперь, что А к тому же положительно определена.