Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности (1051191), страница 5
Текст из файла (страница 5)
Подставив в формулу (4.4) соотношение для базисов W = AV, получимy = (VTATAV)−1 VTATr0 .Это означает, что рассматриваемая ситуация эквивалентна выборуL = K для симметризованной системы ATAx = ATb. Учитывая соотношение (4) и применяя к такой системе предыдущую теорему, получаем сформулированное в условии утверждение.2324.4.Подпространства КрыловаПри построении и реализации проекционных методов важную рольиграют так называемые подпространства Крылова, часто выбираемые в качестве K.Определение 4.4.1 Подпространством Крылова размерности m,порожденным вектором v и матрицей A называется линейное пространствоdefnoKm (v, A) = span v, Av, A2 v, .
. . , Am−1 v .(4.10)В качестве вектора v обычно выбирается невязка начального приближения r0 ; тогда выбор подпространства L и способ построения базисов подпространств полностью определяет вычислительную схемуметода6 .К идее использования подпространств Крылова можно прийти,например, следующим образом.При построении релаксационных методов (см. §3.1) использовалось представление матрицы A в виде A = D − E − F. Было также показано, что методы Якоби и Гаусса-Зейделя являются частными случаями класса методов, основанного на расщеплении A в видеразности (3.9) двух матриц K и R.
Тогда исходная система (1) можетбыть записана в видеKx = b + Rx = b + (K − A)x ,что позволяет построить итерационный процессKxk+1 = Kxk + (b − Axk ) ,или, что то же самое,xk+1 = xk + K−1 rk .(4.11)Выберем K = I и R = I − A, тогда процесс (4.11) будет сведен квидуxk+1 = xk + rk ,(4.12)откуда следуетxk = x0 + r0 + r1 + . .
. + rk−1 .(4.13)Умножив обе части (4.12) слева на (−A) и прибавив к ним b, получимb − Axk+1 = b − Axk − Ark = rk − Ark ,6Иногда, впрочем, могут использоваться эквивалентные подходы, основанныена теоремах 4.3.1 и 4.3.2. См., напр., построение метода GMRES в §5.3.33что позволяет найти выражение для невязки на k-ой итерации черезневязку начального приближения:rk = (I − A)rk−1 = (I − A)k r0 .(4.14)После подстановки (4.14) в (4.13) получаемnxk = x 0 + k−1Xj=0o(I − A)j r0 .т.е. δx ∈ span r0 , Ar0 , .
. . , Ak−1 r0 = Kk (r0 , A).Непосредственно из определения вытекают следующие свойстваподпространств Крылова. Пусть q — полином, такой что q(A)v = 0,причем имеет степень deg q = µ, минимальную среди всех таких полиномов. Тогда• ∀m (m ≥ µ) : Km (v, A) = Kµ (v, A).Более того, Kµ инвариантно относительно A.• ∀m :dim(Km ) = m⇐⇒m ≤ µ.Кроме того, из (4.14) следует что в методах, использующих подпространства Крылова, невязка на k-ой итерации выражается черезначальную невязку некоторым матричным полиномом (см. приложение B).4.5.Базис подпространства Крылова.Ортогонализация АрнольдиДля построения базиса в пространстве Крылова K m (v1 , A) можноприменить следующий подход.Найдем сначала векторы w1 = v1 , w2 = Aw1 , w3 = A2 w1 == Aw2 ,.
. . , wm = Am−1 w1 = Awm−1 . По определению (4.10),Km (v1 , A) = span{w1 , w2 , . . . , wm }.Перейдем теперь от {w1 , w2 , . . . , wm } к {v1 , v2 , . . . , vm }, применивпроцедуру ортогонализацииvk+1 = wk+1 −kXαi vii=1и затем пронормировав полученные векторы.34(4.15)Предположим, что предыдущие k векторов уже построены, т.е.∀i, j (1 ≤ i, j ≤ k) :(vi , vj ) =(0, i 6= j1, i = j(4.16)Тогда (4.15) можно записать какvk+1 = Avk −kXαi vi .i=1Для выполнения условия ортогональности v k+1 ко всем предыдущимвекторам умножим это равенство скалярно на v j (j ≤ k) и приравняем результат к нулю(Avk , vj ) −kXαi (vi , vj ) = 0 .(4.17)i=1С учетом (4.16) отсюда легко получить выражение для коэффициентов αjαj = (Avk , vj ) .Описанный метод может быть оформлен в виде следующей процедуры, называемой ортогонализацией Арнольди 7 [4].Входные данные: v1 , такой что kv1 k2 = 1; A; mВыполнять для j = 1, mw := AvjВыполнять для i = 1, jhij := (w, vi )w := w − hij viувеличить ihj+1,j := kwk2Если hj+1,j = 0, то КОНЕЦ.
Базис построен.wvj+1 :=hj+1,jувеличить jЗамечание. Для коэффициентов ортогонализации здесь введенадвойная индексация, с учетом которой внутренний цикл алгоритмаалгебраически записывается какhj+1,j vj+1 = Avj −jXhij vi .(4.18)i=17Данный вариант построения базиса основан на модифицированной ортогонализации Грама-Шмидта, примененной к векторам v1 , .
. . , Am−1 v1 . Возможны и другиеварианты, например ортогонализация Арнольди-Хаусхолдера [13].35Коэффициенты ортогонализации hi,j можно объединить в видематрицы H, дополнив в ней недостающие позиции нулями. Приэтом, как видно из алгоритма, для заданной размерности пространства m генерируется m + 1 векторов. Последний вектор v m+1 (возможно, нулевой) в матричном представлении означает расширениебазиса V одним дополнительным столбцом, т.е.Vm = [v1 |v2 | . . . |vm ] ;Vm+1 = [v1 |v2 | . . . |vm |vm+1 ] .Соответствующий вектору vm+1 коэффициент hm+1,m означает расширение матрицы H одной дополнительной строкой (возможно, нулевой).Пусть Hm — это (m + 1) × m матрица коэффициентов h ij , дополненная последней строкой за счет hm+1,m , а Hm — та же самая матрица без последней строки, имеющая размерность m × m. Тогда непосредственно из описания алгоритма Арнольди и (4.18) следует чтоматрица Hm является матрицей в верхней форме Хессенберга и длянее справедливы соотношенияAVm = Vm+1 Hm = Vm Hm + wm eTm ;TVmAVm(4.19)(4.20)= Hm .Кроме того, вследствие ортонормальности базиса {v j } имеет место равенствоVTvk = ek .(4.21)П РИМЕР.
Пусть m = 2, а матрица A и вектор v1 равны1 1 0A = 0 1 1 ;0 0 10v1 = 1 .0Подпространство Крылова K2 (v1 , A) естьspan {v1 , Av1 } = {u | u = α1 (0, 1, 0)T + α2 (1, 1, 0)T , αi ∈ R} == {u | u = (β1 , β2 , 0)T , βi ∈ R}.Применение алгоритма Арнольди дает36v1 = (0, 1, 0)Tпри j = 1w := Av1 = (1, 1, 0)Tпри i = 1h11 := 1w := (1, 1, 0)T − (0, 1, 0)T = (1, 0, 0)Th21 := kwk2 = 1v2 := w = (1, 0, 0)Tпри j = 2w := Av2 = (1, 0, 0)Tпри i = 1h12 := 0w := (1, 0, 0)Tпри i = 2h22 := 1w := (1, 0, 0)T − (1, 0, 0)T = (0, 0, 0)Th32 := kwk2 = 0v3 := w = (0, 0, 0)TРезультатом является базис из векторов v 1 = (0, 1, 0)T и v2 == (1, 0, 0)T , а во введенных ранее матричных обозначениях:Vm+1Hm0 1 0= 1 0 0 ;0 0 01 0= 1 1 ;0 0Vm0 1= 1 0 ;0 0Hm =1 01 1!.Непосредственным вычислением нетрудно проверить для этихматриц ранее приведенные соотношения (4.19) и (4.20).4.6.Биортогонализация ЛанцошаАлгоритм ортогонализации Арнольди, предложенный в предыдущем параграфе, для построения каждого нового вектора v k требуетнахождения (k − 1) скалярных произведений и столько же операцийлинейного комбинирования.
Однако, как оказывается, при отказеот требования ортогональности базиса в пользу некоторого болееобщего условия можно построить процедуру меньшей сложности.mОпределение 4.6.1 Системы векторов {x i }mi=1 и {yi }i=1 называютсябиортогональными, если скалярное произведение (x i , yi ) обращается в ноль при i 6= j .Очевидно, что в случае xi ≡ yi условие биортогональности сводится к обычному условию ортогональности.Теорема 4.6.1 Пусть векторы v1 и w1 таковы, что (v1 , w1 ) 6= 0 иmпусть системы векторов {vi }mi=1 и {wi }i=1 определяются соотношениями:vi+1 = Avi − αi vi − βi vi−1 ,37v0 ≡ 0 ;(4.22)wi+1 = ATwi − αi wi − βi wi−1 ,w0 ≡ 0 ;(Avi , wi )αi =;(vi , wi )(vi , wi )βi =,β1 ≡ 0 .(vi−1 , wi−1 )(4.23)(4.24)(4.25)Тогдаm• системы {vi }mi=1 и {wi }i=1 являются биортогональными;m• каждая из систем {vi }mi=1 и {wi }i=1 является линейно независимой и образует базис в Km (v1 , A) и Km (w1 , AT ) соответственно.Доказательство. Первое утверждение теоремы доказывается по индукции.
Действительно, пара векторов v 1 и w1 удовлетворяет условию биортогонализации. Предположим теперь, что уже построеныбиортогональные наборы {v1 , . . . , vi } и {w1 , . . . , wi }, и далее покажем,что для вектора vi+1 , определяемого соотношением (4.22), имеет место (vi+1 , wk ) = 0, k = 1, i.Умножим (4.22) скалярно на wk(vi+1 , wk ) = (Avi , wk ) − αi (vi , wk ) − βi (vi−1 , wk ) .Если k = i, то по предположению индукции последнее скалярноепроизведение обращается в ноль и(vi+1 , wi ) = (Avi , wi ) − αi (vi , wi ) =(Avi , wi )= (Avi , wi ) −(vi , wi ) = 0 .(vi , wi )Если же k < i, то(vi+1 , wk ) = (vi , ATwk ) − βi (vi−1 , wk ) == (vi , wk+1 + αk wk + βk wk−1 ) − βi (vi−1 , wk ) == (vi , wk+1 ) + αk (vi , wk ) + βk (vi , wk−1 ) − βi (vi−1 , wk ) .По предположению индукции, при k < i − 1 все четыре скалярныепроизведения обращаются в ноль; при k = i − 1 равны нулю скалярные произведения во втором и третьем слагаемых, и тогда(vi+1 , wi−1 ) = (vi , wi ) − βi (vi−1 , wi−1 )(vi , wi )= (vi , wi ) −(vi−1 , wi−1 ) = 0 .(vi−1 , wi−1 )Аналогичным образом доказывается, что (v k , wi+1 ) = 0 для k == 1, i.38Чтобы доказать второе утверждение теоремы, заметим, что непосредственно из (4.22) следует span{v 1 , .
. . , vm } = Km (v1 , A). Остаетсялишь показать линейную независимость векторов v 1 , . . . , vm . Предположим от противного, что существуют коэффициенты γ k , для которыхγ1 v1 + γ 2 v2 + . . . + γ m vm = 0 .Составляя скалярные произведения с векторами w k , k = 1, m, получимγk (vk , wk ) = 0 ,k = 1, m ,а так как по ранее доказанной биортогональности (v k , wk ) 6= 08 , товсе коэффициенты γk должны быть нулевыми.Аналогичные рассуждения для w1 , . .
. , wm завершают доказательство теоремы.2Процедура построения векторов v и w согласно формулам (4.22)–(4.25) называется биортогонализацией Ланцоша.Очевидно, что (4.22) является частным случаем (4.18), где извсей суммы оставлены только два последних слагаемых; точно также (4.23) является частным случаем (4.18), записанной для матрицы AT и вектора w1 .