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

Em matemática, o método dos resíduos mínimos generalizados (GMRES) é um método iterativo para a solução numérica de um sistema não simétrico indefinido de equações lineares . O método aproxima a solução pelo vetor em um subespaço de Krylov com resíduo mínimo . A iteração de Arnoldi é usada para encontrar esse vetor.

O método GMRES foi desenvolvido por Yousef Saad e Martin H. Schultz em 1986. É uma generalização e melhoria do método MINRES devido a Paige e Saunders em 1975. O método MINRES requer que a matriz seja simétrica, mas tem a vantagem de requer apenas o manuseio de três vetores. GMRES é um caso especial do método DIIS desenvolvido por Peter Pulay em 1980. DIIS é aplicável a sistemas não lineares.

O método

Denote a norma euclidiana de qualquer vetor v por . Denote o sistema (quadrado) de equações lineares a ser resolvido por ${\ displaystyle \ | v \ |}$

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

A matriz A é considerada invertível de tamanho m- by- m . Além disso, assume-se que b é normalizado, ou seja, isso . ${\ displaystyle \ | b \ | = 1}$

O n -ésimo subespaço Krylov para este problema é

${\ displaystyle K_ {n} = K_ {n} (A, r_ {0}) = \ operatorname {span} \, \ {r_ {0}, Ar_ {0}, A ^ {2} r_ {0}, \ ldots, A ^ {n-1} r_ {0} \}. \,}$

onde é o erro inicial dado uma estimativa inicial . Claramente se . ${\ displaystyle r_ {0} = b-Ax_ {0}}$${\ displaystyle x_ {0} \ neq 0}$${\ displaystyle r_ {0} = b}$${\ displaystyle x_ {0} = 0}$

GMRES aproxima a solução exata de pelo vetor que minimiza a norma euclidiana do resíduo . ${\ displaystyle Ax = b}$${\ displaystyle x_ {n} \ in K_ {n}}$ ${\ displaystyle r_ {n} = b-Ax_ {n}}$

Os vetores podem ser quase linearmente dependentes , portanto, em vez dessa base, a iteração de Arnoldi é usada para encontrar vetores ortonormais que formam uma base para . Em 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}}$

Portanto, o vetor pode ser escrito como com , onde é a matriz m -by- 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}}$

O processo de Arnoldi também produz uma matriz ( ) -by- Upper Hessenberg com ${\ displaystyle n + 1}$${\ displaystyle n}$ ${\ displaystyle {\ tilde {H}} _ {n}}$

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

Para matrizes simétricas, uma matriz tri-diagonal simétrica é realmente obtida, resultando no método de minres .

Como as colunas de são ortonormais, temos ${\ 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} {\ til {H}} _ {n} y_ {n} \ | = \ | Q_ {n + 1} (\ beta e_ {1} - {\ tilde {H}} _ {n} y_ {n}) \ | = \ | \ beta e_ {1} - {\ tilde {H}} _ {n} y_ {n} \ |, \,}$

Onde

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

é o primeiro vetor na base padrão de , e ${\ displaystyle \ mathbb {R} ^ {n + 1}}$

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

${\ displaystyle x_ {0}}$sendo o primeiro vetor de teste (geralmente zero). Portanto, pode ser encontrado minimizando a norma euclidiana do resíduo ${\ displaystyle x_ {n}}$

${\ displaystyle r_ {n} = {\ tilde {H}} _ {n} y_ {n} - \ beta e_ {1}.}$

Este é um problema de mínimos quadrados lineares de tamanho n .

Isso produz o método GMRES. Na -ésima iteração: ${\ displaystyle n}$

1. calcular com o método de Arnoldi;${\ displaystyle q_ {n}}$
2. encontre o que minimiza ;${\ displaystyle y_ {n}}$${\ displaystyle \ | r_ {n} \ |}$
3. compute ;${\ displaystyle x_ {n} = x_ {0} + Q_ {n} y_ {n}}$
4. repita se o resíduo ainda não for pequeno o suficiente.

A cada iteração, um produto matriz-vetor deve ser calculado. Isso custa sobre operações de ponto flutuante para matrizes densas gerais de tamanho , mas o custo pode diminuir para matrizes esparsas . Em adição ao produto de matriz-vector, as operações de vírgula flutuante deve ser calculado no N iteração -ésimo. ${\ displaystyle Aq_ {n}}$${\ displaystyle 2m ^ {2}}$ ${\ displaystyle m}$${\ displaystyle O (m)}$${\ displaystyle O (nm)}$

Convergência

