Gradiente Conjugado

Gradiente conjugado
Miguel Vargas-Flix
[email protected]
http://www.cimat.mx/~miguelvargas
CIMAT, August 26, 2015
1/24
Mtodo de gradiente conjugado

Es un mtodo iterativo para minimizar funciones cuadrticas convexas f : n de la forma
1
f ( x )= x T A xxT b,
2
donde x , b n y A nn es una matriz simtrica positiva definida.
Es un mtodo de descenso de gradiente.
Ejemplo para n=2. Descenso de gradiente (izquierda), gradiente conjugado (derecha)

2/24
Para minimizar f ( x ) calculamos primero su gradiente,
f ( x )=
1 T
x A xxT b =A xb,
2
buscamos minimizar, por tanto igualamos a cero

A x=b,
es decir, buscamos resolver un sistema lineal de ecuaciones.
El mtodo de descenso de gradiente es:
Input: A , x 0, b,
k 0
repetir
r k bA x k
T
rk rk
k T
rk A rk
x k +1 x k + k r k
k k +1
hasta que r k<
x0 es una coordenada inicial.
3/24
Si A es simtrica positiva definida, entonces

T
x A x > 0 , para todo x0.
As, podemos a definir un producto interno
x , y A = x T A y.
Decimos que un vector x es conjugado a otro vector y con respecto a una matriz A si
x , y A =0, con x y.
La idea del algoritmo es utilizar direcciones conjugadas para el descenso en la bsqueda del punto
*
ptimo x , es decir
*
x =1 p1 + 2 p2 ++ n p n,
los coeficientes estn dados a partir de la combinacin lineal
*
A x =1 A p1 +2 A p 2 ++ n A pn =b.
A partir de una matriz A de rango n slo se pueden definir n vectores A-conjugados, por lo tanto
el algoritmo de gradiente conjugado garantiza la obtencin de una solucin en un mximo de n
iteraciones.
4/24
Definamos el residual r k como
r k =A xk b, k =1, 2,
la idea es lograr de forma iterativa que el residual tienda a ser cero, buscando cada vez mejores x k .
El procedimiento ser de forma iterativa, la forma de actualizacin ser
x k +1 =x k + p k ,
tomando p k como una direccin de descenso.
El tamao de paso k que minimiza la funcin f ( x ) a lo largo de la direccin x k + k p k es
T
r k pk
.
k = T
pk A pk
5/24
Si definimos p k +1 como la direccin ms cercana al residual r k bajo la restriccin de ser

conjugado.
Esta direccin est dada por la proyeccin de r k en el espacio ortogonal a p k con respecto al
producto interno inducido por A, as
T
pk A rk
p k +1 =r k + T
pk.
pk A p k
6/24
El algoritmo es el siguiente [Noce06 p108]:

Input: A , x 0, b,
r 0 A x0 b
p0 r 0
k 0
mientras r k>
T
r k pk
k T
pk A pk
x k +1 x k + k p k
r k +1 A x k +1b
T
r k +1 A pk
k +1
pTk A p k
p k +1 r k +1 +k +1 p k
k k +1
x0 es una coordenada inicial (puede ser igual a 0).
Generalmente no es necesario realizar las n iteraciones, se puede definir la precisin deseada
limitando la convergencia con una tolerancia .
7/24
El algoritmo de gradiente conjugado mejorado es el siguiente [Noce06 p112], con la variacin de

utilizar slo una multiplicacin matriz-vector:
Input: A , x 0, b,
r 0 A x0 b
p0 r 0
k 0
mientras r k>
w A pk
T
rk r k
k T
pk w
x k +1 x k + k p k
r k +1 r k + w
T
r k +1 r k +1
k +1
r Tk r k
p k +1 r k +1 +k +1 p k
k k +1
El proceso ms lento del algoritmo es la multiplicacin matriz-vector.
Se puede implementar eficientemente usando almacenamiento con compresin por renglones.
8/24
Compressed Row Storage

