Método generalizado de residuos mínimos

En matemáticas , el método generalizado de residuos mínimos (GMRES) es un método iterativo para la solución numérica de un sistema de ecuaciones lineales no simétrico e indefinido . El método aproxima la solución mediante un vector en un subespacio de Krylov con residuo mínimo . La iteración de Arnoldi se utiliza para hallar este vector.

El método GMRES fue desarrollado por Yousef Saad y Martin H. Schultz en 1986. ^{[ 1 ]} Es una generalización y mejora del método MINRES debido a Paige y Saunders en 1975. ^{[ 2 ] [ 3 ]} El método MINRES requiere que la matriz sea simétrica , pero tiene la ventaja de que solo requiere el manejo de tres vectores. GMRES es un caso especial del método DIIS desarrollado por Peter Pulay en 1980. DIIS es aplicable a sistemas no lineales .

Denotemos la norma euclidiana de cualquier vector v por${\displaystyle \|v\|}$. Denotemos el sistema (cuadrado) de ecuaciones lineales que se va a resolver mediante$${\displaystyle Ax=b.}$$ Se supone que la matriz A es invertible y de tamaño m x m . Además, se supone que b está normalizada, es decir, que${\displaystyle \|b\|=1}$.

El n -ésimo subespacio de Krylov para este problema es$${\displaystyle K_{n}=K_{n}(A,r_{0})=\operatorname {span} \,\{r_{0},Ar_{0},A^{2}r_{0},\ldots ,A^{n-1}r_{0}\}.\,}$$ dónde${\displaystyle r_{0}=b-Ax_{0}}$es el residuo inicial dada una suposición inicial${\displaystyle x_{0}\neq 0}$. Claramente${\displaystyle r_{0}=b}$si${\displaystyle x_{0}=0}$.

GMRES se aproxima a la solución exacta de${\displaystyle Ax=b}$por el vector${\displaystyle x_{n}\in x_{0}+K_{n}}$que minimiza la norma euclidiana del residuo${\displaystyle r_{n}=b-Ax_{n}}$.

Los vectores${\displaystyle r_{0},Ar_{0},\ldots A^{n-1}r_{0}}$podría ser casi linealmente dependiente , por lo que en lugar de esta base, se utiliza la iteración de Arnoldi para encontrar vectores ortonormales.${\displaystyle q_{1},q_{2},\ldots ,q_{n}\,}$que constituyen la base para${\displaystyle K_{n}}$. En particular,${\displaystyle q_{1}=\|r_{0}\|_{2}^{-1}r_{0}}$.

Por lo tanto, el vector${\displaystyle x_{n}\in x_{0}+K_{n}}$se puede escribir como${\ Displaystyle x_ {n} = x_ {0} + Q_ {n} y_ {n}}$con${\displaystyle y_{n}\in \mathbb {R} ^{n}}$, dónde${\displaystyle Q_{n}}$es la matriz m por n formada por${\displaystyle q_{1},\ldots ,q_{n}}$. En otras palabras, encontrar la n -ésima aproximación de la solución (es decir,${\displaystyle x_{n}}$) se reduce a encontrar el vector${\displaystyle y_{n}}$, que se determina minimizando el residuo como se describe a continuación .

El proceso de Arnoldi también construye${\displaystyle {\tilde {H}}_{n}}$, un (${\displaystyle n+1}$)-por-${\displaystyle n}$matriz de Hessenberg superior que satisface$${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}\,}$$ una igualdad que se utiliza para simplificar el cálculo de${\displaystyle y_{n}}$(véase §  Resolución del problema de mínimos cuadrados ). Nótese que, para matrices simétricas, se obtiene una matriz tridiagonal simétrica, lo que da lugar al método MINRES .

Porque columnas de${\displaystyle Q_{n}}$son ortonormales, tenemos$${\displaystyle {\begin{aligned}\left\|r_{n}\right\|&=\left\|b-Ax_{n}\right\|\\&=\left\|b-A(x_{0}+Q_{n}y_{n})\right\|\\&=\left\|r_{0}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-Q_{n+1}{\tilde {H}}_{n}y_{n}\right\|\\&=\left\|Q_{n+1}(\beta e_{1}-{\tilde {H}}_{n}y_{n})\right\|\\&=\left\|\beta e_{1}-{\tilde {H}}_{n}y_{n}\right\|\end{aligned}}}$$ dónde$${\displaystyle e_{1}=(1,0,0,\ldots ,0)^{T}\,}$$es el primer vector en la base estándar de${\displaystyle \mathbb {R} ^{n+1}}$, y$${\displaystyle \beta =\|r_{0}\|\,,}$$${\displaystyle r_{0}}$siendo el vector residual del primer ensayo (generalmente${\displaystyle b}$). Por eso,${\displaystyle x_{n}}$se puede encontrar minimizando la norma euclidiana del residuo$${\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.}$$ Este es un problema de mínimos cuadrados lineales de tamaño n .

