| | | |

Matemática Numérica III

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

1.6 Método GMRES

O GMRES (do inglês, Generalized Minimal Residual Method2525endnote: 25Desenvolvido por Yousef Saad e H. Schultz, 1986. Fonte: Wikipedia.) é um método de subespaço de Krylov2626endnote: 26Alexei 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 𝒙(m)n no subespaço afim 𝒙(0)+𝒦m de dimensão mn, impondo-se a condição de Petrov2727endnote: 27Georgi Iwanowitsch Petrov, 1912 - 1987, engenheiro soviético. Fonte: Wikipedia.-Galerkin2828endnote: 28Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia.

bAxmm, (1.180)

onde m também é um subespaço de dimensão m. Quando 𝒦m é um subespaço de Krylov, i.e.

𝒦m(A,𝒓(0))=span{𝒓(0),A𝒓(0),A2𝒓(0),,Am1𝒓(0)}, (1.181)

temos um método de subespaço de Krylov, onde 𝒓(0) é o resíduo inicial, i.e.

𝒓(0)=bA𝒙(0), (1.182)

sendo 𝒙(0) uma aproximação inicial para a solução do sistema. Notamos que com isso, temos que a aproximação calculada é tal que

A1𝒃𝒙(m)=𝒙(0)+qm1(A)𝒓(0), (1.183)

onde qm1 é um dado polinômio de grau m1. No caso particular de 𝒙(0)=𝟎, temos

A1𝒃qm1(A)b. (1.184)

Diferentes versões deste método são obtidas pelas escolhas do subespaço 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 m=A𝒦m e

𝒦m=𝒦m(A,𝒗(0))=span{𝒗(0),A𝒗(0),,Am1𝒗(0)}, (1.185)

onde 𝒗(0)=𝒓(0)/𝒓(0) é o vetor normalizado do resíduo inicial 𝒓(0)=bA𝒙(0) para uma dada aproximação inicial 𝒙(0) da solução do sistema Ax=b.

Vamos derivar o método observando que qualquer vetor x em 𝒙(0)+𝒦m pode ser escrito como segue

𝒙=𝒙(0)+Vm𝒚 (1.186)

onde, Vm=[𝒗(0),,𝒗(m1)] é a matriz n×m cujas colunas formam uma base ortonormal {𝒗(0),,𝒗(m1)} de 𝒦m e 𝒚Rm. Para computar esta base, podemos usar o método de Gram2929endnote: 29Jørgen Pedersen Gram, 1850 - 1916, matemático dinamarquês. Fonte: Wikipédia.-Schmidt3030endnote: 30Erhard Schmidt, 1876 - 1959, matemático alemão. Fonte: Wikipédia. Arnoldi3131endnote: 31Walter Edwin Arnoldi, 1917 - 1995, engenheiro americano estadunidense. Fonte: Wikipédia.-modificado [6, Subseção 6.3]:

  1. 1.

    Dado v1 de norma 1

  2. 2.

    Para j=1,,m:

    1. (a)

      wjAvj

    2. (b)

      Para i=1,,j:

      1. i.

        hi,j(wj,vi)

      2. ii.

        wjwjhi,jvi

    3. (c)

      hj+1,jwj2

    4. (d)

      Se hj+1,j=0, então pare.

    5. (e)

      vj+1wj/hj+1,j

Seja, então, H¯m=[hi,j]i,j=1m+1,m a matriz de Hessenberg3232endnote: 32Karl 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

J(𝒚)=𝒃A𝒙2, (1.187)
=𝒃A(𝒙(0)+Vm𝒚)2, (1.188)
=𝒓(0)AVm𝒚2, (1.189)
=β𝒗(0)Vm+1H¯m𝒚2, (1.190)
=Vm+1(β𝒆(0)H¯m𝒚), (1.191)

onde 𝒆(0) é o vetor canônico (1,0,,0)Tm+1 e β=𝒓(0)2. Uma vez que Vm+1 é uma matriz ortonormal, temos

J(𝒚)=β𝒆(0)H¯m𝒚2. (1.192)

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

𝒙(m)=𝒙(0)+Vm𝒚(m), (1.193)

onde

𝒚(m)=argmin𝒚mβ𝒆(0)H¯m𝒚2. (1.194)

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

Código 10: gmres_basic.py
1from numpy.linalg import norm, lstsq
2
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)
8
9 info = 0
10 # inicializa a base de Arnoldi
11 V = np.zeros((n, m+1))
12 H = np.zeros((m+1, m))
13
14 r = b - A @ x
15 beta = norm(r)
16 V[:,0] = r / beta
17
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]
23
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
30
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)
35
36 # atualização da solução
37 x = x0 + V[:,:j+1] @ y
38
39 if norm(b - A@x) <= max(rtol*norm_b, atol):
40 info = 1
41 break
42
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

ϵΔu+𝒂u=f,em Ω=(0,1)×(0,1), (1.195)
u=0,em Ω, (1.196)

