Otimização/Método de gradientes conjugados

Origem: Wikilivros, livros abertos por um mundo aberto.

< Otimização

Ir para: navegação, pesquisa

Otimização

A Wikipédia tem mais sobre este assunto:

Método do gradiente conjugado

Índice

1 Algumas considerações históricas
2 O método
- 2.1 Exemplos de conjuntos convexos e côncavos
3 Esquema do método de descida
4 Algoritmo de Hestenes-Stiefel
- 4.1 Exemplos
- 4.2 Implementação em Scilab
5 Algoritmo de Fletcher-Reeves
6 Algoritmo de Polak-Ribière
7 Algoritmo auxiliar

[editar] Algumas considerações históricas

Este método foi originalmente proposto por Hestenes e Stiefel, em 1952.
Seu objetivo inicial foi a resolução de problemas quadráticos sem restrições, mas logo o mesmo foi estendido para casos mais gerais.

[editar] O método

Este método pode ser considerado sob dois pontos de vista:

Como um método de descida, com busca linear exata;
Como um método de resolução de sistema linear, baseado em um processo de ortogonalização.

Definição

Um conjunto não vazio $D \subset \mathbb{R}^n$ é dito convexo quando $\forall x,y \in D$ e $t \in \left[0,1\right]$ vale

$tx+(1-t)y \in D$

[editar] Exemplos de conjuntos convexos e côncavos

Este é um conjunto convexo, pois todo segmento com extremidades no conjunto está totalmente contido no conjunto.

Este é um conjunto côncavo, pois existe um segmento com extremidades no conjunto que não está totalmente contido no conjunto.

Definição

Uma função $f: D \mapsto \mathbb{R}$ é dita convexa quando $D \subset \mathbb{R}^n$ é convexo e $\forall x,y \in D$ e $\forall t \in \left[0,1\right]$ vale

$f(tx+(1-t)y) \le tf(x) + (1-t) f(y)$

Definição

Dado um conjunto convexo $D \subset \mathbb{R}^n$ , uma função $f: D \mapsto \mathbb{R}$ é dita fortemente convexa quando existe uma constante $a > 0$ tal que $f(x)-a\|x\| ^2$ é convexa.

Exercício

Verifique que uma função quadrática $f: \mathbb{R}^n \mapsto \mathbb{R}$ é fortemente convexa se existe uma matriz simétrica definida positiva $A$ , um vetor $a \in \mathbb{R}^n$ e um escalar $\alpha \in \mathbb{R}\$ de modo que $f(x)=\frac{1}{2}x^\top A x + a^\top x + \alpha$ .

Resolução

Sendo uma função quadrática, tem-se $f(x)=\frac{1}{2}x^\top A x + a^\top x + \alpha$ . A matriz

A

pode ser suposta simétrica, pois caso não seja, toma-se $B = \frac{A+A^\top}{2}$ (simétrica), e segue $x^\top B x = x^\top A x$ (verifique).

Além disso, se $f$ é uma função fortemente convexa, então é estritamente convexa. Como $f$ é duas vezes diferenciável (por ser uma função quadrática), a convexidade estrita implica que $\nabla^2 f = A$ é definida positiva.

Nota: Uma matriz é definida positiva se, e somente se, todos os seus autovalores são positivos.

Tem-se:

$\nabla f : \mathbb{R}^n \mapsto \mathbb{R}$

$\nabla^2 f : \mathbb{R}^n \mapsto (\mathbb{R}^n)^2$

Sendo $f(x)=\frac{1}{2}x^\top A x + a^\top x + \alpha$ , segue em particular que $\nabla f = A x + a$ e $\nabla^2 f = A = P^\top\Lambda P$ , onde $Λ$ é uma matriz diagonal cujos elementos da diagonal são os autovalores de $A$ e $P$ é uma matriz onde as colunas são os autovetores correspondentes aos autovalores.

Note que $A$ é uma matriz simétrica, pois é a matriz Hessiana de uma função com segundas derivadas parciais contínuas, e consequentemente vale $\frac{\partial^2 f}{\partial x \partial y} = \frac{\partial^2 f}{\partial y \partial x}$ .

Para introduzir o método de direções conjugadas, serão consideradas somente funções quadráticas.