O n th iteração minimiza o residual no subespaço Krylov . Como cada subespaço está contido no próximo subespaço, o resíduo não aumenta. Após m iterações, onde m é o tamanho da matriz A , o espaço de Krylov K _m é o todo de R ^m e, portanto, o método GMRES chega à solução exata. No entanto, a ideia é que após um pequeno número de iterações (em relação a m ), o vetor x _n já é uma boa aproximação da solução exata. ${\ displaystyle K_ {n}}$

Isso não acontece em geral. De fato, um teorema de Greenbaum, Pták e Strakoš afirma que para cada sequência não crescente a ₁ ,…, a _{m −1} , a _m = 0, pode-se encontrar uma matriz A tal que || r _n || = a _n para todo n , onde r _n é o resíduo definido acima. Em particular, é possível encontrar uma matriz para a qual o resíduo permanece constante por m  - 1 iterações, e só cai para zero na última iteração.

Na prática, porém, o GMRES geralmente tem um bom desempenho. Isso pode ser comprovado em situações específicas. Se a parte simétrica de A , isto é , é definida positiva , então ${\ 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} \ |,}$

onde e denotam o menor e o maior autovalor da matriz , respectivamente. ${\ displaystyle \ lambda _ {\ mathrm {min}} (M)}$${\ displaystyle \ lambda _ {\ mathrm {max}} (M)}$${\ displaystyle M}$

Se A é simétrico e definido positivo, então temos até

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

onde denota o número de condição de A na norma euclidiana. ${\ displaystyle \ kappa _ {2} (A)}$

No caso geral, onde A não é definido positivo, temos

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

onde P _n indica o conjunto de polinómios de grau no máximo n com p (0) = 1, V é a matriz que aparecem no decomposição espectral de um , e σ ( A ) é o espectro de Uma . Grosso modo, isso diz que a convergência rápida ocorre quando os autovalores de A estão agrupados longe da origem e A não está muito longe da normalidade .

Todas essas desigualdades limitam apenas os resíduos em vez do erro real, ou seja, a distância entre a iteração atual x _n e a solução exata.

Extensões do método

Como outros métodos iterativos, GMRES é geralmente combinado com um método de pré - condicionamento para acelerar a convergência.

O custo das iterações cresce como O ( n ² ), onde n é o número da iteração. Portanto, o método às vezes é reiniciado após um número, digamos k , de iterações, com x _k como estimativa inicial. O método resultante é denominado GMRES ( k ) ou GMRES reiniciado. Para matrizes definidas não positivas, este método pode sofrer estagnação na convergência, pois o subespaço reiniciado está frequentemente próximo ao subespaço anterior.

As deficiências do GMRES e do GMRES reiniciado são resolvidas pela reciclagem do subespaço Krylov nos métodos do tipo GCRO, como GCROT e GCRODR. A reciclagem de subespaços de Krylov no GMRES também pode acelerar a convergência quando sequências de sistemas lineares precisam ser resolvidas.

Comparação com outros solucionadores

A iteração de Arnoldi se reduz à iteração de Lanczos para matrizes simétricas. O método do subespaço de Krylov correspondente é o método residual mínimo (MinRes) de Paige e Saunders. Ao contrário do caso assimétrico, o método MinRes é dado por uma relação de recorrência de três termos . Pode-se mostrar que não existe um método de subespaço de Krylov para matrizes gerais, que é dado por uma relação de recorrência curta e ainda minimiza as normas dos resíduos, como faz o GMRES.

Outra classe de métodos baseia-se na iteração assimétrica de Lanczos , em particular o método BiCG . Estes usam uma relação de recorrência de três termos, mas não atingem o resíduo mínimo e, portanto, o resíduo não diminui monotonicamente para esses métodos. A convergência nem mesmo é garantida.

A terceira classe é formada por métodos como CGS e BiCGSTAB . Eles também funcionam com uma relação de recorrência de três termos (portanto, sem otimização) e podem até terminar prematuramente sem atingir a convergência. A ideia por trás desses métodos é escolher os polinômios geradores da sequência de iteração de maneira adequada.

Nenhuma dessas três classes é a melhor para todas as matrizes; sempre há exemplos em que uma classe supera a outra. Portanto, vários solucionadores são testados na prática para ver qual é o melhor para um determinado problema.

Resolvendo o problema dos mínimos quadrados

Uma parte do método GMRES é encontrar o vetor que minimiza ${\ displaystyle y_ {n}}$

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

Observe que é uma matriz ( n  + 1) -by- n , portanto, ela fornece um sistema linear excessivamente restrito de n +1 equações para n incógnitas. ${\ displaystyle {\ tilde {H}} _ {n}}$

