Método de residuo mínimo generalizado - Generalized minimal residual method

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 . 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. Es una generalización y mejora del método MINRES debido a Paige y Saunders en 1975. 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

Denote la norma euclidiana de cualquier vector v por . Denote el sistema (cuadrado) de ecuaciones lineales a resolver por ${\ Displaystyle \ | v \ |}$

${\ Displaystyle Ax = b. \,}$

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 . ${\ 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} \}. \,}$

donde es el error inicial dada una suposición inicial . Claramente si . ${\ 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 residual . ${\ Displaystyle Ax = b}$${\ Displaystyle x_ {n} \ in K_ {n}}$ ${\ Displaystyle r_ {n} = b-Ax_ {n}}$

Los vectores pueden ser casi dependientes linealmente , 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 ,. ${\ 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 . ${\ Displaystyle x_ {n} \ in 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}}$

El proceso de Arnoldi también produce una matriz ( ) -por- superior de Hessenberg con ${\ Displaystyle n + 1}$${\ Displaystyle n}$ ${\ Displaystyle {\ tilde {H}} _ {n}}$

${\ Displaystyle AQ_ {n} = Q_ {n + 1} {\ tilde {H}} _ {n}. \,}$

Para matrices simétricas, en realidad se logra una matriz tri-diagonal simétrica, lo que resulta en el método minres .

Debido a que las columnas de son ortonormales, tenemos ${\ Displaystyle Q_ {n}}$

${\ Displaystyle \ | r_ {n} \ | = \ | b-Ax_ {n} \ | = \ | bA (x_ {0} + Q_ {n} y_ {n}) \ | = \ | r_ {0} -AQ_ {n} y_ {n} \ | = \ | \ beta q_ {1} -AQ_ {n} y_ {n} \ | = \ | \ beta q_ {1} -Q_ {n + 1} {\ tilde {H}} _ {n} y_ {n} \ | = \ | Q_ {n + 1} (\ beta e_ {1} - {\ tilde {H}} _ {n} y_ {n}) \ | = \ | \ beta e_ {1} - {\ tilde {H}} _ {n} y_ {n} \ |, \,}$

dónde

${\ Displaystyle e_ {1} = (1,0,0, \ ldots, 0) ^ {T} \,}$

es el primer vector en la base estándar de , y ${\ Displaystyle \ mathbb {R} ^ {n + 1}}$

${\ Displaystyle \ beta = \ | r_ {0} \ | \ ,,}$

${\ Displaystyle x_ {0}}$siendo el primer vector de prueba (normalmente cero). Por tanto, se puede encontrar minimizando la norma euclidiana del residuo ${\ Displaystyle x_ {n}}$

${\ 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 produce el método GMRES. En la -ésima iteración: ${\ Displaystyle n}$

1. calcular con el método Arnoldi;${\ Displaystyle q_ {n}}$
2. encontrar 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 sobre las operaciones de punto flotante para matrices densas de tamaño en general , 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. ${\ Displaystyle Aq_ {n}}$${\ Displaystyle 2m ^ {2}}$ ${\ Displaystyle m}$${\ Displaystyle O (m)}$${\ Displaystyle O (nm)}$

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. ${\ Displaystyle K_ {n}}$

Esto no sucede en general. 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 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 decir , es definida positiva , entonces ${\ 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} \ |,}$

donde y denotan el valor propio más pequeño y más grande de la matriz , respectivamente. ${\ Displaystyle \ lambda _ {\ mathrm {min}} (M)}$${\ Displaystyle \ lambda _ {\ mathrm {max}} (M)}$${\ Displaystyle M}$

Si A es simétrica y definida positiva, incluso tenemos

${\ Displaystyle \ | r_ {n} \ | \ leq \ left ({\ frac {\ kappa _ {2} (A) ^ {2} -1} {\ kappa _ {2} (A) ^ {2} }} \ right) ^ {n / 2} \ | r_ {0} \ |.}$

donde denota el número de condición de A en la norma euclidiana. ${\ Displaystyle \ kappa _ {2} (A)}$

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 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 .

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 como O ( n ² ), donde n es el número de iteración. 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. El reciclaje de subespacios de Krylov en GMRES también puede acelerar la convergencia cuando es necesario resolver secuencias de sistemas lineales.

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 ${\ Displaystyle y_ {n}}$

${\ Displaystyle \ | {\ tilde {H}} _ {n} y_ {n} - \ beta e_ {1} \ |. \,}$

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. ${\ Displaystyle {\ tilde {H}} _ {n}}$

El mínimo se puede calcular usando una descomposición QR : encuentre una ( n  + 1) -por- ( n  + 1) matriz ortogonal Ω _n y una ( n  + 1) -por- n matriz triangular superior tal que ${\ Displaystyle {\ tilde {R}} _ {n}}$

${\ 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 cero. Por lo tanto, se puede descomponer como

${\ displaystyle {\ tilde {R}} _ {n} = {\ begin {bmatrix} R_ {n} \\ 0 \ end {bmatrix}},}$

donde es una matriz triangular n- por- n (por tanto, cuadrada). ${\ 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 por 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 premultiplicando la matriz de Hessenberg con Ω _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 {\ mbox {y}} \ 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}} \ quad {\ text {with}} \ quad r_ {n + 1, n + 1} = {\ sqrt {\ rho ^ {2} + \ sigma ^ {2 }}}}$