Esto da como resultado el método GMRES. En el${\displaystyle n}$-ª iteración:

1. calcular${\displaystyle q_{n}}$con el método Arnoldi;
2. encontrar el${\displaystyle y_{n}}$lo cual minimiza${\displaystyle \|r_{n}\|}$;
3. calcular${\displaystyle x_{n}=x_{0}+Q_{n}y_{n}}$;
4. Repita el proceso si el residuo aún no es lo suficientemente pequeño.

En cada iteración, un producto matriz-vector${\displaystyle Aq_{n}}$debe calcularse. Esto cuesta aproximadamente${\displaystyle 2m^{2}}$operaciones de punto flotante para matrices densas generales de tamaño${\displaystyle m}$, pero el costo puede disminuir a${\displaystyle O(m)}$para matrices dispersas . Además del producto matriz-vector,${\displaystyle O(nm)}$Las operaciones de punto flotante deben calcularse en la n -ésima iteración.

La enésima iteración minimiza el residuo en el subespacio de Krylov.${\displaystyle K_{n}}$Dado que cada subespacio está contenido en el siguiente, el residuo no aumenta. Tras m iteraciones, donde m es el tamaño de la matriz A , el espacio de Krylov K _m es todo R ^m y, por lo tanto, el método GMRES alcanza la solución exacta. Sin embargo, la idea es que, tras un número reducido de iteraciones (en relación con m ), el vector x _n ya constituye una buena aproximación a la solución exacta.

Esto no ocurre en general. De hecho, un teorema de Greenbaum, Pták y Strakoš establece que para toda secuencia no creciente a ₁ , ..., a _{m − 1} , a _m = 0, se puede encontrar una matriz A tal que ‖ r _n ‖ = a _n para todo n , donde r _n es el residuo definido anteriormente. En particular, es posible encontrar una matriz para la cual el residuo permanece constante durante m − 1 iteraciones, y solo se reduce a cero en la última iteración.  

En la práctica, sin embargo, GMRES suele funcionar bien. Esto se puede demostrar en situaciones específicas. Si la parte simétrica de A , es decir${\displaystyle (A^{T}+A)/2}$, es definida positiva , entonces $${\displaystyle \|r_{n}\|\leq \left(1-{\frac {\lambda _{\min }^{2}(1/2(A^{T}+A))}{\lambda _{\max }(A^{T}A)}}\right)^{n/2}\|r_{0}\|,}$$ dónde${\displaystyle \lambda _{\mathrm {min} }(M)}$y${\displaystyle \lambda _{\mathrm {max} }(M)}$denotan el valor propio más pequeño y el más grande de la matriz.${\displaystyle M}$, respectivamente. ^{[ 4 ]}

Si A es simétrica y definida positiva, entonces incluso tenemos $${\displaystyle \|r_{n}\|\leq \left({\frac {\kappa _{2}(A)^{2}-1}{\kappa _{2}(A)^{2}}}\right)^{n/2}\|r_{0}\|.}$$ dónde${\displaystyle \kappa _{2}(A)}$denota el número de condición de A en la norma euclidiana.