Uma condição necessária de primeira ordem para que $x$ seja um ponto de mínimo para a função $f$ é que $\nabla f(x) = 0$ . Para o presente caso, a função $f$ é convexa, então, a condição necessária $\nabla f(x) = 0$ também é suficiente.

Exercício

Prove que se $A$ é uma matriz simétrica definida positiva, então $f: \mathbb{R}^n \mapsto \mathbb{R}$ dada por $f(x)=x^\top A x + a^\top x + \alpha$ possui um único ponto de mínimo.

Resolução

Uma vez que

A

é simétrica definida positiva, a função

f

é fortemente convexa. Mas toda função fortemente convexa, definida em um conjunto fechado não vazio possui um único minimizador, pois:

Os conjuntos de nível de uma função fortemente convexa são compactos;
Toda função contínua definida em um compacto tem algum minimizador (pelo teorema de Weierstrass);
Os minimizadores de uma função convexa são globais;
Funções fortemente convexas não podem ter mais de um minimizador.

Em particular, $f$ possui um único ponto de mínimo.

No caso de uma função quadrática, tem-se $\nabla f(x) = 0 \Leftrightarrow A x + a = 0$ , ou seja, $x$ é solução do sistema linear $A x = - a$ .

A resolução de um sistema linear nem sempre pode ser feita numericamente de forma eficiente. Por exemplo, se a matriz do sistema é:

$A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 10^{20}+1 \end{bmatrix}$

A solução do sistema linear corresponde à interseção entre duas retas quase paralelas, e os erros de truncamento podem causar imprecisão na solução obtida computacionalmente.

Analiticamente, o sistema $A x = - a$ tem $x = - A - 1 a$ como solução. Então alguém poderia se perguntar: qual o problema em resolver esse sistema linear, se basta calcular a inversa da matriz $A$ e multiplicar pelo vetor $- a$ ? A resposta é que o calculo da inversa de uma matriz em geral é impraticável computacionalmente, por ter custo muito alto. Por isso, nas situações práticas, onde as matrizes tem ordem bem maior do que 2 (digamos 1000), o cálculo de matrizes inversas não é uma opção.

Assim, com o intuito de desenvolver um método computacional para o cálculo de minimizadores, é preciso utilizar outras técnicas. Considere o seguinte:

Em um método de descida tem-se sempre uma sequencia $\{ x_k, t_k, d_k \}\in \mathbb{N}$ , com $x k + 1 = x k + t k d k$ e $t k$ é um minimizador de ${f(x_k + td_k): t\in \mathbb{R}}$

$\nabla f(x) = Ax + a$

e

$0 = \nabla f(x_{k+1}) = A x_{k+1} + a = A(x_k + t_k d_k) + a = Ax_k + t_k A d_k + a$

Logo, $t k A d k = - (A x k + a)$ e multiplicando por $d_k^\top$ obtem-se $t_k d_k^\top A d_k = -(d_k^\top A x_k + d_k^\top a)$ . Consequentemente, o valor de $t k$ é dado por

$t_k= \frac{-(d_k^\top A x_k + a^\top d_k)}{ d_k^\top A d_k }$

Deste modo, o método consistirá de escolher em cada etapa $k$ uma direção $d k$ , e calcular o coeficiente $t k$ pela fórmula anterior, para gerar o próximo ponto $x k + 1$ . Mas como escolher a direção $d k$ ?

Dado $x k$ e escolhido $d k$ , defina $\theta: R \mapsto \mathbb{R}$ como $θ(t) = f (x k + t d k)$ , ou seja, $θ$ é a restrição da função $f$ à reta que passa pelo ponto $x k$ e que tem direção $d k$ . Logo, derivando a expressão de $θ$ em relação a $t$ , obtem-se

$\theta^\prime(t) = \nabla f(x_k + t d_k)^\top d_k$

Então, no ponto de mínimo, $x k + 1$ , tem-se

$0 = \nabla f(x_{k+1})^\top d_k$

Ou seja, a direção $d k$ a ser seguida a partir do ponto $x k$ é ortogonal ao gradiente da função $f$ , no ponto $x k + 1$ .

[editar] Esquema do método de descida

