Matemática numérica III

1 Sistemas Lineares 1.5 Métodos de projeção 1.7 Precondicionamento

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

1.6 Método GMRES

O GMRES (do inglês, Generalized Minimal Residual Method²⁵²⁵25Desenvolvido por Yousef Saad e H. Schultz, 1986. Fonte: Wikipedia.) é um método de subespaço de Krylov²⁶²⁶26Alexei Nikolajewitsch Krylov, 1863 - 1945, engenheiro e matemático russo. Fonte: Wikipédia: Alexei Krylov. e é considerado uma das mais eficientes técnicas para a resolução de sistemas lineares gerais e de grande porte (esparsos).

Métodos de subespaço de Krylov

A ideia básica é resolver o sistema linear

Ax=b

(1.190)

por um método de projeção (lembremos da Seção 1.5). Isto é, buscamos uma solução aproximada $\bm{x}^{(m)}\in\mathbb{R}^{n}$ no subespaço afim $\bm{x}^{(0)}+\mathcal{K}_{m}$ de dimensão $m\leq n$ , impondo-se a condição de Petrov²⁷²⁷27Georgi Iwanowitsch Petrov, 1912 - 1987, engenheiro soviético. Fonte: Wikipedia: Georgi Iwanowitsch Petrov.-Galerkin²⁸²⁸28Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia: Boris Galerkin.

b-Ax_{m}\perp\mathcal{L}_{m},

(1.191)

onde $\mathcal{L}_{m}$ também é um subespaço de dimensão $m$ . Quando $\mathcal{K}_{m}$ é um subespaço de Krylov, i.e.

\mathcal{K}_{m}\left(A,\bm{r}^{(0)}\right)=\operatorname{span}\left\{\bm{r}^{(% 0)},A\bm{r}^{(0)},A^{2}\bm{r}^{(0)},\dotsc,A^{m-1}\bm{r}^{(0)}\right\},

(1.192)

temos um método de subespaço de Krylov, onde $\bm{r}^{(0)}$ é o resíduo inicial, i.e.

\bm{r}^{(0)}=b-A\bm{x}^{(0)},

(1.193)

sendo $\bm{x}^{(0)}$ uma aproximação inicial para a solução do sistema. Notamos que com isso, temos que a aproximação calculada é tal que

A^{-1}\bm{b}\approx\bm{x}^{(m)}=\bm{x}^{(0)}+q_{m-1}(A)\bm{r}^{(0)},

(1.194)

onde $q_{m-1}$ é um dado polinômio de grau $m-1$ . No caso particular de $\bm{x}^{(0)}=\bm{0}$ , temos

A^{-1}\bm{b}\approx q_{m-1}(A)b.

(1.195)

Diferentes versões deste método são obtidas pelas escolhas do subespaço $\mathcal{L}_{m}$ e formas de precondicionamento do sistema.

1.6.1 O GMRES básico

O GMRES é o método de subespaço de Krylov em que se assume $\mathcal{L}_{m}=A\mathcal{K}_{m}$ e

\mathcal{K}_{m}=\mathcal{K}_{m}(A,\bm{v}^{(0)})=\operatorname{span}\{\bm{v}^{(% 0)},A\bm{v}^{(0)},\dotsc,A^{m-1}\bm{v}^{(0)}\},

(1.196)

onde $\bm{v}^{(0)}=\bm{r}^{(0)}/\left\|\bm{r}^{(0)}\right\|$ é o vetor normalizado do resíduo inicial $\bm{r}^{(0)}=b-A\bm{x}^{(0)}$ para uma dada aproximação inicial $\bm{x}^{(0)}$ da solução do sistema $Ax=b$ .

Vamos derivar o método observando que qualquer vetor $x$ em $\bm{x}^{(0)}+\mathcal{K}_{m}$ pode ser escrito como segue

\bm{x}=\bm{x}^{(0)}+V_{m}\bm{y}

(1.197)

