Metoda stabilizowana gradientem dwukoniugatu

W numerycznej algebrze liniowej metoda stabilizowanego gradientu biconjugate , często określana skrótem BiCGSTAB , jest iteracyjną metodą opracowaną przez HA van der Vorsta do numerycznego rozwiązywania niesymetrycznych układów liniowych . Jest to wariant metody gradientu dwukoniugatu (BiCG) i ma szybszą i płynniejszą zbieżność niż oryginalny BiCG, a także inne warianty, takie jak metoda kwadratu gradientu sprzężonego (CGS). Jest to podprzestrzeni Kryłowa . W przeciwieństwie do oryginalnej metody BiCG, nie wymaga ona mnożenia przez transpozycję macierzy systemu.

Kroki algorytmiczne

Niekondycjonowany BiCGSTAB

Aby rozwiązać układ liniowy Ax = b , BiCGSTAB zaczyna od wstępnego przypuszczenia x ₀ i postępuje w następujący sposób:

1. ₀ r = b − Ax ₀
2. Wybierz dowolny wektor r̂ ₀ taki, że ₀₀ ( r̂ , r ) ≠ 0 , np. ₀ r̂ = r ₀ . ( x , y ) oznacza iloczyn skalarny wektorów ( x , y ) = x ^T y
3. ₀₀ ρ = α = ω = 1
4. ₀₀ v = p = 0
5. Dla i = 1, 2, 3, …
   1. ₀ ρ _ja = ( r̂ , r _{ja −1} )
   2. β = ( ρ _ja / ρ _{ja -1} )( α / ω _{ja -1} )
   3. p _ja = r _{ja -1} + β ( p _{ja -1} - ω _{ja -1} v _{ja -1} )
   4. v _ja = Ap _ja
   5. ₀ α = ρ _ja /( r̂ , v _ja )
   6. h = x _{ja −1} + α p _ja
   7. Jeśli h jest wystarczająco dokładne, ustaw x _i = h i wyjdź
   8. s = r _{ja -1} - α v _ja
   9. t = Jak
   10. ω _ja = ( t , s )/( t , t )
   11. x _ja = h + ω _ja s
   12. Jeśli x _i jest wystarczająco dokładne, to zakończ
   13. r _ja = s - ω _ja t

Wstępnie kondycjonowany BiCGSTAB

Warunki wstępne są zwykle używane do przyspieszenia zbieżności metod iteracyjnych. Aby rozwiązać układ liniowy Ax = b z warunkiem wstępnym K = K ₁ K ₂ ≈ A , BiCGSTAB zaczyna od wstępnego przypuszczenia x ₀ i postępuje w następujący sposób:

1. ₀ r = b − Ax ₀
2. Wybierz dowolny wektor r̂ ₀ taki, że ₀₀ ( r̂ , r ) ≠ 0 , np. ₀ r̂ = r ₀
3. ₀₀ ρ = α = ω = 1
4. ₀₀ v = p = 0
5. Dla i = 1, 2, 3, …
   1. ₀ ρ _ja = ( r̂ , r _{ja −1} )
   2. β = ( ρ _ja / ρ _{ja -1} )( α / ω _{ja -1} )
   3. p _ja = r _{ja -1} + β ( p _{ja -1} - ω _{ja -1} v _{ja -1} )
   4. y =   K.
      -1 2    K.
      - 1 1   p _ja
   5. v _ja = tak
   6. ₀ α = ρ _ja /( r̂ , v _ja )
   7. h = x _{ja -1} + α y
   8. Jeśli h jest wystarczająco dokładne, to x _i = h i wyjdź
   9. s = r _{ja -1} - α v _ja
   10. z =   K.
       −1 2    K.
       −1 1   s
   11. t = Az
   12. ω _ja = (   K.
       -1 1   t ,   K.
       -1 1   s )/(   K.
       -1 1   t ,   K.
       -1 1   t )
   13. x _ja = h + ω _ja z
   14. Jeśli x _i jest wystarczająco dokładne, zakończ
   15. r _ja = s - ω _ja t