El mtodo Compressed Row Storage (compresin por renglones) [Saad03 p362], se guardan las
entradas no cero de cada rengln de A por separado.
8 4
8
0
2
0
0
4
0
0
9
0
0
1
1
3
0
0
3
0
0
0
0
0
7
0
0
0
0
0
1
5
1 2
1 3
3 4
2 1 7
1 3 5
9 3 1
V 3= { 2, 1, 7 }
J 3= {1, 3, 5 }
2 3 6
5
6
Con este mtodo, por cada rengln de la matriz se guardan dos arreglos para las entradas distintas
de cero de cada rengln. Un vector J i conteniendo los ndices y otro V i los valores, con i=1,, m.
Este esquema es adecuado cuando todos (o casi todos) los renglones tienen al menos una entrada
distinta de cero.
Este tipo de esquema se puede crear utilizando un vector de vectores.
9/24
Multiplicacin matriz vector

En la multiplicacin matriz-vector c=A b el rden de bsqueda es O ( 1 ), esto es porque no se hace
una bsqueda de las entradas del rengn, se toman las entradas una tras otra.
Sea J i el conjunto de ndices de las entradas no cero del rengln i de A.
|Ji| es el nmero de entradas no cero del rengln i de A. Ntese que |Vi|=|J i|.
(( )
c1
8 4
c2
0 0
c3
= 2 0
c4
0 9
c5
0 0
c6
0
1
1
3
0
0
3
0
0
0
ci = ai j b j
i =1
0
0
7
0
0
0
0
0
1
5
)( )
5
4
0
1
2
9
8 4
1 2
c1
c2
c3
c4
c5
c6
1 3
5
4
0
1
2
9
3 4
2 1 7
1 3 5
9 3 1
2 3 6
5
6
|J i|
ci = Vik b J
k =1
k
i
La ventaja de utilizar compresin por renglones es que los datos de cada rengln de la matriz de
rigidez son accesados en secuencia uno tras otro, esto producir una ventaja de acceso al entrar el
bloque de memoria de cada rengln en el cache del CPU.
10/24
El pseudocdigo es:
ConjugateGradient ( A , x , b , , step max )
A=( Vi , J i ), i=1, 2, , n CRS matrix
n
x Initial value
n
b Right side vector
Tolerance
step max Maximum number of steps
n
r Residual
n
p Descent direction
n
w Result of matrix*vector
j J i
k
sum sum+V i p j
w i sum
dot_pw dot_pw + p i w i
dot_rr 0
for i1, 2, , n
sum 0
for k 1, 2, ,V i
k
j J i
k
sum sum+V i x j
r i sumb i
pi r i
dot_rr dot_rr +r i r i
step 0
while ( step< stepmax ) and
dot_pw 0
for i1, 2, , n
sum 0
for k 1, 2, ,V i
dot_rr
dot_pw
new_dot_rr 0
for i1, 2, , n
x i x i + pi
r i r i + w i
new_dot_rr new_dot_rr+ r i r i
new_dot_rr
dot_rr
for i1, 2, , n
pi p i r i
( dot_rr> )
dot_rr new_dot_rr
step step+1
11/24
Nmero de condicin
El siguiente sistema de ecuaciones [Kind07], A x=b1
)( ) ( )
0.1 3.0 3.0 4.0

0.4 12.2 20.1 26.1
0.1 3.1 7.0 9.0
0.2 6.1 10.1 13.0
x1
10.1
x2
58.8
, tiene solucin
=
19.2
x3
29.4
x4
( )( )
x1
1.0
x2
1.0
,
=
1.0
x3
1.0
x4
b1
ntese que b12=69.227523428.

Si perturbamos el vector de la derecha un poco
0.1 3.0 3.0 4.0 x 1
10.10.005
10.095
0.4 12.2 20.1 26.1 x 2
58.8+0.005
58.805
, entonces
=
=
0.1 3.1 7.0 9.0 x 3
19.20.005
19.195
0.2 6.1 10.1 13.0 x 4
29.40.005
29.395
)( ) ( ) ( ) ( ) ( )
b1 +
x1
351.45
x2
11.00
,
=
1.05
x3
1.20
x4
b2
ahora b22 =69.227531373.

