Método generalizado de residuos mínimos

Método para la solución numérica de ciertos sistemas de ecuaciones

En matemáticas, el método de residuos mínimos generalizados (GMRES) es un método iterativo para la solución numérica de un sistema indefinido no simétrico de ecuaciones lineales . El método aproxima la solución por el vector en un subespacio de Krylov con residuos mínimos . La iteración de Arnoldi se utiliza para encontrar 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.

El método

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

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

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

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

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

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

Como las columnas de son ortonormales, tenemos donde es el primer vector en la base estándar de , y es el primer vector residual de prueba (usualmente ). Por lo tanto, se puede encontrar al minimizar la norma euclidiana del residual. Este es un problema de mínimos cuadrados lineal de tamaño n . ${\displaystyle Q_{n}}$ $${\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}}}$$ $${\displaystyle e_{1}=(1,0,0,\ldots ,0)^{T}\,}$$ ${\displaystyle \mathbb {R} ^{n+1}}$ $${\displaystyle \beta =\|r_{0}\|\,,}$$ ${\displaystyle r_{0}}$ ${\displaystyle b}$ ${\displaystyle x_{n}}$ $${\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.}$$

Esto genera el método GMRES. En la iteración -ésima: ${\displaystyle n}$

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

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

Convergencia

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

En general, esto no sucede. De hecho, un teorema de Greenbaum, Pták y Strakoš establece que para cada 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 se mantenga constante durante m  − 1 iteraciones, y solo caiga 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 , es definida positiva , entonces donde y denotan el valor propio más pequeño y más grande de la matriz , respectivamente. ^[4] ${\displaystyle (A^{T}+A)/2}$ $${\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}\|,}$$ ${\displaystyle \lambda _{\mathrm {min} }(M)}$ ${\displaystyle \lambda _{\mathrm {max} }(M)}$ ${\displaystyle M}$

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

En el caso general, donde A no es definida positiva, tenemos 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 dice que la convergencia rápida ocurre cuando los valores propios de A se agrupan lejos del origen y A no está demasiado lejos de la normalidad . ^[5] $${\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 )|,\,}$$

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.

Extensiones del método

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

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

Las deficiencias de GMRES y GMRES reiniciado se solucionan 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]

Comparación con otros solucionadores

La iteración de Arnoldi se reduce a la iteración de Lanczos para matrices simétricas. El método de 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 está dado por una relación de recurrencia de tres términos . Se puede demostrar que no existe un método de subespacio de Krylov para matrices generales, que esté dado por una relación de recurrencia corta y, sin embargo, minimice las normas de los residuos, como lo hace GMRES.

Otra clase de métodos se basa en la iteración asimétrica de Lanczos, en particular el método BiCG . Estos utilizan una relación de recurrencia de tres términos, pero no alcanzan el residuo mínimo y, por lo tanto, el residuo no disminuye de manera monótona para estos métodos. La convergencia ni siquiera está garantizada.

La tercera clase está formada por 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 lograr 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, se prueban varios solucionadores en la práctica para ver cuál es el mejor para un problema determinado.

Resolver el problema de mínimos cuadrados

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

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

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: 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: Esta sería triangular si σ es cero. Para remediar esto, se necesita la rotación de Givens donde Con esta rotación de Givens, formamos De hecho, es una matriz triangular con . $${\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},}$$ $${\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}}}$$ $${\displaystyle G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}}$$ $${\displaystyle c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\text{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.}$$ $${\displaystyle \Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.}$$ $${\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}}}$$ ${\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 Denotando el vector por con g _n ∈ R ⁿ y γ _n ∈ R , esto es El vector y que minimiza esta expresión está dado por Nuevamente, los vectores son fáciles de actualizar. ^[8] $${\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}}}$$ ${\displaystyle \beta \Omega _{n}e_{1}}$ $${\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}$$ $${\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}}}$$ $${\displaystyle y_{n}=R_{n}^{-1}g_{n}.}$$ ${\displaystyle g_{n}}$

Código de ejemplo

GMRES regular (MATLAB / GNU Octave)

