Método estabilizado por gradiente biconjugado

En las siguientes secciones, ( x , y ) = x ^T y denota el producto escalar de vectores. Para resolver un sistema lineal Ax = b , BiCGSTAB comienza con una estimación inicial x ₀ y procede de la siguiente manera:

1. r ₀ = b − Ax ₀
2. Elija un vector arbitrario r̂ ₀ tal que ( r̂ ₀ , r ₀ ) ≠ 0 , por ejemplo, r̂ ₀ = r ₀
3. ρ ₀ = ( r̂ ₀ , r ₀ )
4. p ₀ = r ₀
5. Para i = 1, 2, 3, …
   1. v = Ap _{i −1}
   2. α = ρ _{i −1} /( r̂ ₀ , v )
   3. h = x _{i −1} + α p _{i −1}
   4. s = r _{i −1} − α v
   5. Si h es suficientemente preciso, es decir, si s es suficientemente pequeño, entonces establezca x _i = h y salga.
   6. t = Como
   7. ω = ( t , s )/( t , t )
   8. x _i = h + ω s
   9. r _i = s − ω t
   10. Si x _i es suficientemente preciso, es decir, si r _i es suficientemente pequeño, entonces salir
   11. ρ _i = ( r̂ ₀ , r _i )
   12. β = ( ρ _yo / ρ _{yo −1} )( α / ω )
   13. p _i = r _i + β ( p _{i −1} − ω v )

En algunos casos, elegir el vector r̂ ₀ aleatoriamente mejora la estabilidad numérica. ^{[ 1 ]}

Los precondicionadores se utilizan habitualmente para acelerar la convergencia de los métodos iterativos. Para resolver un sistema lineal Ax = b con un precondicionador K = K ₁ K ₂ ≈ A , el método BiCGSTAB precondicionado comienza con una estimación inicial x ₀ y procede de la siguiente manera:

1. r ₀ = b − Ax ₀
2. Elija un vector arbitrario r̂ ₀ tal que ( r̂ ₀ , r ₀ ) ≠ 0 , por ejemplo, r̂ ₀ = r ₀
3. ρ ₀ = ( r̂ ₀ , r ₀ )
4. p ₀ = r ₀
5. Para i = 1, 2, 3, …
   1. y = K −1 2   K −1 1   p _{i −1}
   2. v = Ay
   3. α = ρ _{i −1} /( r̂ ₀ , v )
   4. h = x _{i −1} + α y
   5. s = r _{i −1} − α v
   6. Si h es suficientemente preciso, entonces x _i = h y salir.
   7. z = K −1 2   K −1 1   s
   8. t = Az
   9. ω = ( K −1 1   t , K −1 1   s )/( K −1 1   t , K −1 1   t )
   10. x _i = h + ω z
   11. r _i = s − ω t
   12. Si x _i es suficientemente preciso, entonces salir
   13. ρ _i = ( r̂ ₀ , r _i )
   14. β = ( ρ _yo / ρ _{yo −1} )( α / ω )
   15. p _i = r _i + β ( p _{i −1} − ω v )

Esta formulación es equivalente a aplicar BiCGSTAB sin preacondicionamiento al sistema explícitamente preacondicionado.

Ãx̃ = b̃

con à = K −1 1   A K −1 2   , x̃ = K ₂ x y b̃ = K −1 1   b . En otras palabras, con esta formulación son posibles tanto el precondicionamiento izquierdo como el derecho.

En BiCG, las direcciones de búsqueda p _i y p̂ _i y los residuos r _i y r̂ _i se actualizan utilizando las siguientes relaciones de recurrencia :

p _i = r _{i −1} + β _i p _{i −1} ,

p̂ _yo = r̂ _{yo −1} + β _yo p̂ _{yo −1} ,

r _i = r _{i −1} − α _i Ap _i ,

r̂ _yo = r̂ _{yo −1} − α _yo A ^T p̂ _yo .

Las constantes α _i y β _i se eligen de la siguiente manera:

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

β _i = ρ _i / ρ _{i −1}

donde ρ _i = ( r̂ _{i −1} , r _{i −1} ) de modo que los residuos y las direcciones de búsqueda satisfacen la biorthogonalidad y la biconjugación, respectivamente, es decir, para i ≠ j ,

( r̂ _i , r _j ) = 0 ,

( p̂ _i , Ap _j ) = 0 .

Es sencillo demostrar que

r _i = P _i ( A ) r ₀ ,

r̂ _i = P _i ( A ^T ) r̂ ₀ ,

p _{i +1} = T _i ( A ) r ₀ ,

p̂ _{i +1} = T _i ( A ^T ) r̂ ₀

donde P _i ( A ) y T _i ( A ) son polinomios de grado i en A . Estos polinomios satisfacen las siguientes relaciones de recurrencia:

P _i ( A ) = P _{i −1} ( A ) − α _i A T _{i −1} ( A ) ,

T _i ( A ) = P _i ( A ) + β _{i +1} T _{i −1} ( A ) .

No es necesario realizar un seguimiento explícito de los residuos y las direcciones de búsqueda de BiCG. En otras palabras, las iteraciones de BiCG se pueden realizar implícitamente. En BiCGSTAB, se desea tener relaciones de recurrencia para

r̃ _i = Q _i ( A ) P _i ( A ) r ₀

