Método de Gradiente Conjugado

INTRODUCCIÓN

Si A es una matriz simétrica y definida positiva, se cumple que

1
min 𝜙(𝑥) = 2 𝑥 𝑇 𝐴𝑥 − 𝑥 𝑇 𝑏

Es equivalente a resolver
𝐴𝑥 = 𝑏

Se sabe que una función de dos o más variables tiene un mínimo local en aquellos
puntos donde se cumplen las siguientes condiciones:

∇ 𝜙(𝑥) = 0

det⁡(∇2 𝜙 𝑥 ) > 0
INTRODUCCIÓN

Para
1
𝜙(𝑥) = 2 𝑥 𝑇 𝐴𝑥 − 𝑥 𝑇 𝑏

Se tiene que
∇ 𝜙(𝑥) = 𝐴𝑥 − 𝑏

∇2 𝜙 𝑥 = 𝐴

Siempre y cuando A sea simétrica. Entonces, la condición ∇ 𝜙(𝑥) = 0 resulta

equivalente a encontrar la solución del sistema.
Considerando el sistema lineal

𝑎11 𝑎12 𝑥1 𝑏1
𝑎21 𝑎22 𝑥2 =
𝑏2

1 𝑇 1 𝑎11 𝑎12 𝑥1 𝑏1
𝜙 𝑥 = 𝑥 𝐴𝑥 − 𝑥 𝑇 𝑏 = 𝑥1 𝑥2
𝑎21
𝑥
𝑎22 𝑥2 − 1
𝑥2
2 2 𝑏2

1
𝜙 𝑥 = 𝑎11 𝑥12 + 𝑎12 𝑥1 𝑥2 + 𝑎21 𝑥1 𝑥2 + 𝑎22 𝑥22 − 𝑥1 𝑏1 − 𝑥2 𝑏2
2

𝜕𝜙(𝑥)
𝜕𝑥1
∇ 𝜙(𝑥) =
𝜕𝜙(𝑥)
𝜕𝑥2

𝜕𝜙(𝑥)
= 𝑎11 𝑥1 + 12𝑎12 𝑥2 + 12𝑎21 𝑥2 − 𝑏1
𝜕𝑥1

𝜕𝜙(𝑥) 1
= 2𝑎12 𝑥1 + 12𝑎21 𝑥1 + 𝑎22 𝑥2 − 𝑏2
𝜕𝑥2
Si la matriz es simétrica, 𝑎12 = 𝑎21 y se tiene

𝜕𝜙(𝑥)
= 𝑎11 𝑥1 + 𝑎12 𝑥2 − 𝑏1
𝜕𝑥1

𝜕𝜙(𝑥)
= 𝑎21 𝑥1 + 𝑎22 𝑥2 − 𝑏2
𝜕𝑥2

𝜕𝜙(𝑥)
𝜕𝑥1 𝑎11 𝑎12 𝑥1 𝑏1
= 𝑎 𝑎22 𝑥2 −
𝜕𝜙(𝑥) 21 𝑏2
𝜕𝑥2

Por tanto,

∇ 𝜙(𝑥) = 𝐴𝑥 − 𝑏
 𝜕 2 𝜙(𝑥) 𝜕
2 = 𝜕𝑥 (𝑎11 𝑥1 + 𝑎12 𝑥2 − 𝑏1 ) ⁡ = 𝑎11
𝜕𝑥1 1

𝜕 2 𝜙(𝑥) 𝜕
= (𝑎 𝑥 + 𝑎12 𝑥2 − 𝑏1 ) ⁡ = 𝑎12
𝜕𝑥1 𝜕𝑥2 𝜕𝑥2 11 1

𝜕 2 𝜙(𝑥) 𝜕
= (𝑎 𝑥 + 𝑎22 𝑥2 − 𝑏1 ) ⁡ = 𝑎21
𝜕𝑥2 𝜕𝑥1 𝜕𝑥1 21 1

𝜕 2 𝜙(𝑥) 𝜕
= (𝑎 𝑥 + 𝑎22 𝑥2 − 𝑏1 ) ⁡ = 𝑎22
𝜕𝑥12 𝜕𝑥2 21 1

𝜕 2 𝜙(𝑥) 𝜕 2 𝜙(𝑥)
𝜕𝑥12 𝜕𝑥1 𝜕𝑥2 𝑎11 𝑎12
∇2 𝜙 𝑥 = 2 2 = 𝑎 𝑎22
𝜕 𝜙(𝑥) 𝜕 𝜙(𝑥) 21