12/24
Si en vez de perturbar el vector de la derecha, perturbamos la matriz
0.100 3+0.005
3.000
4.000
0.400 12.200
20.100
26.100
0.100 3.100
70.005
9.000
0.200 6.100
10.100 130.005
)( ) ( ) ( ) ( )
x1
10.1
x2
58.8
, entonces
=
19.2
x3
29.4
x4
x1
3.425
x2
1.150
.
=
1.053
x3
0.957
x4
Un sistema de ecuaciones A x=b es considerado

bien condicionado si un pequeo cambio en los valores de A
o un pequeo cambio en b resulta en un pequeo cambio en x.
Un sistema de ecuaciones A x=b es considerado
mal condicionado si un pequeo cambio en los valores de A
o un pequeo cambio en b resulta en un cambio grande en x.
El nmero de condicin indica que tan sensible es una funcin a pequeos cambios en la entrada.
13/24
El nmero de condicin de una matriz A no singular, para una norma est dado por
( A )=AA1.
Para la norma 2,
max ( A )
( A )=AA =
,
min ( A )
1
donde son los valores singulares de la matriz.

Para una matriz A simtrica positiva definida,
max ( A )
( A )=
,
min ( A )
donde son los eigenvalores de A.
As, matrices con un nmero de condicin cercano a 1 se dicen que estn bien condicionadas.
Al reducir el nmero de condicin de un sistema, se puede acelerar la velocidad de
convergencia del gradiente conjugado.
14/24
Gradiente conjugado precondicionado

Entonces, en vez de resolver el problema
A xb=0,
se resuelve el problema
1
M ( A xb )=0,
1
con M
una matriz cuadrada, la cual recibe el nombre de precondicionador.

1
El mejor precondicionador sera claro M =A , as x=M b, y el gradiente conjugado

convergira en un paso.
1
Al igual que la matriz A, el precondicionador M
tiene que ser simtrico positivo definido.

1
Hay dos tipos de precondicionadores, implcitos M y explcitos M .

1
En algunos casos es costoso calcular M , en general se utilizan precondicionadores con inversas

fciles de calcular o precondicionadores implicitos que se puedan factorizar (con Cholesky, por
ejemplo M=L LT).
15/24
El algoritmo es el siguiente:
Input: A , x 0, b,
r 0 A x0 b
1
q0 M r0
p0 q 0
k 0
mientras r k>
w A pk
T
rk qk
k T
pk w
x k +1 x k + k p k
r k +1 r k + w
1
q k +1 M r k +1, or solve M q k + 1 r k +1
T
r k +1 q k +1
k
r Tk q k
p k +1 q k +1 +k +1 p k
k k +1
Ntese que ahora el algoritmo require aplicar el precondicionador en cada paso.
16/24
Gradiente conjugado + precondicionador Jacobi

El precondicionador Jacobi es el ms sencillo, consiste en hacer
M=diag ( A ),
de esta forma el precondicionador es una matriz diagonal, cuya inversa es fcil de calcular
1
( M1)i j = Ai i
0
si i= j
si i j
No se almacena todo M 1, lo usual es guardar un vector con slo los elementos de la diagonal.
El pseudocdigo es...
17/24
ConjugateGradientJacobi ( A , x , b , , step max )

A=( Vi , J i ), i=1, 2, , n CRS matrix
n
x Initial value
n
b Right side vector
Tolerance
step max Maximum number of steps
n
r Residual
p n Descent direction
n
M Preconditioner stored as vector
n
q Applied preconditioner
wn Result of matrix*vector
dot_rr 0
dot_rq0
for i1, 2, , n
sum 0
for k 1, 2, ,V i
k
j J i
k
sum sum+V i x j
r i sumb i
M i Ai i
qi r i / M i
pi q i
dot_rr dot_rr +r i r i
dot_rq dot_rq+r i qi
step 0
while ( step< stepmax ) and ( dot_rr> )
dot_pw 0
for i1, 2, , n
sum 0
for k 1, 2, ,V i
k
j J i
k
sum sum+V i p j
w i sum
dot_pw dot_pw + p i w i
dot_rq

