O método dos
gradientes conjugados
Francisco A. M. Gomes
MT402 – Matrizes – junho de 2008
Funções quadráticas

Chamamos de quadrática uma função
f(x):n  que pode ser escrita na
forma
f ( x )  x Ax  b x  c
1
2
T
T
onde Ann, bn e c.
Exemplo de função quadrática
3 2 
f ( x)  x 
x  2  8 x  c

2 6 
1
2
T
200
150
100
50
0
-50
5
6
0
4
2
-5
0
-2
-10
-4
Curvas de nível da quadrática

50
45
Repare que as curvas
de nível de f, neste
40
35
exemplo, são elipses.
30
25

20
15
Em cada curva do
gráfico, encontramos
10
5
5
10
15
20
25
30
35
40
45
50
pontos com o mesmo
valor de f.
Gradiente

O gradiente de uma função f(x) é dado
por

 x1 f ( x )
  f ( x )
x

f ' (x)   2
  
 

 x n f ( x )

Para um dado ponto x, o gradiente
fornece a direção de maior crescimento
de f(x).
Exemplo de gradiente



O gradiente é
ortogonal à curva de
nível.
No exemplo, as setas
apontam para o lado
de fora das elipses,
de modo que o
ponto (2,-2) é o
mínimo de f(x).
Em (2, -2), temos
f’(x) = 0.
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-4
-3
-2
-1
0
1
2
3
4
5
6
Quadráticas e sistemas lineares

De uma forma geral, escrevemos o gradiente
de uma função quadrática f(x) na forma:
f ' ( x )  21 A Tx  21 Ax  b

Se a matriz A é simétrica, temos
f ' ( x)  Ax  b


Podemos tentar minimizar uma função
quadrática f(x) igualando o gradiente a zero.
Isso equivale a resolver o sistema Ax = b.
Mínimo ou máximo?

Entretanto, dependendo de A, o ponto
crítico de f(x) pode um máximo local.
30
100
25
0
-100
20
-200
15
-300
10
-400
40
40
30
30
20
5
20
10
10
0
0
5
10
15
20
25
30
Pontos de sela

O ponto crítico de f(x) também pode ser
um ponto de sela.
30
400
25
300
200
20
100
0
15
-100
10
-200
40
30
40
30
20
5
20
10
10
0
0
5
10
15
20
25
30
Matrizes definidas positivas

8
6
4
2
0
-2
-4
-6
-8
-4
-3
-2
-1
0
1
2
3
4
5
6
Se a matriz A, além
de simétrica, é
definida positiva,
então a função
quadrática f(x) é um
parabolóide voltado
para cima, e a
solução de f’(x) = 0,
ou de Ax = b, é o
ponto de mínimo de
f(x).
O método da máxima descida

Também chamado de método do
gradiente.



Começa em um ponto x0 arbitrário
Gera uma seqüência de pontos x1, x2, ...
Caminha para o ponto de mínimo do
parabolóide, x*, seguindo a direção do
vetor –f´(x).
Definições importantes




Erro em xi: ei = x* – xi.
Resíduo em xi: ri = b – Axi.
O resíduo pode ser visto como o erro
transformado por A:
ri = b – Axi = A(x* – xi) = Aei.
O resíduo também pode ser visto como a
direção de máxima descida de f(x) em xi:
ri = –f’(xi).
Idéia do método


Dado x0, podemos encontrar um ponto
x1 que diminua o resíduo (ou f(x)),
dando um passo na direção de –f ’(x0).
Em outras palavras, podemos encontrar
x1 definido como
x1 = x0 + r0.

O parâmetro  é denominado
comprimento do passo.
Definindo o tamanho do passo

Para encontrar o melhor valor de ,
minimizamos f(x) ao longo da direção
definida por r0, ou seja, fazemos uma busca
linear:

 minimiza f(x0 + r0)
 é determinado igualando a zero a
derivada de f(x0 + r0) com relação a :
r0Tr0
 T
r0 Ar 0
Algoritmo
ri  b  Ax i
riTri
i  T
ri Ar i
x i 1  x i   i ri

Este procedimento é repetido por várias
iterações, até que xi esteja
suficientemente próximo de x*.
Taxa de convergência


Infelizmente, o método do gradiente não
possui, em geral, boa taxa de convergência.
Na verdade, podemos provar apenas que
f ( x i )  f ( x*)
   1


f ( x 0 )  f ( x*)    1

2i
onde  é o número de condição de A.
Assim, se A tem um número de condição
alto, pouco podemos esperar deste método.
Hiperelipsóides alongados


Quando (A) é alto, as curvas de nível
de f(x) são hiperelipsóides muito
alongados.
Para exemplificar isso, vamos resolver o
sistema Ax = b com
100  10 
A


10
1
.
001



19 
b 
 1
Neste caso, (A) = 1.0201e+005.
Hiperelipsóides alongados, parte 2

Para esse exemplo, o método da máxima
descida fez 3825 iterações, usando
||r||/||b||<0,00001 como critério de parada
e partindo do ponto [0, 0]T.
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
O método dos gradientes conjugados


É um método que usa a mesma fórmula
recursiva do método do gradiente, ou
seja,
xk+1 = xk + kdk,
Mas que não usa como vetor direção o
resíduo no ponto xk:
dk  rk = b – Axk.
Idéia geral do método

1.
2.
Vamos tentar encontrar, a cada iteração
k, um vetor dk que
seja linearmente independente dos
anteriores d0, ..., dk-1;
faça com que xk+1 minimize f(x) no
espaço gerado pelos vetores d já
definidos.
(Naturalmente, esse espaço é o
subespaço de Krylov K(A,d0,k))
Idéia geral do método, parte 2

Assim, xk+1 deve resolver o problema
mink 1 f ( xk  Dk w )
w 