función [x, e] = gmres ( A, b, x, máx_iteraciones, umbral ) n = longitud ( A ); m = máx_iteraciones ;        % utiliza x como vector inicial r = b - A * x ;        b_norm = norma ( b ); error = norma ( r ) / b_norm ;        % inicializa los vectores 1D sn = ceros ( m , 1 ); cs = ceros ( m , 1 ); %e1 = ceros(n, 1); e1 = ceros ( m + 1 , 1 ); e1 ( 1 ) = 1 ; e = [ error ]; r_norm = norm ( r ); Q (:, 1 ) = r / r_norm ; % Nota: este no es el escalar beta de la sección "El método" anterior sino % el escalar beta multiplicado por e1 beta = r_norm * e1 ; para k = 1 : m                                       % ejecutar arnoldi [ H ( 1 : k + 1 , k ), Q (:, k + 1 )] = arnoldi ( A , Q , k ); % eliminar el último elemento en la fila H i y actualizar la matriz de rotación [ H ( 1 : k + 1 , k ), cs ( k ), sn ( k )] = apply_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % actualizar el vector residual beta ( k + 1 ) = - sn ( k ) * beta ( k ); beta ( k ) = cs ( k ) * beta ( k ); error = abs ( beta ( k + 1 )) / b_norm ;                                         % guardar el error e = [ e ; error ];     si ( error <= umbral ) break ; fin fin % si no se alcanza el umbral, k = m en este punto (y no m+1) % calcular el resultado y = H ( 1 : k , 1 : k ) \ beta ( 1 : k ); x = x + Q (:, 1 : k ) * y ; fin                       %----------------------------------------------------% % Función de Arnoldi % %----------------------------------------------------% función [h, q] = arnoldi ( A, Q, k ) q = A * Q (:, k ); % Vector de Krylov para i = 1 : k % Gram-Schmidt modificado, manteniendo la matriz de Hessenberg h ( i ) = q ' * Q (:, i ); q = q - h ( i ) * Q (:, i ); fin h ( k + 1 ) = norm ( q ); q = q / h ( k + 1 ); fin                                     %---------------------------------------------------------------------% % Aplicando rotación de dados a la columna H % %---------------------------------------------------------------------% function [h, cs_k, sn_k] = apply_givens_rotation ( h, cs, sn, k ) % aplicar para la i-ésima columna para i = 1 : k - 1 temp = cs ( i ) * h ( i ) + sn ( i ) * h ( i + 1 ); h ( i + 1 ) = - sn ( i ) * h ( i ) + cs ( i ) * h ( i + 1 ); h ( i ) = temp ; fin                                 % actualiza los próximos valores de seno coseno para la rotación [ cs_k , sn_k ] = givens_rotation ( h ( k ), h ( k + 1 ));        % eliminar H(i + 1, i) h ( k ) = cs_k * h ( k ) + sn_k * h ( k + 1 ); h ( k + 1 ) = 0.0 ; fin                %%----Calcular la matriz de rotación de Givens----%% función [cs, sn] = givens_rotation ( v1, v2 ) % si (v1 == 0) % cs = 0; % sn = 1; % de lo contrario t = sqrt ( v1 ^ 2 + v2 ^ 2 ); % cs = abs(v1) / t; % sn = cs * v2 / v1; cs = v1 / t ; % ver http://www.netlib.org/eispack/comqr.f sn = v2 / t ; % fin fin                 

Véase también

• Método del gradiente biconjugado

Referencias

1. Saad, Youcef; Schultz, Martin H. (1986). "GMRES: Un algoritmo generalizado de residuos mínimos para resolver sistemas lineales no simétricos". Revista SIAM de computación científica y estadística . 7 (3): 856–869. doi :10.1137/0907058. ISSN  0196-5204.
2. Paige y Saunders, "Solución de sistemas indefinidos dispersos de ecuaciones lineales", SIAM J. Numer. Anal., vol. 12, página 617 (1975) https://doi.org/10.1137/0712047
3. Nifa, Naoufal (2017). Solveurs performants pour l'optimisation sous contraintes en identificación de paramètres [ Solucionadores eficientes para optimización restringida en problemas de identificación de parámetros ] (Tesis) (en francés).
4. Eisenstat, Elman y Schultz 1983, Thm 3.3. Nota: todos los resultados para GCR también son válidos para GMRES, cf. Saad y Schultz 1986
5. Trefethen, Lloyd N.; Bau, David, III. (1997). Álgebra lineal numérica . Filadelfia: Sociedad de Matemáticas Industriales y Aplicadas. Teorema 35.2. ISBN 978-0-89871-361-9.{{cite book}}: CS1 maint: multiple names: authors list (link)
6. Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). "Reciclaje de subespacios de Krylov para aplicaciones de CFD y un nuevo solucionador de reciclaje híbrido". Journal of Computational Physics . 303 : 222. arXiv : 1501.03358 . Bibcode :2015JCoPh.303..222A. doi :10.1016/j.jcp.2015.09.040. S2CID  2933274.
7. Gaul, André (2014). Métodos de reciclaje del subespacio de Krylov para secuencias de sistemas lineales (Ph.D.). TU Berlin. doi :10.14279/depositonce-4147.
8. Stoer, Josef; Bulirsch, Roland (2002). Introducción al análisis numérico . Textos de matemáticas aplicadas (3.ª ed.). Nueva York: Springer. §8.7.2. ISBN 978-0-387-95452-3.

• Maestro, Andreas; Vömel, Christof (2005). Numerik linearer Gleichungssysteme . Wiesbaden: Vieweg. ISBN 978-3-528-13135-7.
• Saad, Y. (2003). Métodos iterativos para sistemas lineales dispersos (2.ª ed.). Filadelfia: SIAM. ISBN 978-0-89871-534-7.
• Eisenstat, Stanley C.; Elman, Howard C.; Schultz, Martin H. (1983). "Métodos iterativos variacionales para sistemas no simétricos de ecuaciones lineales". Revista SIAM de análisis numérico . 20 (2): 345–357. doi :10.1137/0720023. ISSN  0036-1429.
• Dongarra et al., Plantillas para la solución de sistemas lineales: bloques de construcción para métodos iterativos, 2.ª edición, SIAM, Filadelfia, 1994
• Imankulov, Timur; Lebedev, Danil; Matkerim, Bazargul; Daribayev, Beimbet; Kassymbek, Nurislam (8 de octubre de 2021). "Simulación numérica del flujo multifásico y multicomponente en medios porosos: análisis de la eficiencia del método basado en Newton". Fluidos . 6 (10): 355. doi : 10.3390/fluids6100355 . ISSN  2311-5521.