En el caso general, donde A no es definida positiva, tenemos $${\displaystyle {\frac {\|r_{n}\|}{\|b\|}}\leq \inf _{p\in P_{n}}\|p(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|,\,}$$ donde P _n denota el conjunto de polinomios de grado como máximo n con p (0) = 1, V es la matriz que aparece en la descomposición espectral de A , y σ ( A ) es el espectro de A. En términos generales, esto significa que la convergencia rápida ocurre cuando los autovalores de A se agrupan lejos del origen y A no está demasiado lejos de la normalidad . ^{[ 5 ]}

Todas estas desigualdades limitan únicamente los residuos en lugar del error real, es decir, la distancia entre la iteración actual x _n y la solución exacta.

Al igual que otros métodos iterativos, GMRES se suele combinar con un método de preacondicionamiento para acelerar la convergencia.

El coste de las iteraciones crece como O( n² ), donde n es el número de iteración. Por lo tanto, a veces el método se reinicia después de un número, digamos ^k , de iteraciones, con x _k como estimación inicial. El método resultante se denomina GMRES( k ) o GMRES reiniciado. Para matrices no definidas positivas, este método puede sufrir estancamiento en la convergencia, ya que el subespacio reiniciado suele estar cerca del subespacio anterior.

Las deficiencias de GMRES y GMRES reiniciado se abordan mediante el reciclaje del subespacio de Krylov en los métodos de tipo GCRO, como GCROT y GCRODR. ^{[ 6 ]} El reciclaje de subespacios de Krylov en GMRES también puede acelerar la convergencia cuando se necesitan resolver secuencias de sistemas lineales. ^{[ 7 ]}

La iteración de Arnoldi se reduce a la iteración de Lanczos para matrices simétricas. El método del subespacio de Krylov correspondiente es el método de residuos mínimos (MinRes) de Paige y Saunders. A diferencia del caso asimétrico, el método MinRes se define mediante una relación de recurrencia de tres términos . Se puede demostrar que no existe un método del subespacio de Krylov para matrices generales que, a pesar de tener una relación de recurrencia corta, minimice las normas de los residuos, como lo hace GMRES.

Otro tipo de métodos se basa en la iteración asimétrica de Lanczos , en particular el método BiCG . Estos métodos utilizan una relación de recurrencia de tres términos, pero no alcanzan el residuo mínimo, por lo que el residuo no disminuye monótonamente. Ni siquiera se garantiza la convergencia.

La tercera clase la conforman métodos como CGS y BiCGSTAB . Estos también trabajan con una relación de recurrencia de tres términos (por lo tanto, sin optimalidad) e incluso pueden terminar prematuramente sin alcanzar la convergencia. La idea detrás de estos métodos es elegir adecuadamente los polinomios generadores de la secuencia de iteración.

Ninguna de estas tres clases es la mejor para todas las matrices; siempre hay ejemplos en los que una clase supera a la otra. Por lo tanto, en la práctica se prueban varios algoritmos de resolución para determinar cuál es el más adecuado para un problema dado.

Una parte del método GMRES consiste en encontrar el vector${\displaystyle y_{n}}$lo cual minimiza $${\displaystyle \left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|.}$$ Tenga en cuenta que${\displaystyle {\tilde {H}}_{n}}$es una matriz de ( n  +  1) por n , por lo tanto da un sistema lineal sobredeterminado de n + 1 ecuaciones para n incógnitas.

El mínimo se puede calcular utilizando una descomposición QR : encontrar una matriz ortogonal Ω _n de ( n  +  1) × ( n  +  1) y una matriz triangular superior de ( n + 1) × n.  ${\displaystyle {\tilde {R}}_{n}}$de tal manera que $${\displaystyle \Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.}$$ La matriz triangular tiene una fila más que columnas, por lo que su fila inferior consta de ceros. Por lo tanto, se puede descomponer como: $${\displaystyle {\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},}$$ dónde${\displaystyle R_{n}}$es una matriz triangular de n por n (por lo tanto, cuadrada).

La descomposición QR se puede actualizar de forma económica de una iteración a la siguiente, porque las matrices de Hessenberg difieren solo en una fila de ceros y una columna: $${\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},}$$ donde h _n+1 = ( h _{1, n +1} , ..., h _{n +1, n +1} ) ^T . Esto implica que premultiplicar la matriz de Hessenberg por Ω _n , aumentada con ceros y una fila con identidad multiplicativa, produce casi una matriz triangular: $${\displaystyle {\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&\rho \\0&\sigma \end{bmatrix}}}$$ Esto sería triangular si σ es cero. Para remediar esto, se necesita la rotación de Givens.$${\displaystyle G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}}$$ dónde $${\displaystyle c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\text{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.}$$ Con esta rotación de Givens, formamos $${\displaystyle \Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.}$$ En efecto, $${\displaystyle \Omega _{n+1}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&r_{n+1,n+1}\\0&0\end{bmatrix}}}$$es una matriz triangular con${\textstyle r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2}}}}$.

Dada la descomposición QR, el problema de minimización se resuelve fácilmente observando que $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|\Omega _{n}({\tilde {H}}_{n}y_{n}-\beta e_{1})\right\|\\&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|.\end{aligned}}}$$ Denotando el vector${\displaystyle \beta \Omega _{n}e_{1}}$por $${\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}$$ con g _n ∈ R ⁿ y γ _n ∈ R , esto es $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|\\&=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y_{n}-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|.\end{aligned}}}$$ El vector y que minimiza esta expresión viene dado por $${\displaystyle y_{n}=R_{n}^{-1}g_{n}.}$$ Nuevamente, los vectores${\displaystyle g_{n}}$son fáciles de actualizar. ^{[ 8 ]}

• Método del gradiente biconjugado