$x_{k+1} = x_k + t_k d_k = (x_{k-1} + t_{k-1} d_{k-1}) + t_k d_k = \ldots = x_1 + t_1 d_1 + \ldots + t_k d_k = x_1 + \sum_{i = 1}^{k}t_i d_i$

Seja $\bar{x}$ o minimizador da função $f$ . Tem-se

$x_{k+1} - \bar{x} = x_1 - \bar{x} + \sum_{i = 1}^{k}t_i d_i$

Mas $0 = \nabla f(\bar{x}) = A\bar{x} + a$ implica que $a = -A\bar{x}$ , logo

$\nabla f(x) = Ax + a = Ax - A\bar{x} = A(x - \bar{x})$

e consequentemente

$A(x_{k+1} - \bar{x}) = A(x_1 - \bar{x}) + \sum_{i = 1}^{k}t_i A d_i$

Donde $\nabla f(x_{k+1}) = \nabla f(x_{1}) + \sum_{i = 1}^{k}t_i A d_i$ . Portanto $0 = \nabla f(x_{k+1})^\top d_k = \nabla f(x_{1})^\top d_k + \sum_{i = 1}^{k}t_i d_i^\top A d_k$ .

Exercício

Provar que se $A$ é uma matriz simétrica, definida positiva, então existe uma matriz simétrica $B$ , de modo que $A = B 2$

Resolução

Sendo uma matriz simétrica, tem-se $A = P^\top \Lambda P$ , com

P

unitária e

$\Lambda =\begin{bmatrix} \lambda_1 & & 0 \\ & \ddots & \\ 0 & & \lambda_n \end{bmatrix} = \begin{bmatrix} \sqrt{\lambda_1} & & 0 \\ & \ddots & \\ 0 & & \sqrt{\lambda_n} \end{bmatrix}^2 = \Lambda_1^2$

Logo $A = P^\top\Lambda_1\Lambda_1 P = (P^\top\Lambda_1 P) (P^\top \Lambda_1 P) = B^2$

Usando o resultado desse exercício, tem-se ainda que $0 = \nabla f(x_{k+1})^\top d_k = \nabla f(x_{1})^\top d_k + \sum_{i = 1}^{k}t_i (B d_i)^\top (B d_k)$

Fazendo $δ = B d$ , o método do gradiente conjugado escolhe as direções de descida tais que $\delta_i^\top d_j = 0, \forall i \not = j$ . Mas quando $i \not = j$ , tem-se na expressão apresentada anteriormente apenas $0 = \nabla f(x_1)^\top d_k + t_k (B d_k)^\top (B d_k) = \nabla f(x_1)^\top d_k + t_k d_k A d_k$

Finalmente, tem-se o algoritmo para este método.

[editar] Algoritmo de Hestenes-Stiefel

Uma comparação da convergência do método de descida do gradiente com tamanho de passo ótimo (em verde) e o método do gradiente conjugado (em vermelho) para a minimização da forma quadrática com um sistema linear dado. O gradiente conjugado, assumindo aritmética exata, converge em no máximo n passos onde n é o tamanho da matriz do sistema (no exemplo, n=2).

Primeiro passo: Escolha 
  Se , então pare: 
  Senão: 
  Calcular 
   $x 1 = x 0 + t 0 d 0$ 


Passo iterativo  $k$ : Dado 
  Se , então pare: 
  Senão: 
  
   $x k + 1 = x k + t k d k$

Pode-se verificar facilmente que $d_{k+1} \perp d_k$ . De fato, como $d_{k+1} = -\nabla f(x_{k+1}) + \frac{\nabla f(x_{k+1})^\top A d_k}{d_k^\top Ad_k}d_k$ , tem-se $Ad_{k+1} = -A\nabla f(x_{k+1}) + \frac{\nabla f(x_{k+1})^\top A d_k}{d_k^\top Ad_k}Ad_k$ . Logo, $d_k^\top Ad_{k+1} = -\nabla f(x_{k+1})^\top A d_k + \frac{\nabla f(x_{k+1})^\top A d_k}{d_k^\top Ad_k} d_k^\top Ad_k = -\nabla f(x_{k+1})^\top A d_k + \nabla f(x_{k+1})^\top A d_k = 0$ .

Exercício

Provar que se $y = B x$ então $\|y_1 - \bar{y}\|^2 = \|y_{k+1} - \bar{y}\|^2 + \sum_{i=1}^{k}t_i^2 \|\delta_i\|^2$ .