onde, $V_{m}=\left[\bm{v}^{(0)},\dotsc,\bm{v}^{(m-1)}\right]$ é a matriz $n\times m$ cujas colunas formam uma base ortonormal $\{\bm{v}^{(0)},\dotsc,\bm{v}^{(m-1)}\}$ de $\mathcal{K}_{m}$ e $\bm{y}\in R^{m}$ . Para computar esta base, podemos usar o método de Gram²⁹²⁹29Jørgen Pedersen Gram, 1850 - 1916, matemático dinamarquês. Fonte: Wikipédia: Jørgen Pedersen Gram.-Schmidt³⁰³⁰30Erhard Schmidt, 1876 - 1959, matemático alemão. Fonte: Wikipédia: Erhard Schmidt. Arnoldi³¹³¹31Walter Edwin Arnoldi, 1917 - 1995, engenheiro americano estadunidense. Fonte: Wikipédia: Walter Edwin Arnoldi.-modificado [11, Subseção 6.3]:

1.

Dado $v_{1}$ de norma 1
2.
Para $j=1,\dotsc,m$ :
1. (a)
  
  $w_{j}\leftarrow Av_{j}$
2. (b)
  Para $i=1,\dotsc,j$ :
  1. i.
    
    $h_{i,j}\leftarrow(w_{j},v_{i})$
  2. ii.
    
    $w_{j}\leftarrow w_{j}-h_{i,j}v_{i}$
3. (c)
  
  $h_{j+1,j}\leftarrow\|w_{j}\|_{2}$
4. (d)
  
  Se $h_{j+1,j}=0$ , então pare.
5. (e)
  
  $v_{j+1}\leftarrow w_{j}/h_{j+1,j}$

Seja, então, $\bar{H}_{m}=[h_{i,j}]_{i,j=1}^{m+1,m}$ a matriz de Hessenberg³²³²32Karl Adolf Hessenberg, 1904 - 1959, engenheiro e matemático alemão. Fonte: Wikipédia: Karl Hessenberg. cujas entradas não nulas são computadas pelo algoritmo acima (Passos 2(a)i-ii). Definimos

	$\displaystyle J(\bm{y})=\\|\bm{b}-A\bm{x}\\|_{2},$		(1.198)
	$\displaystyle\text{}\quad=\left\\|\bm{b}-A(\bm{x}^{(0)}+V_{m}\bm{y})\right\\|_{2},$		(1.199)
	$\displaystyle\text{}\quad=\left\\|\bm{r}^{(0)}-AV_{m}\bm{y}\right\\|_{2},$		(1.200)
	$\displaystyle\text{}\quad=\left\\|\beta\bm{v}^{(0)}-V_{m+1}\bar{H}_{m}\bm{y}% \right\\|_{2},$		(1.201)
	$\displaystyle\text{}\quad=\left\\|V_{m+1}(\beta\bm{e}^{(0)}-\bar{H}_{m}\bm{y})% \right\\|,$		(1.202)

onde $\bm{e}^{(0)}$ é o vetor canônico $(1,0,\dotsc,0)^{T}\in\mathbb{R}^{m+1}$ e $\beta=\|\bm{r}^{(0)}\|_{2}$ . Uma vez que $V_{m+1}$ é uma matriz ortonormal, temos

J(\bm{y})=\left\|\beta\bm{e}^{(0)}-\bar{H}_{m}\bm{y}\right\|_{2}.

(1.203)

A aproximação GMRES é então obtida como

\bm{x}^{(m)}=\bm{x}^{(0)}+V_{m}\bm{y}^{(m)},

(1.204)

onde

\bm{y}^{(m)}=\arg\min_{\bm{y}\in\mathbb{R}^{m}}\left\|\beta\bm{e}^{(0)}-\bar{H% }_{m}\bm{y}\right\|_{2}.

(1.205)

Observamos que este último é um pequeno problema de minimização, sendo que requer a solução de um sistema $(m+1)\times m$ de mínimos quadrados, sendo $m$ normalmente pequeno.

Código 10: gmres_basic.py

⬇

1import numpy as np

2from numpy.linalg import norm, lstsq

4def gmres_basic(A, b, x0, m=50, rtol=1e-5, atol=0.0):

5 m = min(m, b.size)

