Единое окно доступа к образовательным ресурсам

Методы решения СЛАУ большой размерности: Учебное пособие

Голосов: 4

В данном учебном пособии рассматриваются форматы хранения разреженных матриц и реализация операций над ними (в частности, матрично-векторное умножение), предобусловливание (неполное LU-разложение), классические итерационные и проекционные методы на подпространствах Крылова. Приведены основные теоретические обоснования, рассмотрены системы с симметричными и несимметричными матрицами, приведены конкретные примеры. Для желающих использовать при реализации методов готовое программное обеспечение в приложении даны соответствующие Internet-адреса.

Приведенный ниже текст получен путем автоматического извлечения из оригинального PDF-документа и предназначен для предварительного просмотра.
Изображения (картинки, формулы, графики) отсутствуют.
    памяти для реализации алгоритма, второй — к изменению запол-
ненности матрицы. Следовательно, данные методы пригодны лишь
для решения плотных и/или небольших систем; с другой стороны,
их единственным условием сходимости является существование ре-
шения как такового, а потому данный класс пригоден для СЛАУ с
вырожденнными и неквадратными матрицами.
   Непосредственно из определения векторов v k следует, что в дан-
ных методах на k-й итерации поправка к решению вычисляется из
условия обращения k-го уравнения в тождество.
   В простейшем случае, если предполагать, что матрица A ква-
дратная и невырожденная, вычислительная схема имеет следую-
щий вид (приведен вариант с пересчетом матрицы, метод ABR1ORT):

                      Выполнять для i = 1, n
                                       bi T
                            x := x +        a
                                      ai 2 i
                                          2
                            Выполнять для j = i + 1, n
                                         (ai , aj )
                                   α :=
                                            ai 22
                                   aj := aj − αai
                                   bj := bj − αbi
                            увеличить j
                      увеличить i



4.3.       Два важных выбора подпространств
В предыдущем параграфе были рассмотрены методы наискорейшего
спуска (см. пример 1), в котором подпространства K и L были связа-
ны соотношением K = L, и наискорейшего уменьшения невязки (см.
пример 2), основанный на соотношении L = AK. Сами подпростран-
ства являлись одномерными, в качестве базиса K выступал вектор
невязки, и было показано, что в этих случаях задача проектирова-
ния эквивалентна задаче минимизации функционалов.
   Как оказывается, подобные утверждения справедливы и в гораз-
до более общих случаях, которые имеют важное значение при по-
строении более сложных и эффективных методов.

Теорема 4.3.1 Если матрица A симметрична и положительно
определена, то задача проектирования решения СЛАУ (1) на любое
подпространство K ортогонально к нему самому 5 является экви-
валентной задаче минимизации функционала Φ 1 (x) = x − x∗ 2 на
                                                          A
пространстве 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) =
         = (xTA − bT )x0 − bTx0 + 2(xTA − bT )Vy + y T(VTAV)y =
             0                       0
         = y T (VT AV)y − r0 x0 − bTx0 − 2r0 Vy.
                           T               T


   Градиент этого функционала равен

                      Ψ(y) = 2VTAVy − 2VTr0 .

Приравнивая его к нулю, получим

                        y = (VTAV)−1 VTr0 ,

что в точности совпадает с выражением для y из формулы (4.4), если
в ней положить V = W.                                           2


Теорема 4.3.2 Для произвольной невырожденной матрицы A зада-
ча проектирования решения СЛАУ (1) на любое подпространство
K ортогонально к подпространству L = AK является эквивалент-
ной задаче минимизации функционала Φ 2 (x) = rx 2 на простран-
                                                2
стве K.

Доказательство. Подставив в формулу (4.4) соотношение для бази-
сов W = AV, получим

                      y = (VTATAV)−1 VTATr0 .

Это означает, что рассматриваемая ситуация эквивалентна выбору
L = K для симметризованной системы ATAx = ATb. Учитывая соот-
ношение (4) и применяя к такой системе предыдущую теорему, по-
лучаем сформулированное в условии утверждение.               2




                                   32


4.4.    Подпространства Крылова
При построении и реализации проекционных методов важную роль
играют так называемые подпространства Крылова, часто выбирае-
мые в качестве K.

