Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности (1051191), страница 4
Текст из файла (страница 4)
Такая перестановка дает обратный метод последовательной верхней релаксации:(D − ωF)xk+1 = [ωE + (1 − ω)D]xk + ωb .(3.15)Последовательное применение прямого и обратного методов SORдает симметричный метод последовательной верхней релаксации(SSOR5 )(D − ωE)xk+1/2 = [ωF + (1 − ω)D]xk + ωb ;(D − ωF)xk+1 = [ωE + (1 − ω)D]xk+1/2 + ωb .3.3.Связь с предобусловливаниемРанее было показано, что методы Якоби и Гаусса-Зейделя, а такжеSOR и SSOR могут быть представлены в виде (3.8).
Другой формойэтого выражения являетсяxk+1 = Gxk + f ,(3.16)которое получается из (3.8) домножением левой и правой частей наK−1 .Проведя аналогию с (3.1) и (3.2), нетрудно увидеть, что (3.16)можно рассматривать как итерационную схему, построенную длярешения СЛАУ(I − G)x = f ,(3.17)в которойf= K−1 b ;G = K−1 R = K−1 (K − A) = I − K−1 A .Таким образом, система (3.17) эквивалентна системеK−1 Ax = K−1 b ,5Symmetric SOR.25т.е.
предобусловленной СЛАУ (1) с матрицей предобусловливанияK. Для рассмотренных в данной главе методов эти матрицы равны 6 :KJ = D ;KGS = D − E ;1(D − ωE) ;KSOR =ω1(D − ωE)D−1 (D − ωF) .KSSOR =ω(2 − ω)1Скалярные множители ω1 и ω(2−ω)при практических реализацияхSOR и SSOR-предобусловливания обычно не учитываются, так какдают лишь общее масштабирование, практически не влияющее наскорость сходимости.6Выражения для KGS и KSOR приведены для прямых методов. Для обратныхметодов матрицы E и F меняются местами.26Глава 4Проекционные методы.Подпространства Крылова4.1.Общий подход к построению проекционныхметодовРассмотрим систему Ax = b и сформулируем для нее следующую задачу.
Пусть заданы некоторые два подпространства K ⊂ R n и L ⊂Rn . Требуется найти такой вектор x ∈ K, который обеспечивал бы решение исходной системы, «оптимальное относительно подпространства L», т.е. чтобы выполнялось условие∀l ∈ L :(Ax, l) = (b, l) ,называемое условием Петрова-Галеркина. Сгруппировав обе частиравенства по свойствам скалярного произведения и заметив что b −− Ax = rx , это условие можно переписать в виде∀l ∈ L :(rx , l) = 0 ,(4.1)т.е. rx = b − Ax ⊥ L. Такая задача называется задачей проектирования решения x на подпространство K ортогонально к подпространству L.В более общей постановке задача выглядит следующим образом.Для исходной системы (1) известно некоторое приближение x 0 к решению x∗ .
Требуется уточнить его поправкой δ x ∈ K таким образом,чтобы b − A (x0 + δx ) ⊥ L. Условие Петрова-Галеркина в этом случаеможно записать в виде∀l ∈ L :(rx0 +δx , l) = ((b − Ax0 ) − Aδx , l) = (r0 − Aδx , l) = 0 .Пусть dim K = dim L = m. Введем в подпространствах K и L базиmсы {vj }mj=1 и {wj }j=1 соответственно. Нетрудно видеть, что (4.1) вы27полняется тогда и только тогда, когда∀j (1 ≤ j ≤ m) :(r0 − Aδx , wj ) = 0 .(4.2)При введении для базисов матричных обозначений V = [v 1 |v2 | . . .
|vm ]и W = [w1 |w2 | . . . |wm ] можно записать δx = Vy где y ∈ Rm — векторкоэффициентов. Тогда (4.2) может быть записано в видеWT(r0 − AVy) = 0 ,(4.3)откуда WTAVy = WTr0 иy = WTAV−1W Tr 0 .(4.4)Таким образом, решение должно уточняться в соответствии сформулойx1 = x0 + V WTAV−1WTr0 ,(4.5)из которой сразу вытекает важное требование: в практических реализациях проекционных методов подпространства K и L и их базисы должны выбираться так, чтобы матрица W TAV либо быламалой размерности, либо имела простую структуру, удобную дляобращения.Из (4.3) также вытекает соотношениеVy = A−1 r0 = A−1 (b − Ax0 ) = x∗ − x0 ,т.е. Vy = δx представляет собой проекцию на подпространство K разности между точным решением и начальным приближением.Пусть имеется набор пар подпространств hK i , Li iqi=1 таких, чтоK1 ⊕ K2 ⊕ .
. . ⊕ Kq = Rn и L1 ⊕ L2 ⊕ . . . ⊕ Lq = Rn . Тогда очевидно, чтопоследовательное применение описанного ранее процесса ко всем таким парам приведет к решению xq , удовлетворяющему исходной системе Ax = b. Соответственно, в общем виде алгоритм любого методапроекционного класса может быть записан следующим образом:НачалоВыбор пары подпространств K и LПостроение для K и L базисов V = [v1 |v2 | .
. . |vm ]и W = [w1 |w2 | . . . |wm ]r := b − Ax−1y := WTAVW Trx := x + VyПродолжать до достижения сходимости284.2.Случай одномерных подпространств K и LНаиболее простой ситуацией является случай, когда пространства Kи L одномерны. Пусть K = span{v} и L = span{w}.
Тогда (4.5) принимает видxk+1 = xk + γk vk ,(4.6)причем коэффициент γk легко находится из условия ортогональности rk − A(γk vk ) ⊥ wk :(rk − γk Avk , wk ) = (rk , wk ) − γk (Avk , wk ) = 0 ,откуда γk = (rk , wk ) / (Avk , wk ).П РИМЕР 1. Выберем vk = wk = rk . Тогда (4.6) примет видxk+1 = xk +(rk , rk )rk .(Ark , rk )(4.7)Поскольку выражение в знаменателе представляет собой квадратичную форму rkTArk , сходимость процесса гарантирована, если матрица A является симметричной и положительно определенной.Данный метод есть метод наискорейшего спуска 1 ; можно показать что на каждой итерации в нем минимизируется функционалΦ(x) = kx − x∗ k2A в направлении −∇Φ(x).Действительно, в силу положительной определенности A, квадратичная форма (x − x∗ )TA (x − x∗ ) = Φ(x) достигает своего минимума (равного нулю) при x = x∗ и строго выпукла. При этомΦ(x) = (A(x − x∗ ), x − x∗ ) == (Ax, x) − (Ax, x∗ ) − (b, x) + (b, x∗ ) == xTAx − xT∗ Ax − bTx + bTx∗ .В силу симметричности матрицы A справедливо x T∗ Ax = bT A−1 Ax == bTx, и функционал равенΦ(x) = xTAx − 2bTx + bTx∗ ,(4.8)а его градиент ∇Φ(x) = 2Ax − 2b = −2rx ; следовательно, −∇Φ(x) ↑ rx .Покажем теперь, что выбор γ = krx k/(rxTArx ) доставляет минимум Φ(x) в этом направлении.
Подставим в (4.8) выражение x + γr;последним слагаемым bTx∗ при этом можно пренебречь, так как онопостоянно и не влияет на процесс минимизации. Снова будем учитывать симметрию матрицы Af (γ) = Φ(x + γr) = (x + γr)TA(x + γr) − 2bT(x + γr) =1Англ. Steepest Descent Method — SDM.29= xTAx + 2γxTAr + γ 2 r TAr − 2bTx − 2γbTr ==hihixTAx − bTx − bTx + 2γ xTAr − bTr + γ 2 r TAr == −(r + b)Tx − 2γr Tr + γ 2 r TAr .Найдем minγ f (γ); учитывая выпуклость Φ(x), для этого достаточно найти стационарную точку f (γ):f 0 (γ) = 2γr TAr − 2r Tr = 0=⇒r Tr(r, r)=,Tr Ar(Ar, r)γ=что совпадает с (4.7)Замечание. В практических задачах метод наискорейшего спуска обладает достаточно медленной сходимостью, соответствующийанализ приведен в [4, 3].П РИМЕР 2.
Выберем vk = ATrk и wk = Avk . Формула (4.6) принимает видxk+1 = xk +rk , AATrk(AATrk , AATrk )= xk + ATrk , ATrkATAATrk,ATrkATrk = ATrk .(4.9)Данный метод есть метод наискорейшего уменьшения невязки 2 ;условием его сходимости является невырожденность матрицы A.Сравнивая (4.9) и (4.7), нетрудно убедиться что метод наискорейшего уменьшения невязки совпадает с методом наискорейшего спуска, примененным к системе ATAx = AT b.
Тогда на основании соотношения (4) можно утверждать, что в методе (4.9) на каждом шагеминимизируется функционал Φ(x) = kb − Axk 22 по направлению −−∇Φ(x)3 .П РИМЕР 3. Если положить K = L и на различных итерацияхв качестве вектора vk циклически с повторением выбирать единичные орты e1 , e2 , .
. . , en , e1 , . . . то получится рассмотренный в §3.1 метод Гаусса-Зейделя4 .П РИМЕР 4. Выбор на k-й итерации vk = wk = ATek дает ABSкласс методов [1]. Чтобы для различных итераций выполнялосьусловие Ki ⊥ Kj , возможно либо хранение всех vk с их ортогонализацией по мере нахождения (схема Хуанга), либо пересчет матрицыA (схема Абрамова). Первый вариант ведет к увеличению расходов2Англ. Residual norm Steepest Descent — RnSD.Это утверждение можно получить и непосредственным анализом функционала,однако выкладки получаются значительно длиннее, чем в Примере 1.4Обратный порядок выбора соответствует обратному методу Гаусса-Зейделя.330памяти для реализации алгоритма, второй — к изменению заполненности матрицы.
Следовательно, данные методы пригодны лишьдля решения плотных и/или небольших систем; с другой стороны,их единственным условием сходимости является существование решения как такового, а потому данный класс пригоден для СЛАУ свырожденнными и неквадратными матрицами.Непосредственно из определения векторов v k следует, что в данных методах на k-й итерации поправка к решению вычисляется изусловия обращения k-го уравнения в тождество.В простейшем случае, если предполагать, что матрица A квадратная и невырожденная, вычислительная схема имеет следующий вид (приведен вариант с пересчетом матрицы, метод ABR1ORT):Выполнять для i = 1, nbi Tx := x +akai k22 iВыполнять для j = i + 1, n(ai , aj )α :=kai k22aj := aj − αaibj := bj − αbiувеличить jувеличить i4.3.Два важных выбора подпространствВ предыдущем параграфе были рассмотрены методы наискорейшегоспуска (см. пример 1), в котором подпространства K и L были связаны соотношением K = L, и наискорейшего уменьшения невязки (см.пример 2), основанный на соотношении L = AK.
Сами подпространства являлись одномерными, в качестве базиса K выступал векторневязки, и было показано, что в этих случаях задача проектирования эквивалентна задаче минимизации функционалов.Как оказывается, подобные утверждения справедливы и в гораздо более общих случаях, которые имеют важное значение при построении более сложных и эффективных методов.Теорема 4.3.1 Если матрица A симметрична и положительноопределена, то задача проектирования решения СЛАУ (1) на любоеподпространство K ортогонально к нему самому 5 является эквивалентной задаче минимизации функционала Φ 1 (x) = kx − x∗ k2A напространстве K.5Т.е.
ортогонально к пространству L = K.31Доказательство. По условию теоремы K = L, а следовательно, V == W. Выражение для функционала Φ1 было уже найдено (см. (4.8)),а сам функционал по свойствам нормы является строго выпуклым.Таким образом, сформулированная в условии задача минимизациисводится к нахождениюy = arg min Φ1 (x0 + Vy).yРассмотрим эту задачу. В силу выпуклости достаточно найти стационарную точку функционала Ψ(y) = Φ 1 (x0 + Vy), т.е. решить систему ∇Ψ(y) = 0. Пренебрегая постоянным слагаемым в (4.8), имеемΨ(y) = (x0 + Vy)TA(x0 + Vy) − 2bT(x0 + Vy) == (xT0 A − bT )x0 − bTx0 + 2(xT0 A − bT )Vy + y T(VTAV)y == y T (VT AV)y − r0Tx0 − bTx0 − 2r0T Vy.Градиент этого функционала равен∇Ψ(y) = 2VTAVy − 2VTr0 .Приравнивая его к нулю, получимy = (VTAV)−1 VTr0 ,что в точности совпадает с выражением для y из формулы (4.4), еслив ней положить V = W.2Теорема 4.3.2 Для произвольной невырожденной матрицы A задача проектирования решения СЛАУ (1) на любое подпространствоK ортогонально к подпространству L = AK является эквивалентной задаче минимизации функционала Φ 2 (x) = krx k22 на пространстве K.Доказательство.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