6 n = b.size

7 x0 = x0.copy()

8 norm_b = norm(b)

9 r = b - A @ x0

10 beta = norm(r)

11 if beta <= max(rtol*norm_b, atol):

12 return x0, 1, 0

14 V = np.zeros((n, m+1))

15 H = np.zeros((m+1, m))

16 V[:, 0] = r / beta

17 info = 0

19 for j in range(m):

20 w = A @ V[:, j]

21 for i in range(j+1):

22 H[i, j] = np.dot(w, V[:, i])

23 w = w - H[i, j] * V[:, i]

25 H[j+1, j] = norm(w)

26 if H[j+1, j] != 0 and j + 1 < m:

27 V[:, j+1] = w / H[j+1, j]

29 e1 = np.zeros(j+2)

30 e1[0] = beta

31 y, _, _, _ = lstsq(H[:j+2, :j+1], e1, rcond=None)

32 x = x0 + V[:, :j+1] @ y

34 if norm(b - A @ x) <= max(rtol*norm_b, atol):

35 info = 1

36 return x, info, j+1

37 if H[j+1, j] == 0:

38 break

40 return x, info, j+1

Exemplo 1.6.1.(Problema de difusão-advecção 2D)

Consideremos o seguinte problema de difusão-advecção 2D

	$\displaystyle-\epsilon\Delta u+\bm{a}\cdot\nabla u=f,\quad\text{em }\Omega=(0,% 1)\times(0,1),$		(1.206)
	$\displaystyle u=0,\quad\text{em }\partial\Omega,$		(1.207)

onde $\epsilon=1$ é o coeficiente de difusão, $\bm{a}=(20,1)$ é o campo de advecção e a fonte é dada por

