Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности, страница 7

. , hi+1,i )T := Gk (h1i , . . . , hi+1,i )Tувеличить kПостроить Gi так, что [Gi h·,i ]i+1 = 0g := Gi gЕсли |gi+1 | < εтоp := iвыйти из цикла по iувеличить iРешить треугольную СЛАУ H1:p,1:p y = g1:px0 := x0 +pPi=1yi viЕсли p < mто КОНЕЦr0 := b − Ax0 ; β := kr0 k2Продолжать, пока β > εРис. 5.3. GMRES(m) с вращениями Гивенса46С использованием определения вектора v 1 и соотношений (4.19),(4.21) имеемr = b − A(x0 + Vm y) = r0 − AVm y == βv1 − Vm+1 Hm y == βVm+1 e1 − Vm+1 Hm y == Vm+1 (βe1 − Hm y) .(5.4)Поскольку матрица Vm+1 составлена из ортонормированных векторов-столбцов, справедливо равенствоΨ(y) = Φ2 (x0 + Vm y) = kβe1 − Hm yk22 .(5.5)Таким образом, для нахождения коэффициентов линейного комбинирования векторов {vi } в GMRES необходимо решить СЛАУHm y = βe1 ,(5.6)которая является переопределенной (так как матрица H m имеет размерность (m+1)×m) и поэтому, как видно из (5.5), должна решатьсяв смысле наименьших квадратов.Замечание.

Соотношение (5.6), на котором основан метод, можетбыть получено и с помощью непосредственного применения проекционного подхода (4.1)–(4.5). Однако формула (4.5) в этом случаетребует более сложной интерпретации [1].Матрица системы (5.6) имеет верхнюю хессенбергову форму.

Поэтому (5.6) проще всего решать с помощью предварительного приведения к верхнетреугольному виду; последняя строка H m при этомобнуляется. Однако в отличие от FOM, где похожая на (5.6) системаявляется квадратной, в GMRES для приведения к треугольному видуиспользуются ортогональные преобразования 6 .С учетом верхней хессенберговой формы наиболее простым вариантом оказываются вращения Гивенса [2, 5, 15]. Всего необходимо mпоследовательно применяемых вращений G 1 , G2 , .

. . , Gm , причем Gkстроится так, чтобы его применение к (h kk , hk+1,k )T обнуляло второйкомпонент этого вектора.Обозначим Qm = Gm Gm−1 . . . G1 . Будучи произведением матрицвращений, Qm сама является ортогональной матрицей. Последовательное применение вращений приводит к переходу от (5.6) кRm y = g ,6(5.7)См. доказательство теоремы 5.3.1.

Ортогональность преобразований гарантирует неувеличение нормы невязки, которая при этом может быть легко оценена подобно тому, как это делается в FOM.47где Rm = Qm Hm , g = Qm (βe1 ). По аналогии с матрицей Hm , обозначим за Rm матрицу Rm без последней строки (которая по построениювращений является нулевой), а за g — вектор g без последнего коэффициента.Тогда имеет место следующая теорема:Теорема 5.3.1 Вектор y , доставляющий минимум функционалуΨ(y) в методе GMRES, определяется соотношением y = R −1m g . Приэтом норма невязки rm , соответствующей решению xm = x0 + Vm y ,равнаkb − Axm k2 = |gm+1 | .(5.8)Доказательство.

Первая часть утверждения теоремы следует из соотношения (5.4). С учетом ортогональности матрицы Q m имеемkrm k22 = kβe1 − Hm yk22 = kQm (βe1 − Hm y)k22 = kg − Rm yk22 == |gm+1 |2 + kg − Rm yk22 .В правой части первое слагаемое не зависит от y, а следовательно,минимум всей суммы достигается при y = R −1m .Для доказательства второй части теоремы воспользуемся темсвойством ортогональных матриц, что Q Tm Qm = I. То же соотношение (5.4) можно переписать в видеrm = Vm+1 (βe1 − Hm y) = Vm+1 QTm Qm (βe1 − Hm y) == Vm+1 QTm (g − Rm y) .Последний (m + 1)-й компонент вектора R m y равен нулю, а первые m, как только что было доказано, в точности совпадают с первыми m компонентами вектора g.