dot_pw
new_dot_rr 0
new_dot_rq0
for i1, 2, , n
x i x i + pi
r i r i + w i
qi r i / M i
new_dot_rr new_dot_rr+ r i r i
new_dot_rqnew_dot_rq +r i q i
new_dot_rq

dot_rq
for i1, 2, , n
pi p i qi
dot_rr new_dot_rr
dot_rq new_dot_rq
step step+1
18/24
Por qu dan distinto resultado?

Es la misma suma, pero...
// sum1.cpp
// sum2.cpp
float a[] = {1.0, 1.0e-8, 2.0e-8,

3.0e-8, 4.0e-8, 5.0e-8,
1.0e-8, 2.0e-8, 3.0e-8,
4.0e-8, 5.0e-8};
float a[] = {1.0e-8, 2.0e-8, 3.0e-8,

4.0e-8, 5.0e-8, 1.0e-8,
2.0e-8, 3.0e-8, 4.0e-8,
5.0e-8, 1.0};
int main()
{
float sum = 0;
for (int i = 0; i < 11; ++i)
sum += a[i];
int main()
{
float sum = 0;
for (int i = 0; i < 11; ++i)
sum += a[i];
printf("Sum = %1.9f\n", sum);
return 0;
return 0;
}
Sum = 1.000000000
}
Sum = 1.000000358
Las solucin exacta es 1.000000300.

Al calcular los productos punto para el gradiente conjugado se tienen que sumar muchos nmeros
con punto flotante.
19/24
Hay tres solucines para evitar prdida de informacin.
1. Ordenar los nmeros antes de sumarlos

// sum1.cpp
#include <stdio.h>
// sum3.cpp
#include <stdio.h>
float a[] = {1.0, 1.0e-8, 2.0e-8,

3.0e-8, 4.0e-8, 5.0e-8,
1.0e-8, 2.0e-8, 3.0e-8,
4.0e-8, 5.0e-8};
float a[] = {1.0e-8, 1.0e-8, 2.0e-8,

2.0e-8, 3.0e-8, 3.0e-8,
4.0e-8, 4.0e-8, 5.0e-8,
5.0e-8, 1.0};
int main()
{
float sum = 0;
for (int i = 0; i < 11; ++i)
sum += a[i];
int main()
{
float sum = 0;
for (int i = 0; i < 11; ++i)
sum += a[i];
return 0;
return 0;
}
Sum = 1.000000000
}
Sum = 1.000000358
20/24
2. Usar tipos de punto flotante ms grandes para acumular la suma

// sum1.cpp
#include <stdio.h>
// sum4.cpp
#include <stdio.h>
float a[] = {1.0, 1.0e-8, 2.0e-8,

3.0e-8, 4.0e-8, 5.0e-8,
1.0e-8, 2.0e-8, 3.0e-8,
4.0e-8, 5.0e-8};
float a[] = {1.0, 1.0e-8, 2.0e-8,

3.0e-8, 4.0e-8, 5.0e-8,
1.0e-8, 2.0e-8, 3.0e-8,
4.0e-8, 5.0e-8};
int main()
{
float sum = 0;
for (int i = 0; i < 11; ++i)
{
sum += a[i];
}
int main()
{
double sum = 0;
for (int i = 0; i < 11; ++i)
{
sum += a[i];
}
return 0;
return 0;
}
Sum = 1.000000000
}
Sum = 1.000000300
Cuando se suman doubles, se puede acumular el resultado utilizando long double.

Qu pasa cuando no hay puntos flotantes de mayor capacidad?
21/24
3. El algoritmo de sumacin de Kahan