Определение 4.4.1 Подпространством Крылова размерности m,
порожденным вектором v и матрицей A называется линейное про-
странство
                          def
               Km (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) получаем
                                                          
                                          k−1
                          xk = x 0 +           (I − A)j  r0 .
                                          j=0


т.е. δ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 }, применив
процедуру ортогонализации
                                                    k
                             vk+1 = wk+1 −               αi vi            (4.15)
                                                   i=1

и затем пронормировав полученные векторы.

                                          34


      Предположим, что предыдущие k векторов уже построены, т.е.

                                                                  0, i = j
                ∀i, j (1 ≤ i, j ≤ k) :        (vi , vj ) =                   (4.16)
                                                                  1, i = j

      Тогда (4.15) можно записать как
                                                   k
                            vk+1 = Avk −                αi vi .
                                                  i=1

Для выполнения условия ортогональности v k+1 ко всем предыдущим
векторам умножим это равенство скалярно на v j (j ≤ k) и приравня-
ем результат к нулю
                                          k
                        (Avk , vj ) −          αi (vi , vj ) = 0 .           (4.17)
                                         i=1

   С учетом (4.16) отсюда легко получить выражение для коэффи-
циентов αj
                          αj = (Avk , vj ) .
   Описанный метод может быть оформлен в виде следующей про-
цедуры, называемой ортогонализацией Арнольди 7 [4].

            Входные данные: v1 , такой что v1 2 = 1; A; m
            Выполнять для j = 1, m
                  w := Avj
                  Выполнять для i = 1, j
                          hij := (w, vi )
                          w := w − hij vi
                  увеличить i
                  hj+1,j := w 2
                  Если hj+1,j = 0, то КОНЕЦ. Базис построен.
                             w
                  vj+1 :=
                           hj+1,j
            увеличить j

   Замечание. Для коэффициентов ортогонализации здесь введена
двойная индексация, с учетом которой внутренний цикл алгоритма
алгебраически записывается как
                                                        j
                         hj+1,j vj+1 = Avj −                 hij vi .        (4.18)
                                                       i=1
  7
   Данный вариант построения базиса основан на модифицированной ортогонали-
зации Грама-Шмидта, примененной к векторам 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 eT ;
                                             m                        (4.19)
               T
              Vm AVm       = Hm .                                     (4.20)


   Кроме того, вследствие ортонормальности базиса {v j } имеет ме-
сто равенство
                                VTvk = ek .                           (4.21)


   П РИМЕР. Пусть m = 2, а матрица A и вектор v1 равны
                                                       
                        1 1 0                       0
                                                   
                  A =  0 1 1 ;             v1 =  1  .
                        0 0 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}.


   Применение алгоритма Арнольди дает

                                     36


v1 = (0, 1, 0)T
при j = 1       w := Av1 = (1, 1, 0)T
        при i = 1       h11 := 1      w := (1, 1, 0)T − (0, 1, 0)T = (1, 0, 0)T
        h21 := w 2 = 1
        v2 := w = (1, 0, 0)T
при j = 2       w := Av2 = (1, 0, 0)T
        при i = 1       h12 := 0      w := (1, 0, 0)T
        при i = 2       h22 := 1      w := (1, 0, 0)T − (1, 0, 0)T = (0, 0, 0)T
        h32 := w 2 = 0
        v3 := w = (0, 0, 0)T

   Результатом является базис из векторов v 1 = (0, 1, 0)T и v2 =
= (1, 0, 0)T , а во введенных ранее матричных обозначениях:
                                                            
                          0 1 0                        0 1
                                                        
              Vm+1    =  1 0 0 ;            Vm   =  1 0 ;
                          0 0 0                        0 0
                                
                          1 0
                                                       1 0
                 Hm   =  1 1 ;              Hm =                 .
                                                         1 1
                          0 0
   Непосредственным вычислением нетрудно проверить для этих
матриц ранее приведенные соотношения (4.19) и (4.20).


4.6.   Биортогонализация Ланцоша
Алгоритм ортогонализации Арнольди, предложенный в предыду-
щем параграфе, для построения каждого нового вектора v k требует
нахождения (k − 1) скалярных произведений и столько же операций
линейного комбинирования. Однако, как оказывается, при отказе
от требования ортогональности базиса в пользу некоторого более
общего условия можно построить процедуру меньшей сложности.

Определение 4.6.1 Системы векторов {x i }m и {yi }m называются
                                         i=1      i=1