onde Dk é a matriz cujas colunas são os
vetores d0, ..., dk.
Desta forma, ao final de n iterações,
teremos minimizado f(x) em n,
obtendo, assim, a solução do sistema
linear Ax = b.
Encontrando w

Para obter uma expressão computável para
xk+1, substituímos o termo
xk+1 = xk + Dk w
na fórmula de f(xk+1):
f ( x k 1 ) 

1
2
1
2
1
2
( x k  Dk w) T A ( x k  Dk w)  b T ( x k  Dk w)
T
T
x k Axk  x k ADk w 
w TDkT ADk w  b T x k  b TDk w
 f ( x k )  rkTDk w  12 w TDkT ADk w.
Encontrando w, parte 2

Derivando f(xk+1) com relação a w,
obtemos um ponto de mínimo dado por
1
w  (D AD k ) D r
T
k

T
k k
Deste modo,
x k 1  x k  D k w  x k  D k (D kT AD k )1D kTr k .
Vetores ortogonais

Constatamos que rk+1 é ortogonal aos
vetores di (colunas de Dk), pois
D kTr k 1
 D kTb  D kT Ax k 1
1
 D b  D A ( x k  D k (D AD k ) D r )
T
k
T
k k
T
k
T
k
T
k
 D r  D AD k (D kT AD k )1D kTr k
 0

Assim, djTrk+1 = 0, para j = 1, ..., k.
T
k k
Vetores ortogonais, parte 2


Se, nas iterações anteriores, os valores de xj,
j = 1, ..., k, também foram obtidos de modo a
minimizar f(xj-1 + Dj-1 wj-1), concluímos, por
indução, que
djTri = 0 para j < i.
Substituindo este termo na expressão de xk+1,
temos
x k 1  x k  Dk (DkT AD k ) 1 ek ,
onde  = dkTrk e ek corresponde à késima coluna
da identidade.
Usando vetores A-conjugados


Para simplificar a expressão de xk+1 e
facilitar os cálculos, podemos tentar fazer
com que a matriz DkTADk seja diagonal.
Para tanto, basta exigirmos que os
vetores dj sejam A-conjugados, ou seja,
que
diTAdj = 0 para i  j.
Como obter vetores A-conjugados

1.
2.
Lembrando que a iteração do método é
definida por x k 1  x k   k d k ,
podemos obter um conjunto de direções
A-conjugadas:
Escolhendo d0 = r0 (como no método da
máxima descida)
Definindo rk como uma combinação linear
de d0, ..., dk, de modo que
k 1
d k  r k   kj d j .
j 0
(1)
Como obter vetores A-conjugados (2)

Desta última equação, obtemos
k 1
diT Adk  diT Ark   kj diT Ad j .
j 0

Como as direções dj são A-conjugadas,
kj  d Tj Ark / d Tj Ad j

Mas rkT ri1  rkT ri   irkT Ad i, de modo que
rkT A di  (rkT ri  rkTri1 ) / i


Uma vez que, de (1), temos riTrj = 0, i  j,
Concluímos que kj = 0, j = 1,...,k – 2.
Como obter vetores A-conjugados (3)

Apenas k,k-1  0. Chamemos de k este valor

Claramente k  (1 / k 1 ) rkT rk / dkT1 A dk 1





Como k 1  dkT1rk 1 / dkT1 A dk 1
Obtemos k  rkT rk / dkT1rk 1
E como, de (1), temos dkT1rk 1  rkT1rk 1
Chegamos a
2
2
k  rk 2 / rk 1 2
Assim, a expressão de dk torna-se
d k  rk  k d k 1 .
Algoritmo

Dados de entrada:
1.
A matriz A (simétrica, definida positiva).
2.
O vetor b.
3.
4.
Uma aproximação inicial x0 da solução
do sistema.
Os parâmetros  (tolerância do resíduo)
e kmax (número máximo permitido de
iterações).
Algoritmo, parte 2

1
r0  b – Ax0

2
k  0; -1  0; d-1  0

3
Enquanto ||rk||2/||b||2 >  e k  kmax,

3.1
dk  rk + k dk-1

3.2
k  rkTrk/dkTAdk

3.3
xk+1  xk + kdk

3.4
rk+1  rk – kAdk

3.5
k  r k

3.6
kk+1
2
2
/ r k 1
2
2
Convergência do método


Teoricamente, o método converge para a
solução do sistema linear em n iterações.
Entretanto, nem sempre isso acontece em
virtude dos erros de arredondamento e
cancelamento que fazem com que
 o vetor resíduo perca precisão;
 os vetores direção deixem de ser Aconjugados, ou seja,
d iTAd j  0 , i  j

Isso ocorre quando A é mal condicionada.
Exemplo com A bem condicionada








» A = sprand(100,100,0.05,0.5)
+0.1*speye(100);
» A = A'*A;
» condest(A)
ans =
11.085
» b = rand(100,1);
» x = pcg(A,b);
pcg converged at iteration 15 to a solution with
relative residual 9.8e-007
Exemplo com A mal condicionada









» A = sprand(100,100,0.05,0.0001)
+0.1*speye(100);
» A = A'*A;
» condest(A)
ans =
3.7142e+008
» x = pcg(A,b);
pcg stopped at iteration 20 without converging
to the desired tolerance 1e-006...
» x = pcg(A,b,1e-5,1000);
pcg converged at iteration 232 to a solution
with relative residual 2.6e-006
Taxa de convergência
f (x i )  f (x*)    1 

 
f (x 0 )  f (x*)    1 


2i
Assim, a distância entre f(xi) e f(x*) está
limitada por um termo que é próximo de 1
se k(A) é grande.
Felizmente, o método costuma convergir
mais rápido do que podemos prever,
principalmente quando precondicionado.