f(x)=\begin{cases}100,&\text{se }0.4\leq x,y\leq 0.6,\\ 0,&\text{caso contr\'{a}rio}.\end{cases}

(1.208)

Consulte o Exemplo 1.5.1 para mais detalhes.

Assumimos uma malha espacial uniforme de $n\times n$ , com tamanho de malha $h=1/(n-1)$ . Denotamos $u_{i,j}\approx u(x_{i},y_{j})$ , onde $x_{i}=ih$ e $y_{j}=jh$ , para $i,j=1,2,\dotsc,n$ . Aplicando um esquema upwind para a advecção e diferenças finitas centrais para a difusão, obtemos o seguinte esquema de diferenças finitas

	$\displaystyle-\epsilon\left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+\frac{u% _{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}\right)$		(1.209)
	$\displaystyle\text{}\quad+a_{1}\frac{u_{i,j}-u_{i-1,j}}{h}+a_{2}\frac{u_{i,j+1% }-u_{i,j}}{h}=f_{i,j},$		(1.210)

para $i,j=2,3,\dotsc,n-1$ , onde $f_{i,j}=f(x_{i},y_{j})$ . As condições de contorno são dadas por $u_{1,j}=u_{n,1}=u_{i,1}=u_{i,n}=0$ , para $i,j=1,2,\dotsc,n-1$ . Por fim, consideramos a enumeração dos nodos da malha $k=i-1+(j-2)(n-2)$ para obtermos o sistema linear

A\bm{u}=\bm{f},

(1.211)

onde $A$ é uma matriz esparsa de ordem $N=(n-2)^{2}$ , $\bm{u}\in\mathbb{R}^{N}$ é o vetor das incógnitas e $\bm{f}\in\mathbb{R}^{N}$ é o vetor da fonte.

Observando que $A$ é apenas positiva definida (i.e. $A+A^{T}$ é simétrica positiva definida), aplicamos a iteração do mínimo resíduo para resolver o sistema linear. A seguinte tabela mostra o número de iterações necessárias para a convergência do método, para diferentes tamanhos de malha. O critério de parada é dado por $\textrm{rtol}=10^{-5}$ e $\textrm{atol}=0$ e número máximo de iterações $m=50$ . Verifique!

$n$	GMRES	$\bm{r}$
11	24	$8.43\times 10^{-6}$
21	46	$9.52\times 10^{-6}$
41	50	$1.41\times 10^{-3}$

Observação 1.6.1.(Convergência)

Pode-se mostrar que o GMRES converge em no máximo $n$ passos, desconsiderando os erros de arredondamento.

Observação 1.6.2.(GMRES com a ortogonalização de Householder)

No algoritmo acima, o método Arnoldi-modificado de Gram-Schmidt é utilizado. Uma versão numericamente mais eficiente é obtida quando a transformação de Householder³³³³33Alston Scott Householder, 1904 - 1993, matemático americano estadunidense. Fonte: Wikipédia: Alston Scott Householder. é utilizada. Consulte mais em [11, Subseção 6.5.2].

Observação 1.6.3.(GMRES com Reinicialização)

O GMRES com reinicialização é uma variação do método para sistemas que requerem uma aproximação GMRES $x_{m}$ com $m$ grande. Nestes casos, o método original pode demandar um custo muito alto de memória computacional. A alternativa consiste em assumir $m$ pequeno e, caso não suficiente, recalcular a aproximação GMRES com $x_{0}=x_{m}$ . Este algoritmo pode ser descrito como segue.

Exercícios

E. 1.6.1.

Aplique o método GMRES para resolver o sistema linear $A\bm{x}=\bm{b}$ , com

A=\begin{bmatrix}2&1\\ 1&3\end{bmatrix}

(1.212)

e $\bm{b}=(3,4)$ . Usando uma aproximação inicial $\bm{x}^{(0)}=(0.5,0)$ , faça uma análise geométrica das iterações.

Dica: faça um gráfico de contorno da norma $\|\bm{r}(\bm{x})\|_{2}=\|\bm{b}-A\bm{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.6.2.

Aplique o método GMRES para resolver o sistema linear $A\bm{x}=\bm{b}$ , com

A=\begin{bmatrix}2&1\\ 0&3\end{bmatrix}

(1.213)

e $\bm{b}=(3,3)$ . Usando uma aproximação inicial $\bm{x}^{(0)}=(0,0.5)$ , faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno da norma $\|\bm{r}(\bm{x})\|_{2}=\|\bm{b}-A\bm{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.6.3.

Aplique o método GMRES para resolver o sistema linear $A\bm{x}=\bm{b}$ , com

A=\begin{bmatrix}2&1\\ 2.1&2\end{bmatrix}

(1.214)

e $\bm{b}=(3,4.1)$ . Usando uma aproximação inicial $\bm{x}^{(0)}=(0.5,0)$ , faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno da norma $\|\bm{r}(\bm{x})\|_{2}=\|\bm{b}-A\bm{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.6.4.

Seguindo o Código 10, implemente o método GMRES com reinicialização. Teste o método para o sistema linear do Exemplo 1.6.1, $\textrm{rtol}=10^{-5}$ , $\textrm{atol}=0$ e número máximo de iterações $m=20$ . Compare os resultados com o GMRES básico para diferentes tamanhos de malha.

E. 1.6.5.

Considere o problema de difusão-advecção 2D do Exemplo 1.6.1, discretizado com o esquema upwind. Aplique a iteração da descida mais íngrime do resíduo para resolver o sistema linear resultante para os seguintes coeficientes de advecção:

a)

$\bm{a}=(-1,1)$ ;
b)

$\bm{a}=(1,-1)$ ;
c)

$\bm{a}=(-1,-1)$ .

Analise a convergência do método para diferentes tamanhos de malha.

Dica: lembre-se que o esquema upwind deve ser adaptado para cada caso.

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Política de uso de dados

Política de uso de dados

Matemática numérica III

1.6 Método GMRES

Métodos de subespaço de Krylov

1.6.1 O GMRES básico

Exemplo 1.6.1.(Problema de difusão-advecção 2D)

Observação 1.6.1.(Convergência)

Observação 1.6.2.(GMRES com a ortogonalização de Householder)

Observação 1.6.3.(GMRES com Reinicialização)

Exercícios

E. 1.6.1.

E. 1.6.2.

E. 1.6.3.

E. 1.6.4.

E. 1.6.5.

Envie seu comentário