Resolução

Tem-se

$x_{k+1} = x_1 + \sum_{i = 1}^{k}t_i d_i$

Multiplicando ambos os membros por $B$ , e trocando $x 1$ de lugar com $x k + 1$ resulta:

$-Bx_1 = -Bx_{k+1} + \sum_{i = 1}^{k}t_i Bd_i$ ,

ou seja,

$-y_1 = -y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i$ ,

somando $\bar{y}$ em ambos os lados, segue que

$\bar{y}-y_1 = \bar{y}-y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i$ ,

Então

$\begin{array}{lcl} \|\bar{y}-y_1 \|^2 & = & (\bar{y}-y_1)^\top(\bar{y}-y_1) \\ \ & = & \left(\bar{y}-y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i\right)^\top\left(\bar{y}-y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i\right)\\ \ & = & \|\bar{y}-y_{k+1}\|^2 + \sum_{i = 1}^{k}t_i^2 \|\delta_i\|^2 + 2\sum_{i \not= j}^{k}t_it_j \delta_i^\top\delta_j \\ \ & = & \|\bar{y}-y_{k+1}\|^2 + \sum_{i = 1}^{k}t_i^2 \|\delta_i\|^2 \end{array}$

Sendo a última igualdade devida ao fato de $\delta_i^\top\delta_j=0$ para $i \not = j$ .

[editar] Exemplos

Considere $f: \mathbb{R}^2 \mapsto \mathbb{R}$ definida por $f(x,y) = \frac{1}{2}(x^2 + y^2) = \frac{1}{2} \begin{bmatrix} x & y\end{bmatrix}\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}\begin{bmatrix} x\\ y\end{bmatrix}$ . Em outros termos, tomando $u = \begin{bmatrix} x \\ y\end{bmatrix}$ , tem-se $f(u)=\frac{1}{2}u^\top A u$ , onde $A = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} = I_{2\times 2}$ .

Pode-se aplicar o método de direções conjugadas ao seguinte problema

$(P) \left\{\begin{matrix} min f(u)\\ u \in \mathbb{R}^2 \end{matrix}\right.$

Note, desde já, que o conjunto solução é $S = \{ \begin{bmatrix} 0 \\ 0\end{bmatrix}\}$ .

Inicio

Toma-se $x 0$ arbitrário, por exemplo, $x_0 = \begin{bmatrix} 2 \\ 1\end{bmatrix}$ .
Avalia-se o gradiente da função $f$ neste ponto inicial: $\nabla f (x_0) = A x_0 = I_{2\times 2} x_0 = x_0$

Iteração 1

$d_0 = -\nabla f(x_0) = \begin{bmatrix} -2 \\ -1\end{bmatrix}$
$t_0 = \frac{\| \nabla f(x_0) \|^2}{d_0^\top A d_0} = \frac{5}{5} = 1$
$x_1 = x_0 + t_0 d_0 = \begin{bmatrix} 2 \\ 1\end{bmatrix} + 1 \begin{bmatrix} -2 \\ -1\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$

A seguir, verifica-se se o gradiente se anula no novo ponto $x 1$ :

$\nabla f(x_1) = A \begin{bmatrix} 0 \\ 0\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$

Como o gradiente já é nulo, não é preciso fazer a segunda iteração, e o ponto $x 1$ é o (único) minimizador global de $f$ .

Em um caso mais geral, considerando $f: \mathbb{R}^2 \mapsto \mathbb{R}$ definida por $f(x,y) = \frac{\lambda}{2}(x^2 + y^2) = \frac{1}{2} \begin{bmatrix} x & y\end{bmatrix}\begin{bmatrix} \lambda & 0 \\ 0 & \lambda\end{bmatrix}\begin{bmatrix} x\\ y\end{bmatrix}$ , tem-se cálculos muito parecidos em cada passo.

O conjunto solução continua sendo $S = \{ \begin{bmatrix} 0 \\ 0\end{bmatrix}\}$ .

Inicio

Considere $x 0$ como no primeiro exemplo, ou seja, $x_0 = \begin{bmatrix} 2 \\ 1\end{bmatrix}$ .
Avalia-se o gradiente da função $f$ neste ponto inicial: $\nabla f (x_0) = A x_0 = \lambda x_0$

Iteração 1

$d_0 = -\nabla f(x_0) = \lambda \begin{bmatrix} -2 \\ -1\end{bmatrix}$
$t_0 = \frac{\| \nabla f(x_0) \|^2}{d_0^\top A d_0} = \frac{5 \lambda^2}{5 \lambda^3} = \frac{1}{\lambda}$
$x_1 = x_0 + \frac{1}{\lambda} \lambda d_0 = \begin{bmatrix} 2\lambda \\ \lambda\end{bmatrix} + 1 \begin{bmatrix} -2\lambda \\ -\lambda\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$

A seguir, verifica-se se o gradiente se anula no novo ponto $x 1$ :

$\nabla f(x_1) = A \begin{bmatrix} 0 \\ 0\end{bmatrix} = \lambda\begin{bmatrix} 0 \\ 0\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$

Novamente, o gradiente se anula já na primeira iteração, de modo que $x 1$ é o minimizador global de $f$ .

Exercício

Seja $A \in \mathbb{R}^{n\times n}$ uma matriz simétrica definida positiva, cujos autovalores são todos iguais. Então começando de qualquer ponto $x_0 \not = 0$ , o método fornece $x n - 1$ como solução.

Um terceiro exemplo pode ser dado tomando $A = \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix}$ e $f: \mathbb{R}^2 \mapsto \mathbb{R}$ definida por $f(u)=\frac{1}{2}u^\top A u$ . Observe que tal matriz é simétrica e definida positiva:

det(A - λ I) = (2 - λ)(3 - λ) - 1 = λ 2 - 4λ - 3 = (λ - 3)(λ - 1)

Logo, os autovalores de $A$ são $λ = 1$ e $λ = 3$ . Isso também implica que a função é fortemente convexa.

Aplicando o método:

Início

Toma-se um ponto arbitrário no plano, por exemplo $x_0 = \begin{bmatrix} 10 \\ 20\end{bmatrix}$ ;
Verifica-se se tal ponto é o minimizador global, avaliando nele o gradiente da função:

$\nabla f(x_0) = \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 10 \\ 20\end{bmatrix} = \begin{bmatrix} 0 \\ 30\end{bmatrix} \not = \begin{bmatrix} 0 \\ 0\end{bmatrix}$ .

Já que o gradiente não se anulou no chute inicial, é preciso escolher uma direção e um comprimento de passo para determinar a próxima aproximação:

Iteração 1: $d_0 = - \nabla f(x_0) = \begin{bmatrix} 0 \\ -30\end{bmatrix}$; $t_0 = \frac{\| \begin{bmatrix} 0 & 30\end{bmatrix} \|^2}{\begin{bmatrix} 0 & -30\end{bmatrix} \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 0 \\ -30\end{bmatrix}} = \frac{900}{\begin{bmatrix} 0 & -30\end{bmatrix} \begin{bmatrix} -30 \\ -60\end{bmatrix}} = \frac{900}{1800} = \frac{1}{2}$

Feitos esses cálculos, o próximo ponto é dado por

$x_1 = x_0 + t_0 d_0 = \begin{bmatrix} 10 \\ 20\end{bmatrix} + \frac{1}{2}\begin{bmatrix} 0 \\ -30\end{bmatrix} = \begin{bmatrix} 10 \\ 20\end{bmatrix} + \begin{bmatrix} 0 \\ -15\end{bmatrix} = \begin{bmatrix} 10 \\ 5\end{bmatrix}$

Para saber se será necessária uma nova iteração, ou se o minimizador foi encontrado, calcula-se o gradiente da função no ponto:

$\nabla f(x_1) = \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 10 \\ 5\end{bmatrix} = \begin{bmatrix} 15 \\ 0\end{bmatrix} \not = \begin{bmatrix} 0 \\ 0\end{bmatrix}$ .

Novamente, será preciso calcular uma nova direção e um novo comprimento de passo:

Iteração 2: $d_0 = \begin{bmatrix} -15 \\ 0\end{bmatrix} + \beta \begin{bmatrix} 0 \\ -30\end{bmatrix} = \begin{bmatrix} -15 \\ -30\beta\end{bmatrix}$