𝜕𝑥2 𝜕𝑥1 𝜕𝑥12

Entonces
∇2 𝜙 𝑥 = 𝐴
MÉTODO DE DESCENSO SIMPLE

Para el cálculo del mínimo de la función podemos utilizar la ecuación iterativa

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘

Donde 𝛼𝑘 es el tamaño de paso (o distancia) a recorrer en la dirección 𝑝𝑘 .

MÉTODO DE DESCENSO SIMPLE

En la ecuación
𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘

Dada la dirección 𝑝𝑘 , es posible calcular 𝛼𝑘 tal que el punto obtenido corresponda a

un mínimo local en esa dirección.

Expresamos 𝜙(𝑥) como una función de 𝛼 considerando que 𝑥 𝑘 y 𝑝𝑘 son constantes

𝑓(𝛼) = 𝜙(𝑥 𝑘 + 𝛼𝑝𝑘 )

MÉTODO DE DESCENSO SIMPLE

−𝑝𝑘𝑇 𝛻 𝜙(𝑥 𝑘 ) −𝑝𝑘𝑇 (𝐴𝑥 𝑘 − 𝑏)

𝛼𝑘 = =
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘
MÉTODO DE DESCENSO SIMPLE

𝑓(𝛼) = 𝜙(𝑥 𝑘 + 𝛼𝑝𝑘 )

El mínimo de 𝑓 se encuentra en el valor de 𝛼⁡tal que 𝑓′ 𝛼 = 0.

𝑓′ 𝛼 = 𝐷𝑝𝑘 𝜙(𝑥 𝑘 + 𝛼𝑝𝑘 )

= 𝑝𝑘𝑇 ∇𝜙 𝑥 𝑘 + 𝛼𝑝𝑘 = ⁡ 𝑝𝑘𝑇 𝐴 𝑥 𝑘 + 𝛼𝑝𝑘 − 𝑏

𝑓 ′ 𝛼 = 𝑝𝑘𝑇 𝐴𝑥 𝑘 + ⁡𝛼𝑝𝑘𝑇 𝐴𝑝𝑘 − 𝑝𝑘𝑇 𝑏

MÉTODO DE DESCENSO SIMPLE

Igualando a cero esta última ecuación

𝑓 ′ 𝛼 = 𝑝𝑘𝑇 𝐴𝑥 𝑘 + ⁡𝛼𝑝𝑘𝑇 𝐴𝑝𝑘 − 𝑝𝑘𝑇 𝑏 = 0

−𝑝𝑘𝑇 𝐴𝑥 𝑘 + 𝑝𝑘𝑇 𝑏 −𝑝𝑘𝑇 (𝐴𝑥 𝑘 − 𝑏) −𝑝𝑘𝑇 ∇ 𝜙(𝑥 𝑘 )

𝛼= = =⁡
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘

Por lo que se obtiene

−𝑝𝑘𝑇 ∇ 𝜙(𝑥 𝑘 ) −𝑝𝑘𝑇 (𝐴𝑥 𝑘 − 𝑏)

𝛼𝑘 = =
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘
MÉTODO DE DESCENSO SIMPLE

El método hasta ahora está dado por

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘

−𝑝𝑘𝑇 ∇ 𝜙(𝑥 𝑘 ) −𝑝𝑘𝑇 (𝐴𝑥 𝑘 − 𝑏)

𝛼𝑘 = =
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘

Las direcciones 𝑝𝑘 son arbitrarias. Si se eligen en la dirección contraria al gradiente de

la función, se tiene el método de descenso de gradiente.

Definimos al gradiente de la función como el vector residuo

𝑟𝑘 = 𝐴𝑥 𝑘 − 𝑏
MÉTODO DE DESCENSO DE GRADIENTE

Algoritmo

For(k=0 to MAXIT, step 1)

𝑟𝑘 = 𝐴𝑥 𝑘 − 𝑏
IF( 𝑟𝑘 <e1) → Solución
𝑟𝑘𝑇 𝑟𝑘
𝛼𝑘 =
𝑟𝑘𝑇 𝐴𝑟𝑘

