Metodos de gradiente.
Metodos de Krylov
Damian Ginestar Peir
o
Departamento de Matem
atica Aplicada
Universidad Polit
ecnica de Valencia
Curso 2012-2013
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 1 / 48
Indice
1 Metodo de descenso rapido
2 Metodo del gradiente conjugado
3 Metodos de Krylov
Metodo del residuo conjugado
metodo GMRES
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 2 / 48
Metodo de descenso rapido
Resolver Ax = b, con A simetrica y definida positiva (SPD).
Definimos la funcion cuadratica Rn R
1 1
(y ) = (y x)T A(y x) = e T Ae .
2 2
Se tiene (y ) 0 y 6= 0 ( definici
on de matriz SPD).
Error e = y x.
Teorema
La solucion del sistema Ax = b es el mnimo de la funcion (y ).
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 3 / 48
Metodo de descenso rapido
(y ) = 12 (y x)T A(y x) = 12 e T Ae
(yk ) = constant representa un hiperelipsoide en un espacio de
dimension n.
El centro geometrico es la soluci
on x del sistema lineal (mnimo).
Construir una sucesion {yk } tal que limk yk = x.
yk+1 = yk + k pk
Hace falta determinar la direcci
on pk y .
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 4 / 48
Metodo de descenso rapido
Este metodo construye una sucesi on que va hacia el centro del
hiperelipsoide en la direcci
on del gradiente.
El gradiente de en el punto yk es
1 1 T 1
(yk ) = ekT Aek = y Ayk ykT b + x T x = Ayk b = rk
2 2 k 2
Como la direccion del vector gradiente es hacia fuera, la direccion
buscada coincide con el residuo rk en la aproximacion actual.
En consecuencia la nueva aproximaci on es
yk+1 = yk + k rk
donde k es una constante a determinar. C
omo? Minimizando (y )
en la direccion buscada rk .
Es decir, nuestro k es la soluci
on
k = argmin (yk + rk )
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 5 / 48
Metodo de descenso rapido
Desarrollando la funcion (yk + rk ) se tiene un polinomio de segundo
grado en la variable .
(yk + rk ) = (yk + rk x)T A(yk + rk x)
= (yk + rk x)T (Ayk + Ark b)
= (yk + rk x)T (Ark rk )
= (rk ek )T (Ark rk )
= 2 rkT Ark rkT rk + ekT Ark + x T rk
= 2 rkT Ark 2rkT rk + const.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 6 / 48
Metodo de descenso rapido
Como rkT Ark > 0 el mnimo de se alcanza cuando
rkT rk
k =
rkT Ark
Otra forma:
Resolver = 0.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 7 / 48
Metodo de descenso rapido
La k + 1 iteracion se puede representar como
rk = b Axk
r T rk
k = Tk
rk Ark
yk+1 = yk + k rk
Notar que el coste computacional es principalmente dos productos
matriz-vector.
De yk+1 = yk + k rk se sigue que
rk+1 = b Axk+1 = b Axk Ak rk = rk k Ark ,
Los residuos consecutivos rk+1 , rk son ortogonales (demostracion:
Ejercicio).
El error ek+1 es A ortogonal a la direcci
on rk . (demostracion:
Ejercicio).
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 8 / 48
Metodo de descenso rapido
Algoritmo: Descenso rapido
Input: y0 , A, b, kmax , tol
r0 = b Ay0 , k = 0
while krk k > tol kbk and k < kmax do
1 z = Ark
r T rk
2 k = kT
z rk
3 yk+1 = yk + k rk
4 rk+1 = rk k z
5 k =k +1
end while
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 9 / 48
Metodo de descenso rapido
Lema
Sea A simetrica definida positiva y sean 0 < n 2 1 sus
valores propios. Si P(t) es un polinomio real, entonces
||P(A)x||A max |P(j )| ||x||A , x Rn
1jn
donde ||x||A = x T Ax.
Teorema
Sean las mismas condiciones que en el lema anterior. La sucesion {yk } del
metodo de descenso rapido satisface
k
1 n
||yk x||A ||y0 x||A
1 + n
donde x es la solucion exacta del sistema.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 10 / 48
Metodo de descenso rapido
Teorema
q q
(A) 1
(yk ) = ekT Aek = kek kA,2 k ke0 kA,2 , donde =
(A) + 1
Cuando los sistemas vienen de discretizar ecuaciones EDPs, (A)
puede ser muy grande.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 11 / 48
Metodo de descenso rapido
Se estima el n
umero de iteraciones para ganar p digitos en la
aproximacion de la soluci
on:
k
kek kA (A) 1
10p resolviendo 10p
ke0 kA (A) + 1
Tomando logaritmos y usando la aproximaci on de primer orden de
(A) 1 2
Taylor log , se obtiene
(A) + 1 (A) + 1
log 10
k p ((A) + 1)
2
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 12 / 48
Metodo del gradiente conjugado
Es una mejora del Descenso rapido. La sucesi
on de recurrencia es
similar
yk+1 = yk + k pk
Las direcciones se construyen como
p0 = r0
pk = rk + k pk1 , k>0
Se exige que las direcciones sean A conjugadas
T
pk1 Apk = 0 ,
es decir, pk y pk1 son A-ortogonales.
Por tanto, se debe cumplir
rkT Apk1
k = T Ap
pk1 k1
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 13 / 48
Metodo del gradiente conjugado
Como en el metodo de descenso mas rapido, la eleccion de k se
obtiene minimizando (yk+1 ) = (yk + k pk ) dando la expresion
rkT pk
k =
pkT Apk
Residuos consecutivos como en el metodo de descenso mas rapido
satisfacen la relacion de recurrencia
rk+1 = rk +k Apk
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 14 / 48
Metodo del gradiente conjugado
Teorema
Las sucesiones de vectores {ri } y {pi } satisfacen las siguientes relaciones
(i) piT rj = 0, 0 0 i <j k,
(ii) riT rj = 0, i 6= j, 0 i, j k,
(iii) piT Apj = 0, i 6= j, 0 i, j k,
(iv) env{r0 , r1 , . . . , rk } = env{p0 , p1 , . . . , pk } = K(A, r0 , k + 1),
donde K(A, r0 , k + 1) = env{r0 , Ar0 , . . . , Ak r0 }.
Corolario
El metodo del gradiente conjugado obtiene la soluci
on del sistema de n
ecuaciones en como maximo n iteraciones del GC.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 15 / 48
Metodo del gradiente conjugado
Otras relaciones u
tiles
1 pkT rk = rkT rk . Ya que de ekT Apj = 0 se sigue rkT pj = 0 y, por tanto,
pkT rk = (rk + k1 pk1 )T rk = rkT rk
2 rkT Apk = pkT Apk .
3 Combinando 1 y 2, se obtiene una definici
on alternativa de k :
rkT pk r T rk
k = T
= Tk
pk Apk rk Apk
1 1 T
4 on alternativa de k . Como pkT Apk = pkT
Formulaci (rk rk+1 ) = r rk
k k k
T T 1 1 T
rk+1 Apk = rk+1 (rk rk+1 ) = rk+1 rk+1
k k
Por tanto
T T
rk+1 pk rk+1 rk+1
k = =
pkT Apk rkT rk
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 16 / 48
Metodo del gradiente conjugado
Algoritmo: Gradiente conjugado
Input: y0 , A, b, kmax , tol
r0 = p0 = b Ax0 , k = 0
while krk k > tol kbk and k < kmax do
1 z = Apk
p T rk
2 k = Tk
z pk
3 yk+1 = yk + k pk
4 rk+1 = rk k z
r T rk+1
5 k = k+1T
rk rk
6 pk+1 = rk+1 + k pk
7 k =k +1
end while
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 17 / 48
Metodo del gradiente conjugado
Ejercicio
Aplicar el algoritmo del gradiente conjugado para el problema
! ! !
2 1 x1 1
=
1 2 x2 0
on inicial x0 = (0, 0)T .
usando como aproximaci
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 18 / 48
Metodo del gradiente conjugado
Soluci
on:
x0 = (0, 0)T .
p0 = r0 = b = (1, 0)T .
1
r0T r0 0 1
0 = p0t Ap0 = 12 , x1 = x0 + 0 p0 = + 1
2 = 2
0 0 0
1 1 2 0
r1 = r0 0 Ap0 = = 1 , r1T r0 = 0
0 2 1 2
1
r1T r1 1 0 1 1 4
0 = r0T r0
= 4, p1 = r1 + 0 p0 = 1 + 4 = 1
2
0 2
r1T r1 2
1 = p1T AP1
= 3
1 1 2
2 2 4 3
x2 = x1 + 1 p1 = + 3 1 = 1
0 2 3
r2 = 0 soluci
on exacta
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 19 / 48
Metodo del gradiente conjugado
Teorema
Sea A Rnn simetrica y definida positiva. Sea x la solucion exacta del
sistema Ax = b. Entonces la sucesion de vectores del Gradiente Conjugado
{yk } cumple
p !k
2 (A) 1
||x yk ||A 2 p ||x y0 ||A
2 (A) + 1
donde 2 (A) = ||A||2 ||A1 ||2 .
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 20 / 48
Metodos de Krylov
Vimos que el medodo de Richardson construye la sucesion
k
X
xk+1 = x0 + (I A)i r0 .
i=0
as xk+1 span{x0 , r0 , Ar0 , A2 r0 , , Ak r0 }.
Para los residuos se tiene
k
X k
X
rk+1 = b Ax0 A (I A)i r0 = r0 A (I A)i r0 .
i=0 i=0
y as rk+1 span{r0 , Ar0 , A2 r0 , , Ak+1 r0 }.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 21 / 48
Metodos de Krylov
El espacio span{r0 , Ar0 , A2 r0 , , Ak1 r0 } se llama subespacio de Krylov
de dimension k, correspondiente a la matriz A y residuo inicial r0
K k (A; r0 ) = span{r0 , Ar0 , A2 r0 , , Ak1 r0 }
Se puede usar la informacion contenida en K k (A; r0 ) de forma mas
eficiente que con el metodo de Richardson?
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 22 / 48
Metodos de Krylov
Por ejemplo, como ya hemos visto, el metodo de descenso mas rapido se
basa en iteraciones de la forma
xk+1 = xk + k rk
rk+1 = b Axk+1 = rk k Ark .
Para este metodo xk+1 x0 K k+1 (A; r0 )
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 23 / 48
Metodos de Krylov
Los metodos de proyeccion de un paso se basan en una combinacion
optima de los dos u
ltimos vectores base del subespacio de Krylov. Es
posible contruir una combinacion lineal
optima de todos los vectores base
del subespacio de Krylov?
La respuesta en dos pasos:
Primero vemos como construir una base para K k (A; r0 );
Despues veremos como construir una aproximaci on optima como una
combinacion lineal de los vectores base. Primero para matrices
simetricas.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 24 / 48
Metodos de Krylov
La base mas simple para el subespacio de Krylov K k (A; r0 ) es la base:
r0 , Ar0 , A2 r0 , Ak1 r0 .
Pero esta base esta mal condicionada ya que Ak1 r0 apuntara cada vez
mas en la direccion del autovector dominante de A.
Una base estable y ortogonal se puede construir usando el metodo de
Arnoldi.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 25 / 48
Metodos de Krylov
Se elige un vector inicial q1 con kq1 k2 = 1.
FOR k = 1, DO (iteraci
on)
v = Aqk
FOR i = 1, ..., k (ortogonalizaci
on)
T
hi,k = v qi
v = v hi,k qi
END FOR
hk+1,k = kv k2
IF hk+1,k = 0 STOP (subespacio invariante)
qk+1 = v /hk+1,k (nuevo vector)
END FOR
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 26 / 48
Metodos de Krylov
El metodo de Arnoldi se puede resumir en:
h1,1 . . . ... h1,k
. ..
h2,1 . .
.
Hk =
.. .. ..
. . .
O hk,k1 hk,k
y Qk = [q1 q2 qk ] entonces
AQk = Qk Hk + hk+1,k qk+1 ekT
onica de Rk .
donde ek es el k esimo vector de la base can
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 27 / 48
Metodos de Krylov
Si A es simetrica de acuerdo a la relaci
on de Arnoldi
QkT AQk = Hk .
Como A es simetrica
HkT = QkT AT Qk = QkT AQk = Hk .
As Hk es simetrica y Hessenberg superior
luego Hk es tridiagonal.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 28 / 48
Metodos de Krylov
As
h1,1 h2,1 O
.. ..
h2,1 . .
Hk = .
.. ..
. . hk,k1
O hk,k1 hk,k
Con k = hk,k y k = hk1,k el metodo de Arnoldi se simplifica en el
metodo de Lanczos. Con el metodo de Lanczos es posible calcular un
nuevo vector base ortogonal utilizando s
olo los dos vectores base previos.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 29 / 48
Metodos de Krylov
El metodo de Arnoldi y el metodo de Lanczos se propusieron originalmente
como metodos iterativos para calcular autovalores de la matriz A:
QkT AQk = Hk
es casi una transformacion de similaridad. Los autovalores de Hk se
llaman los valores de Ritz de A.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 30 / 48
Metodos de Krylov
El metodo de Lanczos nos proporciona un metodo economico para calcular
vectores base ortogonales para el subespacio de Krylov K k (A; r0 ).
Nuestra aproximacion se escribe
xk = x0 + Qk yk
donde yk se determina de forma que o bien se minimiza
f (xk ) = kxk xk2A = (xk x)T A(xk x)
respecto de la norma inducida por A (simaetica y definida positiva) o bien
g (xk ) = kA(xk x)k22 = rkT rk ,
se minimiza la norma del residuo.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 31 / 48
Metodos de Krylov
Se considera la minimizaci
on del error en la A-norma:
xk = x0 + Qk yk f (xk ) = (x0 + Qk yk x)T A(x0 + Qk yk x) .
f (xk )
Se calcula la derivada respecto de yk y yk = 0. Con lo que
QkT AQk yk = QkT r0
y con QkT AQk = Tk , r0 = kr0 kq1 se tiene
Tk yk = kr0 ke1
con e1 el primer vector de la base can
onica.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 32 / 48
Metodos de Krylov
Es facil ver que los residuos son orgonales a los vectores base
rk = r0 AQk yk QkT rk = QkT r0 QkT AQk yk = 0
Esta condicion es equivalente a minimizar f (xk ) cuando A es SPD. En este
caso se obtiene el metodo del gradiente conjugado.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 33 / 48
Metodo del residuo conjugado
El gradiente conjugado minimiza la A-norma del error. Otra forma de
construir una aproximacion
optima xk es minimizar el residuo
g (xk ) = kA(xk x)k22 = rkT rk
sobre todos los xk {x0 K k (A; b)}.
Definiendo
1 2 0
..
2 2 .
.. .. ..
Tk =
. . .
. ..
0 .. . k
k k
k + 1
El metodo de Lanczos se escribe
AQk = Qk+1 T k
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 34 / 48
Metodo del residuo conjugado
El problema es encontrar xk = x0 + Qk yk tal que krk k es mnimo.
rk = b Axk = r0 AQk yk = kr0 kq1 AQk yk
as hay que minimizar
krk k = kkr0 kq1 AQk yk k
= kkr0 kQk+1 e1 Qk+1 T k yk k
= kkr0 ke1 T k yk k
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 35 / 48
Metodo del residuo conjugado
Resolviendo el sistema sobredeterminado T k yk = kr0 ke1 se tienen las
iteraciones
xk = x0 + Qk yk
que minimizan el residuo. El algoritmo resultatnte se llama MINRES (o
residuo conjugado)
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 36 / 48
Metodo del residuo conjugado
r0 = b Ax0 ; p0 = r0 initializacion
FOR k = 0, 1, , DO
r T Ar
k
k = (Apk )T Ap
k k
xk+1 = xk + k pk
rk+1 = rk k Apk actualizacion residuo
T
rk+1 Ark+1
k = rkT Ark
pk+1 = rk+1 + k pk actualizacion direccion
Apk+1 = Ark+1 + k Apk
END FOR
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 37 / 48
Metodo GMRES
En el caso de A simetrica el metodo de Arnoldi es muy eficiente, los
vectores base nuevos se pueden calcular con una recurrencia de tres
terminos.
Esto permite construir metodos muy eficientes que combinan
recurrencias cortas y una condici
on de optimalidad para el error.
Veremos ahora como se pueden construir metodos para el caso no
simetrico. Estos metodos usan el metodo de Arnoldi.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 38 / 48
Metodo GMRES
Metodo de Arnoldi.
Se elige un vector inicial q1 con kq1 k2 = 1.
FOR k = 1, DO (iteracion)
v = Aqk
FOR i = 1, ..., k ortogonalizacion
T
hi,k = v qi
v = v hi,k qi
END FOR
hk+1,k = kv k2
IF hk+1,k = 0 STOP subespacio invariante
qk+1 = v /hk+1,k nuevo vector base
END FOR
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 39 / 48
Metodo GMRES
En forma compacta
h1,1 . . . ... h1,k
. ..
h2,1 . .
.
Hk =
.. .. ..
. . .
O hk,k1 hk,k
y Qk = [q1 q2 qk ] entonces
AQk = Qk Hk + hk+1,k qk+1 ekT
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 40 / 48
Metodo GMRES
Definiendo
h1,1 . . . ... h1,k
. ..
h2,1 . .
.
Hk =
.. .. ..
. . .
hk,k1 hk,k
O hk+1,k
el metodo de Arnoldi se escribe
AQk = Qk+1 H k
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 41 / 48
Metodo GMRES
El metodo de Anoldi nos da una base ortogonal para el subespacio de
Krylov K k (A; r0 ).
Se tiene la aproximacion
xk = x0 + Qk yk
yk es tal que minimiza el error
f (xk ) = kxk xk2A = (xk x)T A(xk x)
respecto de la A-norma, o bien minimiza el
g (xk ) = kA(xk x)k22 = rkT rk ,
residuo.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 42 / 48
Metodo GMRES
Si A no es SPD, la A-norma no esta bien definida.
Imponiendo ahora que los residuos sean ortogonales a Qk se tiene
QkT (r0 AQk yk ) = 0 kr0 ke1 Hk yk = 0
Resolviendo el sistema peque
no
yk = kr0 kHk1 e1 xk = x0 + Qk yk
Este metodo se llama FOM.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 43 / 48
Metodo GMRES
El FOM es equivalente al GC si A es SPD. Pero tiene inconvenientes:
FOM no trata de forma eficiente la memoria: Qk se tiene que guardar
completa. En cada iteraci
on,un nuevo vector base se tiene que calcular
y guardar. La ortogonalizaci
on ve siendo cada vez mas car con k.
FOM no tiene una propiedad de optimalidad.
El metodo es finito, pero esto s
olo tiene interes teorico ya que se
tendran que calcular n vectores base y guardarlos.
FOM no es robusto ya que Hk podra ser singular.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 44 / 48
Metodo GMRES
El metodo FOM para obtener xk es parte de una familia de tecnicas para extraer una
soluci
on aproximada a partir de un espacio de b
usqueda Qk haciendo que el residuo se
ortogonal a un espacio test Wk .
Esto se formula como: Sea xk = x0 + Qk yk . Encontrar yk tal que
WkT (r0 AQk yk ) = 0
Estas son las condiciones de Petrov-Galerkin.
Si Wk = Qk se llaman condicones de Galerkin.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 45 / 48
Metodo GMRES
Veamos un segundo metodo para obtener un mnimo del residuo.
Minimizar
g (xk ) = kA(xk x)k22 = rkT rk
es un problema bien definido incluso si A es no simetrica.
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 46 / 48
Metodo GMRES
El problema es encontrar xk = x0 + Qk yk tal que krk k is minimo.
rk = b Axk = r0 AQk yk = kr0 kq1 AQk yk
as
krk k = kkr0 kq1 AQk yk k
= kkr0 kQk+1 e1 Qk+1 H k yk k
= kkr0 ke1 H k yk k
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 47 / 48
Metodo GMRES
Resolviendo el problema sobredeterminado
H k yk = kr0 ke1
se tienen las iteraciones
xk = x0 + Qk yk
que minimizan el residuo. El algoritmo resultatnte se llama GMRES.
GMRES, es uno de los metodos mas populares para resolver sistemas no
simetricos
(UPV) M
etodos de gradiente. M
etodos de Krylov Curso 2012-2013 48 / 48