onde $β$ , no algoritmo de Hestenes é dado por:

$\beta = \frac{\begin{bmatrix} 15 & 0\end{bmatrix} \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 0 \\ -30\end{bmatrix}}{\begin{bmatrix} 0 & -30\end{bmatrix} \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 0 \\ -30\end{bmatrix}} = \frac{\begin{bmatrix} 15 & 0\end{bmatrix} \begin{bmatrix} 30 \\ -60\end{bmatrix}}{\begin{bmatrix} 0 & -30\end{bmatrix} \begin{bmatrix} 30 \\ -60\end{bmatrix}} = \frac{15\times 30}{(-30) \times (-60)} = \frac{1}{4}$

Portanto

$d_0 = -15 \begin{bmatrix} 1 \\ 1/2\end{bmatrix}$

Além disso, o tamanho do passo é dado por

$t_1 = \frac{\| \nabla f(x_1) \|^2}{d_0^\top A d_0} = \frac{15^2}{15^2 \begin{bmatrix} 1 & 1/2\end{bmatrix} \begin{bmatrix} 2 & -1 \\ -1 & 2\end{bmatrix} \begin{bmatrix} 1 \\ 1/2\end{bmatrix}} = \frac{1}{\begin{bmatrix} 1 & 1/2\end{bmatrix} \begin{bmatrix} 3/2 \\ 0\end{bmatrix}} = \frac{1}{3/2} = \frac{2}{3}$

Portanto

$x_2 = x_1 + t_1 d_1 = \begin{bmatrix} 10 \\ 5\end{bmatrix} - 15 \frac{2}{3}\begin{bmatrix} 1 \\ 1/2\end{bmatrix} = \begin{bmatrix} 10 \\ 5\end{bmatrix} - 10 \begin{bmatrix} 1 \\ 1/2\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$

Obviamente, este é o minimizador procurado (pois o método tem a propriedade de convergência quadrática, ou seja utiliza no máximo $n$ iterações para chegar a solução quando aplicado a funções quadráticas definidas em $\mathbb{R}^n$ )

Exercício

Implementar o algoritmo de Hestenes-Stiefel em alguma linguagem de programação, por exemplo em Scilab, ou Matlab.

Exercício

Seja $f$ um função quadrática fortemente convexa. Verifique as seguintes igualdades:

$\frac{-\nabla f(x_k)^\top d_k}{d_k^\top A d_k} = \frac{\|\nabla f(x_k) \|^2}{d_k^\top A d_k}$
$\frac{\nabla f(x_k)^\top A d_{k-1}}{d_{k-1}^\top A d_{k-1}} = \frac{\|\nabla f(x_k) \|^2}{\|\nabla f(x_{k-1}) \|^2}$
$\frac{\nabla f(x_k)^\top A d_k}{d_{k-1}^\top A d_{k-1}} = \frac{\nabla f(x_k)^\top (\nabla f(x_k) - \nabla f(x_{k-1}))}{\|\nabla f(x_{k-1})\|^2}$

[editar] Implementação em Scilab

Abaixo é apresentada uma implementação deste algoritmo na linguagem de programação utilizada pelo Scilab.

A = [2 1; 1 2];
 
function [x] = min_gradiente_conjugado(xk)
  //Entrada: xk em R^n, qualquer "chute inicial"
  //  Saída: x, o ponto em que f assume o valor mínimo
 
  k        = 0        //Indica quantos loops já foram feitos
  epsilon  = 0.01
  n        = size(xk,1)
  g        = df(xk)
  seq      = zeros(n,n+1)
  seq(:,1) = xk
  while (norm(g) > epsilon) & (k <= n)
    if (k == 0)
      d = -g
    else
      d = Hestenes(g,d,A)
    end
    t  = busca_linear(g,d,A)
    xk = xk + t*d
    k  = k+1
    seq(:,k+1) = xk
    g  = df(xk)
  end
  plot(seq(1,:),seq(2,:))
  x = xk  