биортогональными, если скалярное произведение (x i , yi ) обращает-
ся в ноль при i = j .

   Очевидно, что в случае xi ≡ yi условие биортогональности сводит-
ся к обычному условию ортогональности.

Теорема 4.6.1 Пусть векторы v1 и w1 таковы, что (v1 , w1 ) = 0 и
пусть системы векторов {vi }m и {wi }m определяются соотноше-
                            i=1      i=1
ниями:
               vi+1 = Avi − αi vi − βi vi−1 ,       v0 ≡ 0 ;           (4.22)

                                     37


               wi+1 = ATwi − αi wi − βi wi−1 ,    w0 ≡ 0 ;                    (4.23)
                      (Avi , wi )
                 αi =              ;                                          (4.24)
                       (vi , wi )
                        (vi , wi )
                 βi =                ,   β1 ≡ 0 .                             (4.25)
                      (vi−1 , wi−1 )
Тогда
   • системы {vi }m и {wi }m являются биортогональными;
                  i=1      i=1

   • каждая из систем {vi }m и {wi }m является линейно независи-
                           i=1       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 ) = 08 , то
все коэффициенты γk должны быть нулевыми.
   Аналогичные рассуждения для w1 , . . . , wm завершают доказа-
тельство теоремы.                                                   2


   Процедура построения векторов v и w согласно формулам (4.22)–
(4.25) называется биортогонализацией Ланцоша.
   Очевидно, что (4.22) является частным случаем (4.18), где из
всей суммы оставлены только два последних слагаемых; точно так
же (4.23) является частным случаем (4.18), записанной для матри-
цы AT и вектора w1 . При этом коэффициенты ортогонализации в
(4.22) и (4.23) одинаковы.
   Из сказанного следует, что можно записать аналог матричной
формулы (4.20)
                            T
                           Wm AVm = Tm ,                   (4.26)
где Tm — симметричная9 трехдиагональная матрица, элементы ко-
торой определяются соотношениями

                                  [Tm ]ii = αi (vi , wi ) ;             (4.27)
                  [Tm ]i+1,i = [Tm ]i,i+1 = βi+1 (vi , wi ) ,           (4.28)

легко выводимыми из (4.22)–(4.23).
   Замечание. Основным недостатком биортогонализации Ланцо-
ша является возможность возникновения ситуации, когда (v i , wi ) =
= 0; при этом продолжение процесса становится невозможным из-за
неопределенности коэффициента βi+1 .



  8
    Иначе способ определения коэффициентов β сделал бы невозможным построе-
ние всех следующих векторов.
  9
    Симметричность следует из того, что для построения vi+1 и wi+1 используются
одни и те же коэффициенты ортогонализации.


                                      39


Глава 5

Методы крыловского типа

В §4.1 был построен класс проекционных методов, основой которо-
го служит условие Петрова-Галеркина (4.1) и формула (4.5). В §4.2
были приведены примеры построения нескольких методов для до-
статочно простого случая одномерных подпространств K и L.
    Однако выбор m = 1, как правило, дает медленную сходимость
и соответствующие методы оказываются пригодными лишь для
СЛАУ небольшой размерности или специальных случаев. В на-
стоящее время широкое распространение получили проекционные
методы, которые в качестве K используют описанные в §4.4 под-
пространства Крылова размерности больше единицы. Такие методы
называются методами крыловского типа 1 . Для упрощения вычис-
лений (WTAV)−1 из формулы (4.5) в них используются соотношения
(4.19)–(4.21).




5.1.       FOM: Метод полной ортогонализации
Наиболее характерным примером применения описанного в §4.1
подхода является метод полной ортогонализации, известный так-
же как метод Арнольди или FOM2 .
   Положим K = L и будем проектировать искомое решение на под-
пространство Крылова Km (v1 , A), где в качестве v1 выступает про-
нормированный вектор начальной невязки

                             v1 = r0 /β ,    β = r0   2   .   (5.1)

Тогда V = W, для построения базиса можно воспользоваться орто-
гонализацией Арнольди (§4.5), а матрица W TAV в силу (4.20) будет
  1
      Англ. Krylov sequence methods.
  2
      Full Orthogonalization Method.


                                            40



    
Яндекс цитирования Яндекс.Метрика