En matemáticas, el método residual mínimo generalizado (GMRES) es un método iterativo para la solución numérica de un sistema no simétrico indefinido de ecuaciones lineales . El método aproxima la solución por el vector en un subespacio de Krylov con un residuo mínimo . Se utiliza la iteración de Arnoldi para encontrar este vector.
El método GMRES fue desarrollado por Yousef Saad y Martin H. Schultz en 1986. [1] GMRES es una generalización y mejora del método MINRES , mientras que este último también requiere que la matriz sea simétrica, pero tiene la ventaja de que solo requiere manipulación. de tres vectores. [2] [3] Fue desarrollado por Chris Paige y Michael Saunders en 1975. GMRES también es un caso especial del método DIIS desarrollado por Peter Pulay en 1980. DIIS también es aplicable a sistemas no lineales.
El método
Denote la norma euclidiana de cualquier vector v por. Denote 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á normalizado, es decir, que.
El n -ésimo subespacio de Krylov para este problema es
dónde es el error inicial dada una suposición inicial . Claramente Si .
GMRES se aproxima a la solución exacta de por el vector que minimiza la norma euclidiana del residuo .
Los vectores podría estar cerca de ser linealmente dependiente , por lo que en lugar de esta base, se usa la iteración de Arnoldi para encontrar vectores ortonormales que forman una base para . En particular,.
Por lo tanto, el vector Se puede escribir como con , dónde es la matriz m -por- n formada por.
El proceso de Arnoldi también produce un ()-por-matriz de Hessenberg superior con
Para matrices simétricas, en realidad se logra una matriz tri-diagonal simétrica, lo que resulta en el método minres .
Porque las columnas de son ortonormales, tenemos
dónde
es el primer vector en la base estándar de, y
siendo el primer vector de prueba (normalmente cero). Por eso, se puede encontrar minimizando la norma euclidiana del residuo
Este es un problema de mínimos cuadrados lineales de tamaño n .
Esto produce el método GMRES. Sobre el-th iteración:
- calcular con el método Arnoldi;
- encuentra el que minimiza ;
- calcular ;
- repita si el residuo aún no es lo suficientemente pequeño.
En cada iteración, un producto matriz-vector debe calcularse. 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,Las operaciones de punto flotante deben calcularse en la n -ésima iteración.
Convergencia
El n º iterate minimiza el residual en el subespacio Krylov. Dado que cada subespacio está contenido en el siguiente subespacio, el residual 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 un pequeño número de iteraciones (relativas am ), el vector x n ya es una buena aproximación a la solución exacta.
Esto no sucede en general. De hecho, un teorema de Greenbaum, Pták y Strakoš establece que para cada secuencia no creciente a 1 ,…, 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 residual permanece constante durante m - 1 iteraciones, y solo cae a cero en la última iteración.
En la práctica, sin embargo, GMRES a menudo funciona bien. Esto se puede demostrar en situaciones específicas. Si la parte simétrica de A , es, es positivo definido , entonces
dónde y denotar el valor propio más pequeño y más grande de la matriz, respectivamente. [4]
Si A es simétrica y definida positiva, incluso tenemos
dónde denota el número de condición de A en la norma euclidiana.
En el caso general, donde A no es definida positiva, tenemos
donde P n denota el conjunto de polinomios de grado a lo sumo 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]
Todas estas desigualdades limitan solo 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 a medida que O ( n 2 ), 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 un estancamiento en la convergencia, ya que el subespacio reiniciado suele estar cerca del subespacio anterior.
Las deficiencias de GMRES y GMRES reiniciadas se abordan mediante el reciclaje del subespacio de Krylov en los métodos de tipo GCRO como GCROT y GCRODR. [6] El reciclaje de los subespacios de Krylov en GMRES también puede acelerar la convergencia cuando es necesario 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 del subespacio de Krylov correspondiente es el método residual mínimo (MinRes) de Paige y Saunders. A diferencia del caso asimétrico, el método MinRes viene dado por una relación de recurrencia de tres términos . Se puede demostrar que no existe un método subespacial de Krylov para matrices generales, que viene dado por una relación de recurrencia corta y, sin embargo, minimiza 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 mínimo residual y, por lo tanto, el residual no disminuye monótonamente 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 funcionan 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, en la práctica se prueban varios solucionadores 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
Tenga en cuenta que es una matriz ( n + 1) -por- n , por lo que da un sistema lineal sobrerestringido de n +1 ecuaciones para n incógnitas.
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) -by- 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
dónde es una matriz triangular n- por- n (por 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 por una fila de ceros y una columna:
donde h n + 1 = ( h 1, n + 1 , ..., h n + 1, n + 1 ) T . Esto implica que premultiplicando la matriz de Hessenberg con Ω n , aumentada con ceros y una fila con identidad multiplicativa, produce casi una matriz triangular:
Esto sería triangular si σ es cero. Para remediar esto, se necesita la rotación de Givens.
dónde
Con esta rotación de Givens, formamos
En efecto,
es una matriz triangular.
Dada la descomposición QR, el problema de minimización se resuelve fácilmente al señalar que
Denotando el vector por
con g n ∈ R n y γ n ∈ R , esto es
El vector y que minimiza esta expresión viene dado por
De nuevo, los vectores son fáciles de actualizar. [8]
Código de ejemplo
GMRES normal (octava MATLAB / GNU)
función [x, e] = gmres ( A, b, x, max_iterations, umbral ) n = longitud ( A ); m = max_iteraciones ; % usa x como vector inicial r = b - A * x ; b_norm = norma ( b ); error = norma ( r ) / b_norm ; % inicializar 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 = norma ( r ); Q (:, 1 ) = r / r_norm ; beta = r_norm * e1 ; % Nota: este no es el escalar beta en la sección "El método" anterior, sino el escalar beta multiplicado por e1 para k = 1 : m % ejecutar arnoldi [ H ( 1 : k + 1 , k ) Q (:, k + 1 )] = arnoldi ( A , Q , k ); % eliminar el último elemento en H ith fila y actualizar la matriz de rotación [ H ( 1 : k + 1 , k ) cs ( k ) sn ( k )] = aplicar_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % actualiza el vector residual beta ( k + 1 ) = - sn ( k ) * beta ( k ); beta ( k ) = cs ( k ) * beta ( k ); error = abs ( beta ( k + 1 )) / b_norm ; % guarda el error e = [ e ; error ]; si ( error <= umbral ) romper ; final final % si no se alcanza el umbral, k = m en este punto (y no m + 1) % calcula el resultado y = H ( 1 : k , 1 : k ) \ beta ( 1 : k ); x = x + Q (:, 1 : k ) * y ; final% ------------------------------------------------- ---%% Función Arnoldi%% ------------------------------------------------- ---%función [h, q] = arnoldi ( A, Q, k ) q = A * Q (:, k ); % Krylov Vector para i = 1 : k % Gram-Schmidt modificado, manteniendo la matriz de Hessenberg h ( i ) = q ' * Q (:, i ); q = q - h ( i ) * Q (:, i ); final h ( k + 1 ) = norma ( q ); q = q / h ( k + 1 ); final% ------------------------------------------------- --------------------%% Aplicando la rotación de Givens a H col%% ------------------------------------------------- --------------------%función [h, cs_k, sn_k] = apply_givens_rotation ( h, cs, sn, k ) % aplica 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 ; final % actualiza los siguientes valores sin cos 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 ; final%% ---- Calcular la matriz de rotación dada ---- %%función [cs, sn] = givens_rotation ( v1, v2 ) % si (v1 == 0)% cs = 0;% sn = 1;% demás t = raíz cuadrada ( 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 ; % finalfinal
Ver también
Referencias
- ^ Y. Saad y MH Schultz
- ^ Paige y Saunders, "Solución de sistemas dispersos indefinidos de ecuaciones lineales", SIAM J. Numer. Anal., Vol 12, página 617 (1975) https://doi.org/10.1137/0712047
- ^ N.Nifa. "Tesis Doctoral" (PDF) .
- ↑ Eisenstat, Elman y Schultz, Thm 3.3. Nota: todos los resultados de GCR también son válidos para GMRES, cf. Saad y Schultz
- ^ Trefethen y Bau, Thm 35.2
- ^ Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). "Reciclaje de subespacios de Krylov para aplicaciones CFD y un nuevo solucionador de reciclaje híbrido". Revista de Física Computacional . 303 : 222. arXiv : 1501.03358 . Código bibliográfico : 2015JCoPh.303..222A . doi : 10.1016 / j.jcp.2015.09.040 .
- ^ Galia, André (2014). Reciclaje de métodos subespaciales de Krylov para secuencias de sistemas lineales (Ph.D.). TU Berlín. doi : 10.14279 / depositonce-4147 .
- ^ Stoer y Bulirsch, §8.7.2
Notas
- A. Meister, Numerik linearer Gleichungssysteme , 2da edición, Vieweg 2005, ISBN 978-3-528-13135-7 .
- Y. Saad, Métodos iterativos para sistemas lineales dispersos , segunda edición, Sociedad de Matemáticas Industriales y Aplicadas , 2003. ISBN 978-0-89871-534-7 .
- Y. Saad y MH Schultz, "GMRES: Un algoritmo residual mínimo generalizado para resolver sistemas lineales no simétricos", SIAM J. Sci. Stat. Computación. , 7 : 856-869, 1986. doi : 10.1137 / 0907058 .
- SC Eisenstat, HC Elman y MH Schultz, "Métodos iterativos variacionales para sistemas no simétricos de ecuaciones lineales", SIAM Journal on Numerical Analysis , 20 (2), 345–357, 1983.
- J. Stoer y R. Bulirsch, Introducción al análisis numérico , tercera edición, Springer, Nueva York, 2002. ISBN 978-0-387-95452-3 .
- Lloyd N. Trefethen y David Bau, III, Álgebra lineal numérica , Sociedad de Matemáticas Industriales y Aplicadas, 1997. ISBN 978-0-89871-361-9 .
- Dongarra y col. , Plantillas para la solución de sistemas lineales: bloques de construcción para métodos iterativos , 2da edición, SIAM, Filadelfia, 1994
- Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). "Reciclaje de subespacios de Krylov para aplicaciones CFD y un nuevo solucionador de reciclaje híbrido". Journal of Computational Physics 303: 222. doi: 10.1016 / j.jcp.2015.09.040