To sformułowanie jest równoważne z zastosowaniem BiCGSTAB bez wstępnego kondycjonowania do jawnie wstępnie kondycjonowanego systemu

Ãx̃ = b̃

gdzie à =   K.
-1 1   ZA   K.
-1 2   , x̃ = K. ₂ x i b̃ =   K.
-1 1   b . Innymi słowy, przy tym sformułowaniu możliwe jest zarówno warunkowanie wstępne w lewo, jak iw prawo.

Pochodzenie

BiCG w postaci wielomianu

i _W BiCG kierunki wyszukiwania pi _i p̂ są _i oraz reszty ri _i r̂ aktualizowane przy użyciu następujących relacji powtarzalności :

p _ja = r _{ja -1} + β _ja p _{ja -1} ,

p̂ _ja = r̂ _{ja -1} + β _ja p̂ _{ja -1} ,

r _ja = r _{ja -1} - α _ja Ap _ja ,

r̂ _ja = r̂ _{ja -1} − α _ja ZA ^T p̂ _ja .

Wybrano stałe α _i oraz β _i

α _ja = ρ _ja / ( p̂ _ja , Ap _ja ) ,

β _ja = ρ _ja / ρ _{ja −1}

gdzie ρ _i = ( r̂ _{i −1} , r _{i −1} ) tak, że reszty i kierunki poszukiwań spełniają odpowiednio biortogonalność i bikonjugację, tj. dla i ≠ j ,

( r̂ _ja , r _jot ) = 0 ,

( p̂ _ja , Ap _jo ) = 0 .

Pokazanie tego jest proste

r _ja = P _ja ( ZA ) r ₀ ,

r̂ _ja = P _ja ( ZA ^T ) r̂ ₀ ,

p _{ja +1} = T _ja ( ZA ) r ₀ ,

p̂ _{ja +1} = T _ja ( ZA ^T ) r̂ ₀

gdzie P _i ( A ) i T _i ( A ) są wielomianami i -tego stopnia w A . Te wielomiany spełniają następujące relacje powtarzalności:

P. _ja ( ZA ) = P. _{ja -1} ( ZA ) - α _ja ZA T _{ja -1} ( ZA ) ,

T _ja ( ZA ) = P _ja ( ZA ) + β _{ja +1} T _{ja -1} ( ZA ) .

Wyprowadzenie BiCGSTAB z BiCG

Nie jest konieczne wyraźne śledzenie pozostałości i kierunków poszukiwań BiCG. Innymi słowy, iteracje BiCG mogą być wykonywane niejawnie. W BiCGSTAB chce się mieć relacje powtarzalności dla

r̃ _ja = Q _ja ( ZA ) P _ja ( ZA ) r ₀

gdzie Q _ja ( A ) = ( I − ω ₁ A )( I − ω ₂ A )⋯( I − ω _i A ) z odpowiednimi stałymi ω _j zamiast r _i = P _ja ( A ) r ₀ w nadziei, że Q _i ( A ) umożliwi szybszą i płynniejszą zbieżność w r̃ _i niż r _i .

Z relacji rekurencji dla P _i ( A ) i T _i ( A ) oraz definicji Q _i ( A ) wynika, że

₀₀₀ Q _ja ( ZA ) P _ja ( ZA ) r = ( ja - ω _ja ZA ) ( Q _{ja -1} ( ZA ) P _{ja -1} ( ZA ) r - α _ja ZA Q _{ja -1} ( ZA ) T _{ja -1} ( ZA ) r ) ,

co pociąga za sobą konieczność relacji rekurencyjnej dla Q _i ( A ) T _i ( A ) r ₀ . Można to również wywnioskować z relacji BiCG:

₀₀ Q _ja ( ZA ) T _ja ( ZA ) r = Q _ja ( ZA ) P _ja ( ZA ) r + β _{ja +1} ( ja - ω _ja ZA ) Q _{ja -1} ( ZA ) P _{ja -1} ( ZA ) r ₀ .

Podobnie jak w przypadku definiowania r̃ _i , BiCGSTAB definiuje

p̃ _{ja +1} = Q _ja ( ZA ) T _ja ( ZA ) r ₀ .

Zapisane w postaci wektorowej relacje powtarzalności dla p̃ _i oraz r̃ _i są

p̃ _ja = r̃ _{ja -1} + β _ja ( ja - ω _{ja -1} ZA ) p̃ _{ja -1} ,

r̃ _ja = ( ja - ω _ja ZA ) ( r̃ _{ja -1} - α _ja ZA p̃ _ja ) .

Aby wyprowadzić relację rekurencyjną dla x _i , zdefiniuj

s _ja = r̃ _{ja -1} - α _ja ZA p̃ _ja .

Relację powtarzalności dla r̃ _i można zatem zapisać jako

r̃ _ja = r̃ _{ja -1} - α _ja ZA p̃ _ja - ω _ja Jako _ja ,

co odpowiada

x _ja = x _{ja -1} + α _ja p̃ _ja + ω _ja s _ja .

Wyznaczanie stałych BiCGSTAB

Teraz pozostaje wyznaczenie stałych BiCG α _i oraz β _i oraz wybranie odpowiedniego ω _i .

W BiCG, β _ja = ρ _ja / ρ _{ja −1} z

₀₀ ρ _ja = ( r̂ _{ja -1} , r _{ja -1} ) = ( P _{ja -1} ( ZA ^T ) r̂ , P _{ja -1} ( ZA ) r ) .

Ponieważ BiCGSTAB nie śledzi jawnie r̂ _i lub r _i , ρ _i nie jest natychmiast obliczalne z tego wzoru. Można to jednak powiązać ze skalarem

₀₀₀₀₀ ρ̃ _ja = ( Q _{ja -1} ( ZA ^T ) r̂ , P _{ja -1} ( ZA ) r ) = ( r̂ , Q _{ja -1} ( ZA ) P _{ja -1} ( ZA ) r ) = ( r̂ , r _{ja - 1} ) .