Поэтому g − R m y = gm+1 em+1 иrm = gm+1 Vm+1 QTm em+1 ,откуда и следует утверждение (5.8).2Как и в случае FOM, основной проблемой в GMRES является выбор размерности m пространства K; хранение базисных векторовэтого пространства определяет основные затраты памяти при программной реализации метода, составляющие O(mn) вещественныхчисел.Наиболее распространенной модификацией GMRES является перезапускаемый GMRES7 или GMRES(m). Выбирается некоторая размерность m < n, определяемая, по существу, лишь доступным объемом памяти, и после обновления решения x m = x0 + Vm y оно принимается за новое начальное приближение, после чего процесс повторяется.7Restarted GMRES.48Такое обновление после построения m базисных векторов и нахождения коэффициентов их комбинирования называется GMRESциклом, тогда как построение одного очередного базисного векторасчитается GMRES-итерацией.Внутри цикла реализуется проверка условия завершения в соответствии с (5.8); такая проверка легко реализуется, если совмещатьпостроение базиса с помощью процедуры Арнольди и применениевращений Гивенса к уже найденным элементам матрицы H m .Соответствующий алгоритм показан на рис.5.3.Очевидно, что при программной реализации нет необходимостив непосредственном хранении матриц вращения G i , так как каждаятакая матрица полностью определяется двумя числами — косинусом и синусом угла поворота — и индексами компонент вектора, ккоторым применяется вращение.Предобусловливание GMRES осуществляется точно так же, какпредобусловливание FOM.5.4.BCG: Метод бисопряженных градиентовВ §4.6 была описана биортогонализация Ланцоша, которая в отличие от ортогонализации Арнольди, использует для построения базиса экономичные трехчленные формулы.

Выберем в качестве пространства K пространство Крылова Km (r0 , A), а в качестве L — пространство Km (r̃0 , AT ), где вектор r̃0 удовлетворяет условию (r0 , r̃0 ) 6== 0.Тогда в соответствии с (4.26), если в методе FOM заменить процедуру Арнольди процедурой Ланцоша, решение будет уточняться поформулеxm = x0 + βVm T−1(5.9)m e1 ,где β = (r0 , r̃0 ). Решение СЛАУ с трехдиагональной матрицей, очевидно, удобнее всего искать с помощью метода прогонки [4].Описанная модификация FOM носит название двойственного 8метода Ланцоша. Однако подобная схема все еще требует одновременного хранения всех базисных векторов либо их перевычисленияпосле нахождения коэффициентов. Поэтому на практике обычноиспользуется другой подход.Запишем LU-разложение для матрицы T mTm = Lm Um ,8Метод называется двойственным, так как с минимальными изменениями,практически не усложняющими его схему, может быть использован для одновременного решения (1) и системы с транспонированной матрицей AT z = b̃.49Выбрать начальное приближение x0r0 := b − Ax0Выбрать вектор r̃0 , удовлетворяющий условию (r0 , r̃0 ) 6= 0p0 := r0p̃0 := r̃0Для j = 1, 2, .

. .αj := (rj , r̃j ) / (Apj , p̃j )xj+1 := xj + αj pjrj+1 := rj − αj Apjr̃j+1 := r̃j − αj AT p̃jβj := (rj+1 , r̃j+1 ) / (rj , r̃j )Если krj+1 k2 < εто КОНЕЦ (решение достигнуто)Если βj = 0то КОНЕЦ (при данном r̃0 метод не сходится)pj+1 := rj+1 + βj pjp̃j+1 := r̃j+1 + βj p̃jувеличить jРис. 5.4. BCGи определим матрицу Pm следующим образом:Pm = Vm U−1m .(5.10)Подставим это выражение в (4.22)−1xm = x0 + βVm U−1m Lm e1 == x0 + βPm L−1m e1 .(5.11)Рассмотрим теперь свойства матрицы Pm .

По аналогии с (5.10) запишемP̃m = Wm (LTm )−1 ,(5.12)Тогда имеет местоT−1−1−1(P̃m )T APm = L−1m Wm AVm Um = Lm Tm Um = Im Dm ,(5.13)где Dm есть диагональная матрица с элементамиdii = (vi , wi ).Если обозначить образующие Pm столбцы векторами pk , а столбцы P̃m — векторами p̃k , k = 1, m, то (5.13) означает, чтоi 6= j =⇒ p̃TiApj = (Api , p̃j ) = 0.50(5.14)mОпределение 5.4.1 Системы векторов {p̃ i }mi=1 и {pi }i=1 , удовлетворяющие условию (5.14), называются A-бисопряженными или, еслине возникает неоднозначности, просто бисопряженными.Очевидно, что бисопряженность эквивалентна биортогональности относительно скалярного A-произведения, а потому для нахождения векторов p̃ и p вместо соотношений (5.10) и (5.12) можнопользоваться биортогонализацией Ланцоша (4.22)–(4.25) с использованием в ней вместо скалярных произведений (2) скалярных Aпроизведений (3).Остается заметить, что присутствующее в (5.11) выражение−1Lm e1 представляет собой первый столбец матрицы L −1m , элементыкоторого могут быть найдены по формуле (A.1).

Нетрудно видеть,что эти элементы связаны рекуррентным соотношением−1[L−1m ]i+1,1 = −γi+1 [Lm ]i,1 .На рис.5.4 приведен алгоритм, построенный в соответствии сописанными соображениями. Такая вычислительная схема носитназвание метода бисопряженных градиентов 9 или BCG10 .В качестве вектора r̃0 на практике обычно выбирается r̃0 = r0 илиr̃0 = ek , где индекс k таков, что [r0 ]k 6= 0. Как и двойственный методЛанцоша, BCG может быть легко приспособлен для решения системы ATz = b̃ одновременно с (1).910Объяснение присутствия в названии слова «градиент» см.