O método dos gradientes conjugados

Francisco A. M. GomesMT402 – Matrizes – junho de 2008

Funções quadráticas

Chamamos de quadrática uma função f(x):n que pode ser escrita na forma

onde Ann, bn e c.

cxbAxx)x(f TT21

Exemplo de função quadrática

-4-2

02

46

-10

-5

0

5-50

0

50

100

150

200

cx82x62

23x)x(f T

21

Curvas de nível da quadrática

5 10 15 20 25 30 35 40 45 50

5

10

15

20

25

30

35

40

45

50 Repare que as

curvas de nível de

f, neste exemplo,

são elipses.

Em cada curva do

gráfico,

encontramos

pontos com o

mesmo valor de f.

Gradiente

O gradiente de uma função f(x) é dado por

Para um dado ponto x, o gradiente fornece a direção de maior crescimento de f(x).

)x(f

)x(f

)x(f

)x('f

n

2

1

x

x

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-6

-5

-4

-3

-2

-1

0

1

2

3

4

Quadráticas e sistemas lineares

De uma forma geral, escrevemos o gradiente de uma função quadrática f(x) na forma:

Se a matriz A é simétrica, temos

Podemos tentar minimizar uma função quadrática f(x) igualando o gradiente a zero.

Isso equivale a resolver o sistema Ax = b.

bAxxA)x('f 21T

21

bAx)x('f

Mínimo ou máximo?

Entretanto, dependendo de A, o ponto crítico de f(x) pode um máximo local.

010

2030

40

0

10

20

30

40-400

-300

-200

-100

0

100

5 10 15 20 25 30

5

10

15

20

25

30

Pontos de sela

O ponto crítico de f(x) também pode ser um ponto de sela.

010

2030

40

0

10

20

30

40-200

-100

0

100

200

300

400

5 10 15 20 25 30

5

10

15

20

25

30

Matrizes definidas positivas

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).

-4 -3 -2 -1 0 1 2 3 4 5 6-8

-6

-4

-2

0

2

4

6

8

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 :

0T0

0T0

Arrrr

Algoritmo

Este procedimento é repetido por várias iterações, até que xi esteja suficientemente próximo de x*.

iii1i

iTi

iTi

i

ii

rxx

Arrrr

Axbr

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

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.

i2

0

i

11

*)x(f)x(f*)x(f)x(f

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

Neste caso, (A) = 1.0201e+005.

1

19b

001.110

10100A

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 0.2 0.4 0.6 0.8 1

0

0.1

0.2

0.3

0.4

0.5

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

Vamos tentar encontrar, a cada iteração k, um vetor dk que

1. seja linearmente independente dos anteriores d0, ..., dk-1;

2. 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

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.

)wDx(fmin kkw 1k

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):

.wADDwwDr)x(f

wDbxbwADDw

wADxAxx

)wDx(b)wDx(A)wDx()x(f

kTk

T21

kTkk

kT

kT

kTk

T21

kT

kkT

k21

kkT

kkT

kk21

1k

Encontrando w, parte 2

Derivando f(xk+1) com relação a w,

obtemos um ponto de mínimo dado por

Deste modo,

kTk

1k

Tk rD)ADD(w

.rD)ADD(DxwDxx kTk

1k

Tkkkkk1k

Vetores ortogonaisVetores ortogonais

Constatamos que rk+1 é ortogonal aos vetores di (colunas de Dk), pois

Assim, djTrk+1 = 0, para j = 1, ..., k.

0

rD)ADD(ADDrD

)rD)ADD(Dx(ADbD

AxDbDrD

kTk

1k

Tkk

Tkk

Tk

kTk

1k

Tkkk

Tk

Tk

1kTk

Tk1k

Tk

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

onde = dkTrk e ek corresponde à késima

coluna da identidade.

,e)ADD(Dxx k1

kTkkk1k

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

Lembrando que a iteração do método é definida por podemos obter um conjunto de direções A-conjugadas:

1. Escolhendo d0 = r0 (como no método da máxima descida)

2. Definindo rk como uma combinação linear de d0, ..., dk, de modo que

(1)

,dxx kkk1k

1k

0jjkjkk .drd

Como obter vetores A-conjugados (2)

Desta última equação, obtemos

Como as direções dj são A-conjugadas,

Mas , de modo que

Uma vez que, de (1), temos riTrj = 0, i j,

Concluímos que kj = 0, j = 1,...,k – 2.

1k

0jj

Tikjk

Tik

Ti .AddArdAdd

jTjk

Tjkj Add/Ard

iTkii

Tk1i

Tk Adrrrrr

i1iTki

Tki

Tk /)rrrr(Adr

Como obter vetores A-conjugados (3)

Apenas k,k-1 0. Chamemos de k este

valor Claramente

Como Obtemos E como, de (1), temos Chegamos a

Assim, a expressão de dk torna-se

1kT

1kkTk1kk Add/rr)/1(

2

21k

2

2kk r/r

.drd 1kkkk

1kT

1k1kT

1k1k Add/rd

1kT

1kkTkk rd/rr

1kT

1k1kT

1k rrrd

Algoritmo

Dados de entrada:

1. A matriz A (simétrica, definida positiva).

2. O vetor b.

3. Uma aproximação inicial x0 da

solução do sistema.

4. 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/dk

TAdk

3.3 xk+1 xk + kdk

3.4 rk+1 rk – kAdk

3.5 3.6 k k + 1

2

21k

2

2kk r/r

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 A-

conjugados, ou seja,

Isso ocorre quando A é mal condicionada.

ji,0Add jTi

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

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.

i2

0

i

11

*)x(f)x(f*)x(f)x(f