endfunction
 
 
function [fu] = f(u)
  fu=(1/2)*(u'*A*u)
endfunction
 
 
function [grf] = df(u)
  grf = A*u
endfunction
 
 
function [d] = Hestenes(g,d,A)
  d=-g + ((g'*A*d)/(d'*A*d))*d
endfunction
 
 
function [t] = busca_linear(g,d,A)
  t=(g'*g)/(d'*A*d)
endfunction

Para facilitar a compreensão do método, pode ser útil exibir as curvas de nível da função. Uma forma de implementar uma função com esse propósito é a seguinte:

function plotar_curvas
  qtd=101
  tam=max(seq)
  x=linspace(-tam,tam,qtd)
  y=x
  z=zeros(qtd,qtd)
  for i=1:qtd
    for j=1:qtd
      z(i,j)=f([x(i);y(j)])
    end
  end
  contour2d(x,y,z,10)
  a=gca();
  a.x_location = "middle"; 
  a.y_location = "middle"; 
endfunction

[editar] Algoritmo de Fletcher-Reeves

Um dos autores deste material sugeriu a adição de uma imagem para ilustrar o método de Fletcher-Reeves.

Esta versão é na verdade uma extensão do algoritmo anterior, permitindo a aplicação no caso de funções que não são quadráticas.

Primeiro passo: Escolha 
 Se , então pare: 
 Senão:  (como em todo método de descida)
 Calcular  $t 0$ , através de uma busca linear
  $x 1 = x 0 + t 0 d 0$ 
Passo iterativo:
 Se , então pare: 
 Senão: 
 Calcular  $t k$ , através de uma busca linear
  $x k + 1 = x k + t k d k$ 
  $k = k + 1$

Um dos autores deste material sugeriu a adição de uma implementação do algoritmo acima em SciLab.

[editar] Algoritmo de Polak-Ribière

Um dos autores deste material sugeriu a adição de uma imagem para ilustrar o método de Fletcher-Reeves.

Uma outra versão é a seguinte:

Primeiro passo: Tomar 
 Se , então pare: 
 Senão:  (como em todo método de descida)
 Calcular  $t 0$ , através de uma busca linear
  $x 1 = x 0 + t 0 d 0$ 
  $k = 1$ 
Passo iterativo:
 Se , então pare: 
 Senão: 
 Calcular  $t k$ , através de uma busca linear
  $x k + 1 = x k + t k d k$ 
  $k = k + 1$

Um dos autores deste material sugeriu a adição de uma implementação do algoritmo acima em SciLab.

Exercício

Verificar que, no caso de uma função $f$ quadrática e fortemente convexa, os algoritmos de Hestenes-Stiefel, de Fletcher-Reeves e de Polak-Ribière são os mesmos.

Exercício

Seja $f (x) = e - x + e x$ . Implemente o método de gradientes conjugados, e utilize o algoritmo para determinar o ponto de mínimo da função $f$ . Note que o espaço é unidimensional, então o método de gradientes conjugados reduz-se ao método dos gradientes, com primeira direção $-\nabla f(x_0)$ . Observe ainda que $f$ é uma função coerciva fortemente convexa.

[editar] Algoritmo auxiliar

Para o caso de funções não quadráticas, é preciso usar algum método de busca linear para a implementação do método dos gradientes conjugados, seja a versão de Fletcher-Reeves ou a de Polak-Ribière. Uma possibilidade é a busca de linear de Armijo (ver Izmailov & Solodov (2007), vol 2, pag. 65), cujo algoritmo é esboçado a seguir:

function busca_linear_Armijo (f, theta, alpha, delta, t0)
  while (alpha * pred > ared)
    t = d * t
  end
endfunction

com:

$p r e d = - t θ$
$θ(t) = f (x + t d)$
$\theta'(t) = \nabla f(x+td)^\top d$

Implementar a regra de Armijo e corrigir o esboço acima.

Otimização/Método de gradientes conjugados

Origem: Wikilivros, livros abertos por um mundo aberto.

Índice

[editar] Algumas considerações históricas

[editar] O método

[editar] Exemplos de conjuntos convexos e côncavos

[editar] Esquema do método de descida

[editar] Algoritmo de Hestenes-Stiefel

[editar] Exemplos

[editar] Implementação em Scilab

[editar] Algoritmo de Fletcher-Reeves

[editar] Algoritmo de Polak-Ribière

[editar] Algoritmo auxiliar

Acessos

Ferramentas pessoais

Pesquisa

Navegação

Projecto

Imprimir/exportar

Ferramentas