Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 73
Текст из файла (страница 73)
Тогда н Т должна быть положительно определенной, а потому невырожденной (см. вопрос 6.13). Пусть х = х+ Я« — другое возможное решение из К и пусть т = 6 — Ах. Нужно показать, что !!т!!и-1 достигает минимума при « = О. Однако !!Р!!гл, — — ттА 'т по определению = (т — АЯ«)тА 1(т — АЯ«) так как т = Ь вЂ” Ах = Ь вЂ” А(х + О«) = т — А(г« = ттА 1тт 2(АГ~«)тА — 1т + (АЯ«) тА — 1(АЯ«) !!т!!г 1 2«тдтт + /!АЯ«!!2А г так как (Аб)«)тА " = «т~т 4А " = «гатт !!т!!л — г + !!АЯ«!!и — 1 так как С3 т = О, позтому значение !!т!!л-1 минимально в том и только том случае, если АО« = О. Так как А невырожденна, а Я имеет полный столбцевой ранг, то соотношение АЯ« = 0 равносильно равенству « = О.
Восстанавливая индексы, покажем, что тг = Щтг!!«да+и Поскольку хг й Кг, имеем тг = 6 — Ахг й Кгг ы Таким образом, тг есть линейная комбинация столбцов матрицы Ягг ы потому что Кг.г1 натянуто на зти столбцы. Учитывая, что Яг тг = О, заключаем, что дг+1 — это единственный столбец матрицы Яг+и т к которому вектор тг не ортогонален. П 6.6.3. Метод сопряженных градиентов Наиболее предпочтительным методом для симметричных положительно определенных матриц является алгоритм сопряженных градиентов. Теорема 6.8 характеризует решение хг, вычисляемое этим алгоритмом.
Хотя М11чВЕЯ мо- 321 6.6. Методы крыловского подпрострапства жет выглядеть более естественным методом, чем СО, поскольку минимизирует ||ге||э, а не ||та||я-ы оказывается, что его реализация более трудоемка, он более подвержен численной неустойчивости, а потому часто дает приближенные решения худшего качества, чем СО. Мы увидим, что СС имеет то особенно привлекательное свойство, что в его реализации предусмотрено одновременное хранение в памяти лишь четырех векторов, а не й векторов оы ..., оь.
Кроме того, в его внутреннем цикле, помимо матрично-векторного произведения, выполняются только два скалярных произведения, три операции типа «эахру» (прибавление кратного одного вектора к другому) и небольшое количество скалярных операций. Таким образом, и память, и работа в методе очень невелики. А теперь мы выведем формулы метода. Это можно сделать несколькими различными способами. Мы начнем с алгоритма Ланцоша (элгоритм 6.10), вычисляющего столбцы ортогональной матрицы Яь и элементы треддиагональной матрицы Ты а также с формулы хь = бдьТе ~е1 ||Ь||г нз теоремы 6.8. Будет показано, как можно вычислить хь непосредственно, пользуясь рекуррентными соотношениями для трех групп векторов.
Одновременно в памяти хранятся лишь самые последние векторы из каждой группы, записываемые на места своих предшественников. Первая группа векторов — это приближенные решения хы Вторая группа — это невязки те = Ь вЂ” Ах»; согласно теореме 6.8, они параллельны векторам Ланцоша аь+ы Третья группа векторов — это сопряженные градиенты ры Эти векторы называются градиентаами по следующей причине: шаг метода СС может интерпретироваться как выбор числа о из условия, чтобы новое приближенное решение хь = хь» + прь минимизировало норму невязки ||т~Д|л- = (тьтА»ть)»~г. Иными словами, рь используются как нанравлени градиентного поиска.
Они называются сопряженными, или, точнее, А-сопряженными, потому что рьтАр = 0 при»' ~ й. Другими словами, векторы рь ортогонэльны по отношению к скалярному произведению, определяемому матрицей А (см. лемму 1.3). Поскольку А — симметричная положительно определенная матрица, то же самое верно для Ть = ЦтьА(;1» (см. вопрос 6.13). Это означает, что к Ть можно применить алгоритм Холесского; в результате получим Ть = ТьР»Т», где т Т,ь — нижняя двухдиагональная матрица с единичной диагональю, а Рь — диагональная матрица. Далее, используя формулу для хь из теоремы 6.8, находим хь = ЯеТь ~е»||Ь!|г = Сгг(1, ~Р„~й ~)е»||Ь|~р = Ягель ~)(Р„йь ~е»||Ь||г) = (Рь)(уь), где Рь = ЯьТ~ и уь = Р„й„е»||Ь||г.
Положим Рь = (рд,...,рь). Сопряженные градиенты р, окажутся параллельны столбцам р, матрицы Ры Теперь мы знаем достаточно для того, чтобы доказать следующую лемму: Лемма 6.8. Столбцьг р; матрицы Рь являются А-сопряоюенными. Иначе го- воря, матрица РьтАР» — диагональная. Глава 6. Итервциоииые методы для линейных систем Выведем теперь простые рекурсии для столбцов матриц Рь и элементов векторов дь. Мы покажем, что уь 1 = [111,...,11ь 1]т совпадает с первыми Й вЂ” 1 компонентами вектоРаУь = [т11,...,Ць 1, тая]т, а Рь 1 совпаДает с пеРвыми Й вЂ” 1 столбцами матрицы Рь.
Поэтому соотношение хь = Рь ° уь = [Рь „рь] ° =Рь 1уь 1+рь11ь = хе 1+рь111 (6.32) УЬ вЂ” 1 01 может служить искомой рекурсией для векторов хь. Рекуррентное соотношение для величин пь выводится следующим образом. Поскольку Ть 1 является ведущей главной подматрицей порядка Й вЂ” 1 в Ть, матрицы Хь 1 и Х1ь 1 также будут ведущими главными подматрицами порядка Й вЂ” 1 соответственно в Хь и .Оь: о А А-1 А,1 ОЬ 1ь-1 1 т -т 1 ° йаь(Х1ь 1, ах) 1, ет где размерность вектора еть, — — [О,...,0,1] равна Й вЂ” 1.
Точно так же, Х1ь и Х..~1 суть ведущие главные подматрицы порядка Й вЂ” 1 соответственно в Ц,1 = 61ай(0 11, 11„1) и Х-1 ~ь Выражения элементов последней строки матрицы Х~1 не существенны для нас. Из сказанного следует, что вектор уь 1 = Х1,, 1 Хь ~1е1[]Ь][1, где е1 имеет размерность Й вЂ” 1, совпадает с первыми Й вЂ” 1 компонентами вектора -1 уь = Х1, 'Х, "е1][Ь][з = А. Пь1 Х,,~ е1[]Ь]]т ~ ~ уь 1 11ь '1ь е1[)Ь|]т Докавательси1во.
Имеем РтЛР =[О Х-т)'А[О Х- ) =Х-'Д Ад,)Х- =Х-1[т,)Х- Х-1[Х Х1 Хт)Х-и Х1 П б.б. Методы крыловского подпространствв 323 Теперь нам нужна рекурсия для столбцов матрицы Рв = (Рм ...,рь]. Поскольку матрица 1ь „— верхняя треугольная, то же самое верно в отношении т матрицы 7~ 1, являющейся в 7 ь ведущей главной подматрицей порцдка Й вЂ” 1. Поэтому Рь в совпадает с первыми 1с — 1 столбцами матрицы — т О ~-т (О ! ( ь, * (г) ~-т 0 1 ~ Соотношение Рк = ЩЕ,, т дает Рай~~ = Яы откуда, приравнивая й-е столбцы обеих частей, получаем рекурсию Рк = Ф вЂ” 1ь-1Рь-ы (6.33) Итак, мы имеем рекурсии для векторов дь (из алгоритма Ланцоша), векторов рь (нз уравнения (6.33)) и приближенных решений яь (из уравнения (6.32)). Все зти рекурсии являются короткими, т.
е. для реализации каждой из них нужны лишь один или два предыдущих вектора. Поэтому вместе они дают способ вычисления приближенного решения хь при хранении в памяти небольшого числа векторов и выполнении во внутреннем цикле малого числа скалярных произведений, операций типа вахру и скалярных операций. Чтобы получить алгоритм СС в окончательном виде, нужно еще немного упростить эти рекурсии. Так как, согласно теореме 6.8, векторы гь и дв.ь~ параллельны, ланцошеву рекурсию для де+1 можно заменить рекурсией гь = Ь вЂ” Ахь или, что эквивалентно, рекурсией гь = гь 1 — пвАрь (получаемой умножением соотношения кь = ть 1 + парь на А и последующим вычитанием из тождества Ь = Ь). В результате получаем векторные рекурсии гь = гк 1 — цьАры хь = хь 1 + перь из уравнения (6.32), Рь = бь — 1ь-1рь-~ из уравнения (6.33). (6.34) (6.35) (6.36) Чтобы исключить здесь ды сделаем подстановки дь = гь-1Дгь-1!!а и Рв = !!гу, в//эры что даст — Арк Чь !!гь-1!!2 — ил Ары + чь !! -!! Рь + пвры !!га д!!т1ь !! ь- !! + рь 'Рь — ы гь = гв (6.37) гк — 1 (6.38) хь Рь = гь-1 (6.39) гь Нам еще нужны формулы для чисел иь и рв.
Как мы увидим, для них имеется несколько математически эквивалентных выражений через скалярные произведения векторов, вычисляемых алгоритмом. Наши окончательные формулы выбраны так, чтобы минимизировать число используемых скалярных произведений, а также потому, что они более устойчивы, чем альтернативные варианты. 324 Глава 6. Итерационные методы для линейных систем ра Арь = р„Ать, + 0 = т„,Арь. т т т (6.40) Далее, умножим обе части уравнения (6.37) слева на тат и используем соотношение тат,ть = 0 (напомним, что векторы т; параллельны столбцам ортогональной матрицы 9).
Получим т„ т т т„, Ара т т ть ,тА,„ (6.41) согласно уравнению (6.40) (Уравнение (6.41) можно было бы вывести и из свойства числа иы указанного в теореме 6.8, а именно: иь минимизирует норму невязки )(ть))л, тт~А ~т~ = (ть 1 — иьАрь) А ~(ть 1 — иаАрь) согласно уравнению (6.37) т -г т 2 Т = та А ть 1 — 2иьрьта 1 + иьра Арь. Это выражение есть квадратичная функция от иь, поэтому минимум легко находится приравниванием производной нулю и решением полученного урав- нения относительно иь. Это дает р,ть т 1тьАра + „„„тт согласно уравнению (6.39) ртАрь ть ть т рь Арь Здесь было использовано соотношение ра,ть 1 = О, вытекающее из того, что т та г ортогонален ко всем векторам из 1Са ы включая вектор ра ы) Чтобы получить формулу для ~ц„умножим обе части уравнения (6.39) слева на рьт А и используем А-сопряженность векторов ра и рь 1 (лемма 6.8).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