onde ϵ=1 é o coeficiente de difusão, 𝒂=(20,1) é o campo de advecção e a fonte é dada por

f(x)={100,se 0.4x,y0.6,0,caso contrário. (1.197)

Consulte o Exemplo 1.5.1 para mais detalhes.

Assumimos uma malha espacial uniforme de n×n , com tamanho de malha h=1/(n1). Denotamos ui,ju(xi,yi), onde xi=ih e yj=jh, para i,j=1,2,,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

ϵ(ui+1,j2ui,j+ui1,jh2+ui,j+12ui,j+ui,j1h2) (1.198)
+a1ui,jui1,jh+a2ui,j+1ui,jh=fi,j, (1.199)

para i,j=2,3,,n1, onde fi,j=f(xi,yj). As condições de contorno são dadas por u1,j=un,1=ui,1=ui,n=0, para i,j=1,2,,n1. Por fim, consideramos a enumeração dos nodos da malha k=i1+(j2)(n2) para obtermos o sistema linear

A𝒖=𝒇, (1.200)

onde A é uma matriz esparsa de ordem N=(n2)2, 𝒖N é o vetor das incógnitas e 𝒇N é o vetor da fonte.

Observando que A é apenas positiva definida (i.e. A+AT é 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 rtol=105 e atol=0 e número máximo de iterações m=50. Verifique!

n GMRES 𝒓
11 24 8.43×106
21 46 9.52×106
41 50 1.41×103
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 Householder3333endnote: 33Alston Scott Householder, 1904 - 1993, matemático americano estadunidense. Fonte: Wikipédia. é utilizada. Consulte mais em [6, 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 xm 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 x0=xm. Este algoritmo pode ser descrito como segue.

Exercícios

E. 1.6.1.

Aplique o método GMRES para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.201)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0.5,0), faça uma análise geométrica das iterações.

Dica: faça um gráfico de contorno da norma 𝒓(𝒙)2=𝒃A𝒙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𝒙=𝒃, com

A=[2103] (1.202)

e 𝒃=(3,3). Usando uma aproximação inicial 𝒙(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 𝒓(𝒙)2=𝒃A𝒙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𝒙=𝒃, com

A=[212.12] (1.203)

e 𝒃=(3,4.1). Usando uma aproximação inicial 𝒙(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 𝒓(𝒙)2=𝒃A𝒙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, rtol=105, 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:

  1. a)

    𝒂=(1,1);

  2. b)

    𝒂=(1,1);

  3. c)

    𝒂=(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!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
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.

Matemática Numérica III

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

1.6 Método GMRES

O GMRES (do inglês, Generalized Minimal Residual Method2525endnote: 25Desenvolvido por Yousef Saad e H. Schultz, 1986. Fonte: Wikipedia.) é um método de subespaço de Krylov2626endnote: 26Alexei 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 𝒙(m)n no subespaço afim 𝒙(0)+𝒦m de dimensão mn, impondo-se a condição de Petrov2727endnote: 27Georgi Iwanowitsch Petrov, 1912 - 1987, engenheiro soviético. Fonte: Wikipedia.-Galerkin2828endnote: 28Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia.

bAxmm, (1.180)

onde m também é um subespaço de dimensão m. Quando 𝒦m é um subespaço de Krylov, i.e.

𝒦m(A,𝒓(0))=span{𝒓(0),A𝒓(0),A2𝒓(0),,Am1𝒓(0)}, (1.181)

temos um método de subespaço de Krylov, onde 𝒓(0) é o resíduo inicial, i.e.

𝒓(0)=bA𝒙(0), (1.182)

sendo 𝒙(0) uma aproximação inicial para a solução do sistema. Notamos que com isso, temos que a aproximação calculada é tal que

A1𝒃𝒙(m)=𝒙(0)+qm1(A)𝒓(0), (1.183)

onde qm1 é um dado polinômio de grau m1. No caso particular de 𝒙(0)=𝟎, temos

A1𝒃qm1(A)b. (1.184)

Diferentes versões deste método são obtidas pelas escolhas do subespaço 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 m=A𝒦m e

𝒦m=𝒦m(A,𝒗(0))=span{𝒗(0),A𝒗(0),,Am1𝒗(0)}, (1.185)

onde 𝒗(0)=𝒓(0)/𝒓(0) é o vetor normalizado do resíduo inicial 𝒓(0)=bA𝒙(0) para uma dada aproximação inicial 𝒙(0) da solução do sistema Ax=b.

Vamos derivar o método observando que qualquer vetor x em 𝒙(0)+𝒦m pode ser escrito como segue

𝒙=𝒙(0)+Vm𝒚 (1.186)

onde, Vm=[𝒗(0),,𝒗(m1)] é a matriz n×m cujas colunas formam uma base ortonormal {𝒗(0),,𝒗(m1)} de 𝒦m e 𝒚Rm. Para computar esta base, podemos usar o método de Gram2929endnote: 29Jørgen Pedersen Gram, 1850 - 1916, matemático dinamarquês. Fonte: Wikipédia.-Schmidt3030endnote: 30Erhard Schmidt, 1876 - 1959, matemático alemão. Fonte: Wikipédia. Arnoldi3131endnote: 31Walter Edwin Arnoldi, 1917 - 1995, engenheiro americano estadunidense. Fonte: Wikipédia.-modificado [6, Subseção 6.3]:

  1. 1.

    Dado v1 de norma 1

  2. 2.

    Para j=1,,m:

    1. (a)

      wjAvj

    2. (b)

      Para i=1,,j:

      1. i.

        hi,j(wj,vi)

      2. ii.

        wjwjhi,jvi

    3. (c)

      hj+1,jwj2

    4. (d)

      Se hj+1,j=0, então pare.

    5. (e)

      vj+1wj/hj+1,j

Seja, então, H¯m=[hi,j]i,j=1m+1,m a matriz de Hessenberg3232endnote: 32Karl 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

J(𝒚)=𝒃A𝒙2, (1.187)
=𝒃A(𝒙(0)+Vm𝒚)2, (1.188)
=𝒓(0)AVm𝒚2, (1.189)
=β𝒗(0)Vm+1H¯m𝒚2, (1.190)
=Vm+1(β𝒆(0)H¯m𝒚), (1.191)

onde 𝒆(0) é o vetor canônico (1,0,,0)Tm+1 e β=𝒓(0)2. Uma vez que Vm+1 é uma matriz ortonormal, temos

J(𝒚)=β𝒆(0)H¯m𝒚2. (1.192)

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

𝒙(m)=𝒙(0)+Vm𝒚(m), (1.193)

onde

𝒚(m)=argmin𝒚mβ𝒆(0)H¯m𝒚2. (1.194)

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

Código 10: gmres_basic.py
1from numpy.linalg import norm, lstsq
2
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)
8
9 info = 0
10 # inicializa a base de Arnoldi
11 V = np.zeros((n, m+1))
12 H = np.zeros((m+1, m))
13
14 r = b - A @ x
15 beta = norm(r)
16 V[:,0] = r / beta
17
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]
23
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
30
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)
35
36 # atualização da solução
37 x = x0 + V[:,:j+1] @ y
38
39 if norm(b - A@x) <= max(rtol*norm_b, atol):
40 info = 1
41 break
42
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