O mínimo pode ser calculado usando uma decomposição QR : encontre uma matriz ortogonal ( n  + 1) -by- ( n  + 1) Ω _n e uma matriz triangular superior ( n  + 1) -by- n de modo que ${\ displaystyle {\ tilde {R}} _ {n}}$

${\ displaystyle \ Omega _ {n} {\ tilde {H}} _ {n} = {\ tilde {R}} _ {n}.}$

A matriz triangular possui uma linha a mais do que colunas, portanto, sua linha inferior consiste em zero. Portanto, pode ser decomposto como

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

onde é uma matriz triangular n -by- n (portanto quadrada). ${\ displaystyle R_ {n}}$

A decomposição QR pode ser atualizada de forma barata de uma iteração para a próxima, porque as matrizes de Hessenberg diferem apenas por uma linha de zeros e uma coluna:

${\ displaystyle {\ tilde {H}} _ {n + 1} = {\ begin {bmatrix} {\ tilde {H}} _ {n} & h_ {n + 1} \\ 0 & h_ {n + 2, n + 1} \ end {bmatrix}},}$

onde h _{n + 1} = ( h _{1, n + 1} , ..., h _{n + 1, n + 1} ) ^T . Isso implica que a pré- multiplicação da matriz de Hessenberg com Ω _n , aumentada com zeros e uma linha com identidade multiplicativa, produz quase uma 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}}}$

Isso seria triangular se σ fosse zero. Para remediar isso, é necessária a rotação Givens

${\ displaystyle G_ {n} = {\ begin {bmatrix} I_ {n} & 0 & 0 \\ 0 & c_ {n} & s_ {n} \\ 0 & -s_ {n} & c_ {n} \ end {bmatrix}}}$

Onde

${\ displaystyle c_ {n} = {\ frac {\ rho} {\ sqrt {\ rho ^ {2} + \ sigma ^ {2}}}} \ quad {\ mbox {and}} \ quad s_ {n} = {\ frac {\ sigma} {\ sqrt {\ rho ^ {2} + \ sigma ^ {2}}}}.}$

Com esta rotação Givens, formamos

${\ displaystyle \ Omega _ {n + 1} = G_ {n} {\ begin {bmatrix} \ Omega _ {n} & 0 \\ 0 & 1 \ end {bmatrix}}.}$

De fato,

${\ 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 {com}} \ quad r_ {n + 1, n + 1} = {\ sqrt {\ rho ^ {2} + \ sigma ^ {2 }}}}$

é uma matriz triangular.

Dada a decomposição QR, o problema de minimização é facilmente resolvido observando 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 o vetor por ${\ displaystyle \ beta \ Omega _ {n} e_ {1}}$

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

com g _n ∈ R ⁿ e γ _n ∈ R , isso é

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

O vetor y que minimiza esta expressão é dado por

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

Novamente, os vetores são fáceis de atualizar. ${\ displaystyle g_ {n}}$

Código de exemplo

GMRES regular (MATLAB / GNU Octave)

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

Veja também

• Método de gradiente biconjugado

Referências

Notas

• A. Meister, Numerik linearer Gleichungssysteme , 2ª edição, Vieweg 2005, ISBN  978-3-528-13135-7 .
• Y. Saad, Iterative Methods for Sparse Linear Systems , 2ª edição, Society for Industrial and Applied Mathematics , 2003. ISBN  978-0-89871-534-7 .
• Y. Saad e MH Schultz, "GMRES: A generalized minimal residual algorithm for resolvendo nonsymmetric linear systems", SIAM J. Sei. Estado. Comput. , 7 : 856–869, 1986. doi : 10.1137 / 0907058 .
• SC Eisenstat, HC Elman e MH Schultz, "métodos iterativos variacionais para sistemas não simétricos de equações lineares", SIAM Journal on Numerical Analysis , 20 (2), 345-357, 1983.
• J. Stoer e R. Bulirsch, Introdução à análise numérica , 3ª edição, Springer, Nova York, 2002. ISBN  978-0-387-95452-3 .
• Lloyd N. Trefethen e David Bau, III, Numerical Linear Algebra , Society for Industrial and Applied Mathematics, 1997. ISBN  978-0-89871-361-9 .
• Dongarra et al. , Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods , 2ª Edição, SIAM, Filadélfia, 1994
• Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). "Reciclagem de subespaços de Krylov para aplicações CFD e um novo solucionador de reciclagem de híbridos". Journal of Computational Physics 303: 222. doi: 10.1016 / j.jcp.2015.09.040