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²⁵²⁵endnote: ²⁵Desenvolvido por Yousef Saad e H. Schultz, 1986. Fonte: Wikipedia.) é um método de subespaço de Krylov²⁶²⁶endnote: ²⁶Alexei Nikolajewitsch Krylov, 1863 - 1945, engenheiro e matemático russo. Fonte: Wikipédia. 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.179)

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

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

(1.180)

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,\boldsymbol{r}^{(0)}\right)=\operatorname{span}\left\{% \boldsymbol{r}^{(0)},A\boldsymbol{r}^{(0)},A^{2}\boldsymbol{r}^{(0)},\dotsc,A^% {m-1}\boldsymbol{r}^{(0)}\right\},

(1.181)

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

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

(1.182)

sendo $\boldsymbol{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}\boldsymbol{b}\approx\boldsymbol{x}^{(m)}=\boldsymbol{x}^{(0)}+q_{m-1}(A% )\boldsymbol{r}^{(0)},

(1.183)

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

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

(1.184)

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,\boldsymbol{v}^{(0)})=\operatorname{span}\{% \boldsymbol{v}^{(0)},A\boldsymbol{v}^{(0)},\dotsc,A^{m-1}\boldsymbol{v}^{(0)}\},

(1.185)

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

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

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

(1.186)

onde, $V_{m}=\left[\boldsymbol{v}^{(0)},\dotsc,\boldsymbol{v}^{(m-1)}\right]$ é a matriz $n\times m$ cujas colunas formam uma base ortonormal $\{\boldsymbol{v}^{(0)},\dotsc,\boldsymbol{v}^{(m-1)}\}$ de $\mathcal{K}_{m}$ e $\boldsymbol{y}\in R^{m}$ . Para computar esta base, podemos usar o método de Gram²⁹²⁹endnote: ²⁹Jørgen Pedersen Gram, 1850 - 1916, matemático dinamarquês. Fonte: Wikipédia.-Schmidt³⁰³⁰endnote: ³⁰Erhard Schmidt, 1876 - 1959, matemático alemão. Fonte: Wikipédia. Arnoldi³¹³¹endnote: ³¹Walter Edwin Arnoldi, 1917 - 1995, engenheiro americano estadunidense. Fonte: Wikipédia.-modificado [9, 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³²³²endnote: ³²Karl Adolf Hessenberg, 1904 - 1959, engenheiro e matemático alemão. Fonte: Wikipédia. cujas entradas não nulas são computadas pelo algoritmo acima (Passos 2(a)i-ii). Definimos

	$\displaystyle J(\boldsymbol{y})=\\|\boldsymbol{b}-A\boldsymbol{x}\\|_{2},$		(1.187)
	$\displaystyle\text{}\quad=\left\\|\boldsymbol{b}-A(\boldsymbol{x}^{(0)}+V_{m}% \boldsymbol{y})\right\\|_{2},$		(1.188)
	$\displaystyle\text{}\quad=\left\\|\boldsymbol{r}^{(0)}-AV_{m}\boldsymbol{y}% \right\\|_{2},$		(1.189)
	$\displaystyle\text{}\quad=\left\\|\beta\boldsymbol{v}^{(0)}-V_{m+1}\bar{H}_{m}% \boldsymbol{y}\right\\|_{2},$		(1.190)
	$\displaystyle\text{}\quad=\left\\|V_{m+1}(\beta\boldsymbol{e}^{(0)}-\bar{H}_{m}% \boldsymbol{y})\right\\|,$		(1.191)

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

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

(1.192)

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

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

(1.193)

onde

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

(1.194)

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

⬇

1from numpy.linalg import norm, lstsq

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

4 m = min(m, b.size)

5 n = b.size

6 x = x0.copy()

7 norm_b = norm(b)

9 info = 0

10 # inicializa a base de Arnoldi

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

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

14 r = b - A @ x

15 beta = norm(r)

16 V[:,0] = r / beta

18 for j in range(m):

19 w = A @ V[:,j]

20 for i in range(j+1):

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

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

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

25 if H[j+1,j] != 0:

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

27 else:

28 m = j+1

29 break

31 # resolução do sistema menor

32 e1 = np.zeros(j+2)

33 e1[0] = beta

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

36 # atualização da solução

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

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

40 info = 1

41 break

43 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+\boldsymbol{a}\cdot\nabla u=f,\quad\text{em }% \Omega=(0,1)\times(0,1),$		(1.195)
	$\displaystyle u=0,\quad\text{em }\partial\Omega,$		(1.196)

onde $\epsilon=1$ é o coeficiente de difusão, $\boldsymbol{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ário}.\end{cases}

(1.197)

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_{i})$ , onde $x_{i}=ih$ e $y_{j}=jh$ , para $i,j=1,2,\dotsc,n$ . Aplicando um espuema 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.198)
	$\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.199)

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\boldsymbol{u}=\boldsymbol{f},

(1.200)

onde $A$ é uma matriz esparsa de ordem $N=(n-2)^{2}$ , $\boldsymbol{u}\in\mathbb{R}^{N}$ é o vetor das incógnitas e $\boldsymbol{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	$\boldsymbol{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 ao menos $n$ passos.

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³³³³endnote: ³³Alston Scott Householder, 1904 - 1993, matemático americano estadunidense. Fonte: Wikipédia. é utilizada. Consulte mais em [9, Subsetion 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\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.201)

e $\boldsymbol{b}=(3,4)$ . Usando uma aproximação inicial $\boldsymbol{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 $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{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\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.202)

e $\boldsymbol{b}=(3,3)$ . Usando uma aproximação inicial $\boldsymbol{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 $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{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\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.203)

e $\boldsymbol{b}=(3,4.1)$ . Usando uma aproximação inicial $\boldsymbol{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 $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{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)

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

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

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

Analise a convergência do método para diferentes tamanho 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