Ze względu na biortogonalność, r _{i −1} = P _{i −1} ( A ) r ₀ jest prostopadła do U _{i −2} ( A ^T ) r̂ ₀ gdzie U _{i −2} ( A ^T ) jest dowolnym wielomianem stopnia i − 2 w A ^T . Stąd _{i −1} ( A ⁾ r liczą się w iloczynach skalarnych ( P _{i −1} ( A ^T ) r̂ , P P _{i −1} ( A ^T ) ₀₀ Q _{i −1} ( A T ) ) tylko wyrazy najwyższego rzędu z i i ₀₀ ( Q _{ja -1} ( ZA ^T ) r̂ , P _{ja -1} ( ZA ) r ) . Wiodące współczynniki P _{i −1} ( A ^T ) i Q _{i −1} ( A ^T ) to (−1) ^{i −1} α ₁ α ₂ ⋯ α _{i −1} oraz (−1) ^{i −1} ω ₁ ω _{2 odpowiednio} ⋯ ω _{ja −1} . Wynika, że

ρ _ja = ( α ₁ / ω ₁ ) ( α ₂ / ω ₂ )⋯( α _{ja -1} / ω _{ja -1} ) ρ̃ _ja ,

a zatem

β _ja = ρ _ja / ρ _{ja -1} = ( ρ̃ _ja / ρ̃ _{ja -1} ) ( α _{ja -1} / ω _{ja -1} ) .

W podobny sposób można wyprowadzić prosty wzór na α _{i .} w BiCG,

₀₀₀₀ α _ja = ρ _ja /( p̂ _ja , Ap _ja ) = ( P _{ja -1} ( ZA ^T ) r̂ , P _{ja -1} ( ZA ) r )/( T _{ja -1} ( ZA ^T ) r̂ , ZA T _{ja - 1} ( A ) r ) .

Podobnie jak w powyższym przypadku, tylko wyrazy najwyższego rzędu P _{i −1} ( A ^T ) i T _{i −1} ( A ^T ) mają znaczenie w iloczynach skalarnych dzięki biortogonalności i bikonjugacji. Zdarza się, że P _{i −1} ( A ^T ) i T _{i −1} ( A ^T ) mają ten sam współczynnik wiodący. W ten sposób można je zastąpić jednocześnie Q _{i −1} ( A ^T ) we wzorze, który prowadzi do

_0000000 α _ja = ( Q _{ja -1} ( ZA ^T ) r̂ , P _{ja -1} ( ZA ) r )/( Q _{ja -1} ( ZA ^T ) r̂ , ZA T _{ja -1} ( ZA ) r ) = ρ̃ _ja / ( r̂ , ZA Q _{ja -1} ( ZA ) T _{ja -1} ( ZA ) r ) = ρ̃ _ja / ( r̂ , Ap̃ _ja ) .

Na koniec BiCGSTAB wybiera ω _i , aby zminimalizować r̃ _i = ( ja - ω _i A ) s _i w normie 2 jako funkcję ω _i . Osiąga się to, gdy

(( ja - ω _ja ZA ) s _ja , Ponieważ _ja ) = 0 ,

dając optymalną wartość

ω _ja = ( Jak _ja , s _ja )/( Jak _ja , Jak _ja ) .

Uogólnienie

BiCGSTAB można postrzegać jako połączenie BiCG i GMRES , gdzie po każdym kroku BiCG następuje krok GMRES( 1 ) (tj. . Jednak ze względu na zastosowanie minimalnych wielomianów resztkowych pierwszego stopnia taka naprawa może nie być skuteczna, jeśli macierz A ma duże zespolone pary własne. W takich przypadkach BiCGSTAB prawdopodobnie znajdzie się w stagnacji, co potwierdzają eksperymenty numeryczne.

Można oczekiwać, że minimalne wielomiany resztkowe wyższego stopnia mogą lepiej poradzić sobie z tą sytuacją. Daje to początek algorytmom obejmującym BiCGSTAB2 i bardziej ogólny BiCGSTAB( l ). W BiCGSTAB( l ) krok GMRES( l ) następuje po każdych l krokach BiCG. BiCGSTAB2 jest równoważne BiCGSTAB( l ) z l = 2 .

Zobacz też

• Metoda gradientu dwukoniugatowego
• Metoda kwadratu gradientu sprzężonego
• Metoda gradientu sprzężonego

• Van der Vorst, HA (1992). „Bi-CGSTAB: szybki i płynnie zbieżny wariant Bi-CG do rozwiązywania niesymetrycznych systemów liniowych”. SIAM J. Sci. Stan. Oblicz. 13 (2): 631–644. doi : 10.1137/0913035 . hdl : 10338.dmlcz/104566 .
•   Saad, Y. (2003). "§7.4.2 BICGSTAB" . Metody iteracyjne dla rzadkich systemów liniowych (wyd. 2). SYJAM. s. 231–234 . ISBN 978-0-89871-534-7 .
• ^ Gutknecht, MH (1993). „Warianty BICGSTAB dla macierzy o złożonym widmie”. SIAM J. Sci. Oblicz. 14 (5): 1020–1033. doi : 10.1137/0914062 .
• Bibliografia   _ Fokkema, DR (listopad 1993). „BiCGstab ( l ) dla równań liniowych obejmujących macierze niesymetryczne o widmie zespolonym” (PDF) . Transakcje elektroniczne w analizie numerycznej . Kent, OH: Kent State University. 1 : 11–32. ISSN 1068-9613 .