es una matriz triangular.

Dada la descomposición QR, el problema de minimización se resuelve fácilmente al señalar que

${\ Displaystyle \ | {\ tilde {H}} _ {n} y_ {n} - \ beta e_ {1} \ | = \ | \ Omega _ {n} ({\ tilde {H}} _ {n} y_ {n} - \ beta e_ {1}) \ | = \ | {\ tilde {R}} _ {n} y_ {n} - \ beta \ Omega _ {n} e_ {1} \ |.}$

Denotando el vector por ${\ Displaystyle \ beta \ Omega _ {n} e_ {1}}$

${\ Displaystyle {\ tilde {g}} _ {n} = {\ begin {bmatrix} g_ {n} \\\ gamma _ {n} \ end {bmatrix}}}$

con g _n ∈ R ⁿ y γ _n ∈ R , esto es

${\ Displaystyle \ | {\ tilde {H}} _ {n} y_ {n} - \ beta e_ {1} \ | = \ | {\ tilde {R}} _ {n} y_ {n} - \ beta \ Omega _ {n} e_ {1} \ | = \ left \ | {\ begin {bmatrix} R_ {n} \\ 0 \ end {bmatrix}} y_ {n} - {\ begin {bmatrix} g_ {n } \\\ gamma _ {n} \ end {bmatrix}} \ right \ |.}$

El vector y que minimiza esta expresión viene dado por

${\ Displaystyle y_ {n} = R_ {n} ^ {- 1} g_ {n}.}$

Nuevamente, los vectores son fáciles de actualizar. ${\ Displaystyle g_ {n}}$

Código de ejemplo

GMRES normal (octava MATLAB / GNU)

function [x, e] = gmres( A, b, x, max_iterations, threshold)
  n = length(A);
  m = max_iterations;

  % use x as the initial vector
  r = b - A * x;

  b_norm = norm(b);
  error = norm(r) / b_norm;

  % initialize the 1D vectors
  sn = zeros(m, 1);
  cs = zeros(m, 1);
  %e1 = zeros(n, 1);
  e1 = zeros(m+1, 1);
  e1(1) = 1;
  e = [error];
  r_norm = norm(r);
  Q(:,1) = r / r_norm;
  beta = r_norm * e1;     %Note: this is not the beta scalar in section "The method" above but the beta scalar multiplied by e1
  for k = 1:m

    % run arnoldi
    [H(1:k+1, k) Q(:, k+1)] = arnoldi(A, Q, k);
    
    % eliminate the last element in H ith row and update the rotation matrix
    [H(1:k+1, k) cs(k) sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);
    
    % update the residual vector
    beta(k + 1) = -sn(k) * beta(k);
    beta(k)     = cs(k) * beta(k);
    error       = abs(beta(k + 1)) / b_norm;

    % save the error
    e = [e; error];

    if (error <= threshold)
      break;
    end
  end
  % if threshold is not reached, k = m at this point (and not m+1) 
  
  % calculate the result
  y = H(1:k, 1:k) \ beta(1:k);
  x = x + Q(:, 1:k) * y;
end

%----------------------------------------------------%
%                  Arnoldi Function                  %
%----------------------------------------------------%
function [h, q] = arnoldi(A, Q, k)
  q = A*Q(:,k);   % Krylov Vector
  for i = 1:k     % Modified Gram-Schmidt, keeping the Hessenberg matrix
    h(i) = q' * Q(:, i);
    q = q - h(i) * Q(:, i);
  end
  h(k + 1) = norm(q);
  q = q / h(k + 1);
end

%---------------------------------------------------------------------%
%                  Applying Givens Rotation to H col                  %
%---------------------------------------------------------------------%
function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
  % apply for ith column
  for 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;
  end

  % update the next sin cos values for rotation
  [cs_k sn_k] = givens_rotation(h(k), h(k + 1));

  % eliminate H(i + 1, i)
  h(k) = cs_k * h(k) + sn_k * h(k + 1);
  h(k + 1) = 0.0;
end

%%----Calculate the Given rotation matrix----%%
function [cs, sn] = givens_rotation(v1, v2)
%  if (v1 == 0)
%    cs = 0;
%    sn = 1;
%  else
    t = sqrt(v1^2 + v2^2);
%    cs = abs(v1) / t;
%    sn = cs * v2 / v1;
    cs = v1 / t;  % see http://www.netlib.org/eispack/comqr.f
    sn = v2 / t;
%  end
end

Ver también

• Método de gradiente biconjugado

Referencias

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 , 2ª 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 , 3ª 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