𝑥 𝑘+1 = 𝑥 𝑘 − 𝛼𝑘 𝑟𝑘
End
→ Número máximo de iteraciones alcanzado
MÉTODO DE DIRECCIONES CONJUGADAS

Otra forma muy conveniente para este problema en particular, es utilizar un conjunto
especial de direcciones llamadas direcciones conjugadas, las cuales cumplen con la
siguiente propiedad

= 0⁡si⁡𝑖 ≠ 𝑗
𝑝𝑖𝑇 𝐴𝑝𝑗 →
≠ 0⁡si⁡𝑖 = 𝑗

Las direcciones conjugadas son un conjunto de n vectores linealmente independientes

entre sí, por lo que forman una base vectorial.
MÉTODO DE DIRECCIONES CONJUGADAS

Si las direcciones de búsqueda corresponden a las direcciones conjugadas, el método

convergerá en n iteraciones sin importar cuál sea el vector inicial 𝑥 0 .

Demostración:
𝑥 ∗ − 𝑥 0 = ⁡ 𝜎0 𝑝0 + 𝜎1 𝑝1 + 𝜎2 𝑝2 + ⋯ + 𝜎𝑛−1 𝑝𝑛−1

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 0 = ⁡ 𝜎0 𝑝𝑘𝑇 𝐴𝑝0 + 𝜎1 𝑝𝑘𝑇 𝐴𝑝1 + 𝜎2 𝑝𝑘𝑇 𝐴𝑝2 + ⋯ + 𝜎𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘 + ⋯ + 𝜎𝑛−1 𝑝𝑘𝑇 𝐴𝑝𝑛−1

Por la propiedad de conjugabilidad

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 0 = ⁡ 𝜎𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 0
𝜎𝑘 =
𝑝𝑘𝑇 𝐴𝑝𝑘
MÉTODO DE DIRECCIONES CONJUGADAS

𝑥 𝑘 = 𝑥 0 + ⁡ 𝛼0 𝑝0 + 𝛼1 𝑝1 + 𝛼2 𝑝2 + ⋯ + 𝛼𝑘−1 𝑝𝑘−1

𝑝𝑘𝑇 𝐴𝑥 𝑘 = 𝑝𝑘𝑇 𝐴𝑥 0 + ⁡ 𝛼0 𝑝𝑘𝑇 𝐴𝑝0 + 𝛼1 𝑝𝑘𝑇 𝐴𝑝1 + 𝛼2 𝑝𝑘𝑇 𝐴𝑝2 + ⋯ + 𝛼𝑘−1 𝑝𝑘𝑇 𝐴𝑝𝑘−1

𝑝𝑘𝑇 𝐴𝑥 𝑘 = 𝑝𝑘𝑇 𝐴𝑥 0

