O método dos gradientes conjugados

1. O método dos gradientes conjugados Francisco A. M. Gomes MT402 – Matrizes – junho de 2008

2. Funções quadráticas • Chamamos de quadrática uma função f(x):n  que pode ser escrita na forma onde Ann, bn e c.

3. Exemplo de função quadrática

4. 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. Curvas de nível da quadrática

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

6. 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. Exemplo de gradiente

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

8. Mínimo ou máximo? • Entretanto, dependendo de A, o ponto crítico de f(x) pode um máximo local.

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

10. 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). Matrizes definidas positivas

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

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

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

14. 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 :

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

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

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

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

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

20. Idéia geral do método • 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))

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

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

23. Encontrando w, parte 2 • Derivando f(xk+1) com relação a w, obtemos um ponto de mínimo dado por • Deste modo,

24. Vetores ortogonais • Constatamos que rk+1 é ortogonal aos vetores di (colunas de Dk), pois • Assim, djTrk+1 = 0, para j = 1, ..., k.

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

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

27. Como obter vetores A-conjugados • Lembrando que a iteração do método é definida por 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 (1)

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

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

30. Algoritmo • Dados de entrada: • A matriz A (simétrica, definida positiva). • O vetor b. • 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).

31. 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 • 3.6 k  k + 1

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

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

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

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