ϵΔu+𝒂u=f,em Ω=(0,1)×(0,1), (1.195)
u=0,em Ω, (1.196)

onde ϵ=1 é o coeficiente de difusão, 𝒂=(20,1) é o campo de advecção e a fonte é dada por

f(x)={100,se 0.4x,y0.6,0,caso contrário. (1.197)

Consulte o Exemplo 1.5.1 para mais detalhes.

Assumimos uma malha espacial uniforme de n×n , com tamanho de malha h=1/(n1). Denotamos ui,ju(xi,yi), onde xi=ih e yj=jh, para i,j=1,2,,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

ϵ(ui+1,j2ui,j+ui1,jh2+ui,j+12ui,j+ui,j1h2) (1.198)
+a1ui,jui1,jh+a2ui,j+1ui,jh=fi,j, (1.199)

para i,j=2,3,,n1, onde fi,j=f(xi,yj). As condições de contorno são dadas por u1,j=un,1=ui,1=ui,n=0, para i,j=1,2,,n1. Por fim, consideramos a enumeração dos nodos da malha k=i1+(j2)(n2) para obtermos o sistema linear

A𝒖=𝒇, (1.200)

onde A é uma matriz esparsa de ordem N=(n2)2, 𝒖N é o vetor das incógnitas e 𝒇N é o vetor da fonte.

Observando que A é apenas positiva definida (i.e. A+AT é 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 rtol=105 e atol=0 e número máximo de iterações m=50. Verifique!

n GMRES 𝒓
11 24 8.43×106
21 46 9.52×106
41 50 1.41×103
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 Householder3333endnote: 33Alston Scott Householder, 1904 - 1993, matemático americano estadunidense. Fonte: Wikipédia. é utilizada. Consulte mais em [6, 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 xm 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 x0=xm. Este algoritmo pode ser descrito como segue.

Exercícios

E. 1.6.1.

Aplique o método GMRES para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.201)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0.5,0), faça uma análise geométrica das iterações.

Dica: faça um gráfico de contorno da norma 𝒓(𝒙)2=𝒃A𝒙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𝒙=𝒃, com

A=[2103] (1.202)

e 𝒃=(3,3). Usando uma aproximação inicial 𝒙(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 𝒓(𝒙)2=𝒃A𝒙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𝒙=𝒃, com

A=[212.12] (1.203)

e 𝒃=(3,4.1). Usando uma aproximação inicial 𝒙(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 𝒓(𝒙)2=𝒃A𝒙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, rtol=105, 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:

  1. a)

    𝒂=(1,1);

  2. b)

    𝒂=(1,1);

  3. c)

    𝒂=(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!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
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.

Pedro H A Konzen
| | | |