Este algoritmo [Gold91] evita perder la contribucin de los nmeros pequeos.
// sum4.cpp
float a[] = {1.0, 1.0e-8, 2.0e-8,
3.0e-8, 4.0e-8, 5.0e-8,
1.0e-8, 2.0e-8, 3.0e-8,
4.0e-8, 5.0e-8};
int main()
{
float sum = 0;
float c = 0;
for (int i = 0; i < 11; ++i)
{
float y = a[i] + c;
float t = sum + y;
c = y - (t - sum);
sum = t;
}
c guardar la compensacin para los valores pequeos

suma la compensacin al valor leido
si sum es grande, los valores de y se pierden
(t - sum) es la parte grande de y, y - (t - sum) es la chica

return 0;
}
1.000000358
http://en.wikipedia.org/wiki/Kahan_summation_algorithm
22/24
Preguntas?
[email protected]
23/24
Referencias
[Noce06] J. Nocedal, S. J. Wright. Numerical Optimization. Springer, 2006.
[Piss84] S. Pissanetzky. Sparse Matrix Technology. Academic Press, 1984.
[Saad03] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003.
[Kind07] U. Kindeln. Resolucin de sistemas lineales de ecuaciones: Mtodo del gradiente
conjugado. Universidad Politcnica de Madrid. 2007
[Gold91] D. Goldberg. What Every Computer Scientist Should Know About Floating-Point
Arithmetic. Computing Surveys. Association for Computing Machinery. 1991.
24/24

Gradiente Conjugado

Cargado por

Copyright:

Formatos disponibles

Gradiente Conjugado

Cargado por

Información del documento

Descripción original:

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Gradiente Conjugado

Cargado por

Copyright:

Formatos disponibles

Gradiente conjugado

Mtodo de gradiente conjugado

Ejemplo para n=2. Descenso de gradiente (izquierda), gradiente conjugado (derecha)

Para minimizar f ( x ) calculamos primero su gradiente,

buscamos minimizar, por tanto igualamos a cero

Si A es simtrica positiva definida, entonces

x A x > 0 , para todo x0.

As, podemos a definir un producto interno

Definamos el residual r k como

Si definimos p k +1 como la direccin ms cercana al residual r k bajo la restriccin de ser

El algoritmo es el siguiente [Noce06 p108]:

El algoritmo de gradiente conjugado mejorado es el siguiente [Noce06 p112], con la variacin de

Compressed Row Storage

Multiplicacin matriz vector

0.1 3.0 3.0 4.0

ntese que b12=69.227523428.

ahora b22 =69.227531373.

Si en vez de perturbar el vector de la derecha, perturbamos la matriz

Un sistema de ecuaciones A x=b es considerado

donde son los valores singulares de la matriz.

Gradiente conjugado precondicionado

una matriz cuadrada, la cual recibe el nombre de precondicionador.

El mejor precondicionador sera claro M =A , as x=M b, y el gradiente conjugado

Al igual que la matriz A, el precondicionador M

tiene que ser simtrico positivo definido.

Hay dos tipos de precondicionadores, implcitos M y explcitos M .

En algunos casos es costoso calcular M , en general se utilizan precondicionadores con inversas

Gradiente conjugado + precondicionador Jacobi

ConjugateGradientJacobi ( A , x , b , , step max )

Por qu dan distinto resultado?

float a[] = {1.0, 1.0e-8, 2.0e-8,

float a[] = {1.0e-8, 2.0e-8, 3.0e-8,

printf("Sum = %1.9f\n", sum);

printf("Sum = %1.9f\n", sum);

Las solucin exacta es 1.000000300.

Hay tres solucines para evitar prdida de informacin.

1. Ordenar los nmeros antes de sumarlos

float a[] = {1.0, 1.0e-8, 2.0e-8,

float a[] = {1.0e-8, 1.0e-8, 2.0e-8,

printf("Sum = %1.9f\n", sum);

printf("Sum = %1.9f\n", sum);

2. Usar tipos de punto flotante ms grandes para acumular la suma

float a[] = {1.0, 1.0e-8, 2.0e-8,

float a[] = {1.0, 1.0e-8, 2.0e-8,

printf("Sum = %1.9f\n", sum);

printf("Sum = %1.9f\n", sum);

Cuando se suman doubles, se puede acumular el resultado utilizando long double.

3. El algoritmo de sumacin de Kahan

c guardar la compensacin para los valores pequeos

printf("Sum = %1.9f\n", sum);

También podría gustarte