donde Q _i ( A ) = ( I − ω ₁ A )( I − ω ₂ A )⋯( I − ω _i A ) con constantes adecuadas ω _j en lugar de r _i = P _i ( A ) r ₀ con la esperanza de que Q _i ( A ) permita una convergencia más rápida y suave en r̃ _i que r _i .

De las relaciones de recurrencia para P _i ( A ) y T _i ( A ) y la definición de Q _i ( A ) se deduce que

Q _i ( A ) P _i ( A ) r ₀ = ( I − ω _i A )( Q _{i −1} ( A ) P _{i −1} ( A ) r ₀ − α _i A Q _{i −1} ( A ) T _{i −1} ( A ) r ₀ ) ,

lo cual implica la necesidad de una relación de recurrencia para Q _i ( A ) T _i ( A ) r ₀ . Esto también se puede derivar de las relaciones BiCG:

Q _i ( A ) T _i ( A ) r ₀ = Q _i ( A ) P _i ( A ) r ₀ + β _{i +1} ( I − ω _i A ) Q _{i −1} ( A ) T _{i −1} ( A ) r ₀ .

De forma similar a la definición de r̃ _i , BiCGSTAB define

p̃ _{i +1} = Q _i ( A ) T _i ( A ) r ₀ .

Escritas en forma vectorial, las relaciones de recurrencia para p̃ _i y r̃ _i son:

p̃ _i = r̃ _{i −1} + β _i ( I − ω _{i −1} A ) p̃ _{i −1} ,

r̃ _i = ( I − ω _i A )( r̃ _{i −1} − α _i A p̃ _i ).

Para derivar una relación de recurrencia para x _i , definimos

s _i = r̃ _{i −1} − α _i A p̃ _i .

La relación de recurrencia para r̃ _i se puede escribir entonces como

r̃ _i = r̃ _{i −1} − α _i A p̃ _i − ω _i As _i ,

que corresponde a

x _i = x _{i −1} + α _i p̃ _i + ω _i s _i .

Ahora queda determinar las constantes BiCG α _i y β _i y elegir un ω _i adecuado .

En BiCG, β _i = ρ _i / ρ _{i −1} con

ρ _i = ( r̂ _{i −1} , r _{i −1} ) = ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) .

Dado que BiCGSTAB no realiza un seguimiento explícito de r̂ _i o r _i , ρ _i no se puede calcular inmediatamente a partir de esta fórmula. Sin embargo, se puede relacionar con el escalar.

ρ̃ _i = ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) = ( r̂ ₀ , Q _{i −1} ( A ) P _{i −1} ( A ) r ₀ ) = ( r̂ ₀ , r _{i −1} ) .

Debido a la biorthogonalidad, r _{i −1} = P _{i −1} ( A ) r ₀ es ortogonal a U _{i −2} ( A ^T ) r̂ ₀ donde U _{i −2} ( A ^T ) es cualquier polinomio de grado i − 2 en A ^T . Por lo tanto, solo los términos de orden más alto de P _{i −1} ( A ^T ) y Q _{i −1} ( A ^T ) importan en los productos escalares ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) y ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) . Los coeficientes principales de P _{i −1} ( A ^T ) y Q _{i −1} ( A ^T ) son (−1) ^{i −1} α ₁ α ₂ ⋯ α _{i −1} y (−1) ^{i −1} ω ₁ ω ₂ ⋯ ω _{i −1} , respectivamente. De ello se deduce que

ρ _yo = ( α ₁ / ω ₁ )( α ₂ / ω ₂ )⋯( α _{yo −1} / ω _{yo −1} ) ρ̃ _yo ,

y por lo tanto

β _yo = ρ _yo / ρ _{yo −1} = ( ρ̃ _yo / ρ̃ _{yo −1} )( α _{yo −1} / ω _{yo −1} ) .

De manera similar, se puede derivar una fórmula simple para α _{i . En BiCG,}

α _i = ρ _i /( p̂ _i , Ap _i ) = ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ )/( T _{i −1} ( A ^T ) r̂ ₀ , A T _{i −1} ( A ) r ₀ ) .

De forma similar al caso anterior, solo los términos de orden superior de P _{i −1} ( A ^T ) y T _{i −1} ( A ^T ) importan en los productos escalares gracias a la biorthogonalidad y la biconjugación. Sucede que P _{i −1} ( A ^T ) y T _{i −1} ( A ^T ) tienen el mismo coeficiente principal. Por lo tanto, pueden ser reemplazados simultáneamente con Q _{i −1} ( A ^T ) en la fórmula, lo que lleva a

α _i = ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ )/( Q _{i −1} ( A ^T ) r̂ ₀ , A T _{i −1} ( A ) r ₀ ) = ρ̃ _i /( r̂ ₀ , A Q _{i −1} ( A ) T _{i −1} ( A ) r ₀ ) = ρ̃ _i /( r̂ ₀ , Ap̃ _i ) .

Finalmente, BiCGSTAB selecciona ω _i para minimizar r̃ _i = ( I − ω _i A ) s _i en norma 2 como función de ω _i . Esto se logra cuando

(( I − ω _i A ) s _i , As _i ) = 0 ,

dando el valor óptimo

ω _i = ( As _i , s _i )/( As _i , As _i ) .