𝑝𝑘𝑇 𝐴(𝑥 𝑘 −𝑥 0 ) = 0

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 0 = 𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 𝑘 + 𝑥 𝑘 − 𝑥 0 = 𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 𝑘 ) + 𝑝𝑘𝑇 𝐴(𝑥 𝑘 −𝑥 0 = 𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 𝑘

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 0 = 𝑝𝑘𝑇 𝐴𝑥 ∗ − 𝐴𝑥 𝑘 = 𝑝𝑘𝑇 𝑏 − 𝐴𝑥 𝑘

𝑝𝑘𝑇 𝐴 𝑥 ∗ − 𝑥 𝑘 −𝑝𝑘𝑇 𝐴𝑥 𝑘 − 𝑏
𝜎𝑘 = = = 𝛼𝑘
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘
MÉTODO DE DIRECCIONES CONJUGADAS

Finalmente
𝑥 ∗ − 𝑥 0 = ⁡ 𝛼0 𝑝0 + 𝛼1 𝑝1 + 𝛼2 𝑝2 + ⋯ + 𝛼𝑛−1 𝑝𝑛−1

𝑥 ∗ = 𝑥 0 + ⁡ 𝛼0 𝑝0 + 𝛼1 𝑝1 + 𝛼2 𝑝2 + ⋯ + 𝛼𝑛−1 𝑝𝑛−1

Que resulta ser la secuencia dada por

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘
MÉTODO DE DIRECCIONES CONJUGADAS

Las direcciones conjugadas están relacionadas con los eigenvectores de la matriz del
sistema. Para un sistema de 2x2 tienen la siguiente interpretación geométrica:
MÉTODO DE DIRECCIONES CONJUGADAS

Las direcciones conjugadas están relacionadas con los eigenvectores de la matriz del
sistema. Para un sistema de 2x2 tienen la siguiente interpretación geométrica:
MÉTODO DE DIRECCIONES CONJUGADAS

Geométricamente, el algoritmo de direcciones conjugadas se representa por

MÉTODO DE DIRECCIONES CONJUGADAS

Algoritmo

For(k=0 to n-1, step 1)

𝑟𝑘 = 𝐴𝑥 𝑘 − 𝑏
IF( 𝑟𝑘 <e1) → Solución
−𝑝𝑘𝑇 𝑟𝑘
𝛼𝑘 =
𝑝𝑘𝑇 𝐴𝑝𝑘

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘
End
→ Número máximo de iteraciones alcanzado
MÉTODO DE GRADIENTE CONJUGADO

El cálculo directo de las direcciones conjugadas es un problema muy costoso

computacionalmente. Sin embargo, no es necesario hacer esto, puesto que se cumple
el siguiente teorema:

𝑝𝑖𝑇 𝑟𝑘 = 0 para 𝑖 = 0,1,2, ⋯ , 𝑘 − 1

Esta propiedad permite calcular la dirección 𝑝𝑘 como una combinación lineal entre el
gradiente de la función 𝑟𝑘 (o residuo) y la dirección conjugada anterior 𝑝𝑘−1 , lo que
permite ahorrar muchas operaciones y memoria. La primera dirección 𝑝0 se escoge
igual a la dirección contraria al gradiente, por lo que 𝑝0 = −𝑟0

𝑝𝑘 = −𝑟𝑘 + 𝛽𝑘 𝑝𝑘−1

𝑝𝑘+1 = −𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘

MÉTODO DE GRADIENTE CONJUGADO

𝑝𝑘+1 = −𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘

𝑝𝑘𝑇 𝐴(𝑝𝑘+1 = −𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘 )

𝑝𝑘𝑇 𝐴𝑝𝑘+1 = −𝑝𝑘𝑇 𝐴𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘𝑇 𝐴𝑝𝑘

0 = −𝑝𝑘𝑇 𝐴𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘𝑇 𝐴𝑝𝑘

𝑝𝑘𝑇 𝐴𝑟𝑘+1
𝛽𝑘+1 = 𝑇
𝑝𝑘 𝐴𝑝𝑘
MÉTODO DE GRADIENTE CONJUGADO

Algoritmo original

𝑟0 = 𝐴𝑥 0 − 𝑏
𝑝0 = −𝑟0
For(k=0 to MAXIT, step 1)
−𝑝𝑘𝑇 𝑟𝑘
𝛼𝑘 =
𝑝𝑘𝑇 𝐴𝑝𝑘
𝑘+1 𝑘
𝑥 = 𝑥 + 𝛼𝑘 𝑝𝑘
𝑟𝑘+1 = 𝐴𝑥 𝑘+1 − 𝑏
IF( 𝑟𝑘+1 <e1) → Solución
𝑝𝑘𝑇 𝐴𝑟𝑘+1
𝛽𝑘+1 =
𝑝𝑘𝑇 𝐴𝑝𝑘
𝑝𝑘+1 = −𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘
End
→ Número máximo de iteraciones alcanzado
MÉTODO DE GRADIENTE CONJUGADO

El algoritmo presentado requiere 3 productos matriz-vector por iteración, teniendo

una complejidad de orden n2. Mediante algunas sustituciones algebraicas, es posible
obtener una versión del algoritmo que emplea un solo producto matriz-vector,
resultado en la forma práctica del algoritmo.

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼𝑘 𝑝𝑘

𝑟𝑘+1 = 𝐴𝑥 𝑘+1 − 𝑏 = 𝐴(𝑥 𝑘 + 𝛼𝑘 𝑝𝑘 ) − 𝑏

= 𝐴𝑥 𝑘 − 𝑏 + 𝛼𝑘 𝐴𝑝𝑘

𝑟𝑘+1 = 𝑟𝑘 + 𝛼𝑘 𝐴𝑝𝑘
MÉTODO DE GRADIENTE CONJUGADO

Se cumple también que 𝑟𝑖𝑇 𝑟𝑘 = 0 para 𝑖 = 0,1,2, ⋯ , 𝑘 − 1 . Esto permite realizar las
siguientes simplificaciones

−𝑝𝑘𝑇 𝑟𝑘 −𝑟𝑘𝑇 𝑝𝑘 −𝑟𝑘𝑇 (−𝑟𝑘 +𝛽𝑘 𝑝𝑘−1 ) 𝑟𝑘𝑇 𝑟𝑘

1. 𝛼𝑘 = = = =
𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘

2. Usando 𝑟𝑘+1 = 𝑟𝑘 + 𝛼𝑘 𝐴𝑝𝑘 → ⁡𝛼𝑘 𝐴𝑝𝑘 = 𝑟𝑘+1 − 𝑟𝑘

𝑝𝑘𝑇 𝐴𝑟𝑘+1 𝑟𝑘+1

𝑇 𝑇
𝐴𝑝𝑘 𝑟𝑘+1 𝑇
𝐴𝑝𝑘 𝛼𝑘 𝑟𝑘+1 (𝑟𝑘+1 − 𝑟𝑘 )
𝛽𝑘+1 = 𝑇 = 𝑇 = 𝑇 ⋅ = 𝑇
𝑝𝑘 𝐴𝑝𝑘 𝑝𝑘 𝐴𝑝𝑘 𝑝𝑘 𝐴𝑝𝑘 𝛼𝑘 𝑝𝑘 (𝑟𝑘+1 − 𝑟𝑘 )⁡

𝑇 𝑇 𝑇
𝑟𝑘+1 𝑟𝑘+1 𝑟𝑘+1 𝑟𝑘+1 𝑟𝑘+1 𝑟𝑘+1
𝛽𝑘+1 = = =
−𝑝𝑘𝑇 𝑟𝑘 − −𝑟𝑘 + 𝛽𝑘 𝑝𝑘−1 𝑇 𝑟𝑘 𝑟𝑘𝑇 𝑟𝑘
MÉTODO DE GRADIENTE CONJUGADO

Las simplificaciones resultantes son las siguientes:

𝑟𝑘+1 = 𝐴𝑥 𝑘+1 − 𝑏 = 𝑟𝑘 + 𝛼𝑘 𝐴𝑝𝑘

−𝑝𝑘𝑇 𝑟𝑘 𝑟𝑘𝑇 𝑟𝑘
𝛼𝑘 = 𝑇 =
𝑝𝑘 𝐴𝑝𝑘 𝑝𝑘𝑇 𝐴𝑝𝑘

𝑝𝑘𝑇 𝐴𝑟𝑘+1 𝑟𝑘+1

𝑇
𝑟𝑘+1
𝛽𝑘+1 = 𝑇 =
𝑝𝑘 𝐴𝑝𝑘 𝑟𝑘𝑇 𝑟𝑘
MÉTODO DE GRADIENTE CONJUGADO

Algoritmo práctico

𝑟0 = 𝐴𝑥 0 − 𝑏
𝑝0 = −𝑟0
For(k=0 to MAXIT, step 1)
𝑟𝑘𝑇 𝑟𝑘
𝛼𝑘 =
𝑝𝑘𝑇 𝐴𝑝𝑘
𝑘+1 𝑘
𝑥 = 𝑥 + 𝛼𝑘 𝑝𝑘
𝑟𝑘+1 = 𝑟𝑘 + 𝛼𝑘 𝐴𝑝𝑘 ⁡
IF( 𝑟𝑘+1 <e1) → Solución
𝑇
𝑟𝑘+1 𝑟𝑘+1
𝛽𝑘+1 =
𝑟𝑘𝑇 𝑟𝑘
𝑝𝑘+1 = −𝑟𝑘+1 + 𝛽𝑘+1 𝑝𝑘
End
→ Número máximo de iteraciones alcanzado
NOTAS SOBRE GRADIENTE CONJUGADO

1. Si no existieran errores numéricos, GC convergería cuando mucho en n iteraciones.

Ya que esto no es posible, hay casos que tomarán más iteraciones en converger
(recomendable usar 1.5n como límite).

2. El método funciona para matrices simétricas, sean éstas o no definidas positivas.

3. La convergencia depende de la distribución de los eigenvalores de la matriz

(asociados a las direcciones conjugadas). Este concepto está relacionado con el
número de condición, por lo que si la matriz está malcondicionada, el método
difícilmente converge.

4. Es fácil de programar y paralelizar.