Mostrar código-fonte de Otimização/Método de gradientes conjugados



{{wikipedia|Método do gradiente conjugado}}

==Algumas considerações históricas==

* Este método foi originalmente proposto por [[Otimização/Bibliografia#Hestenes & Stiefel (1952)|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.

==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 <math>D \subset \mathbb{R}^n</math> é dito convexo quando <math>\forall x,y \in D</math> e <math>t \in \left[0,1\right]</math> vale
:<math>tx+(1-t)y \in D</math>
}}

=== Exemplos de conjuntos convexos e côncavos ===

<center><gallery widths="200px" heights="200px" perrow="2">
Imagem:Conjunto convexo.svg|Este é um conjunto convexo, pois todo segmento com extremidades no conjunto está totalmente contido no conjunto.
Imagem:Conjunto côncavo.svg|Este é um conjunto côncavo, pois existe um segmento com extremidades no conjunto que não está totalmente contido no conjunto.
</gallery></center>


{{Definição
|Uma função <math>f: D \mapsto \mathbb{R}</math> é dita convexa quando <math>D \subset \mathbb{R}^n</math> é convexo e <math>\forall x,y \in D</math> e <math>\forall t \in \left[0,1\right]</math> vale
:<math>f(tx+(1-t)y) \le tf(x) + (1-t) f(y)</math>
}}

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

{{Exercício
|Verifique que uma função quadrática <math>f: \mathbb{R}^n \mapsto \mathbb{R}</math> é ''fortemente convexa'' se existe uma matriz simétrica definida positiva <math>A</math>, um vetor <math>a \in \mathbb{R}^n</math> e um escalar <math>\alpha \in \mathbb{R}\ </math> de modo que <math>f(x)=\frac{1}{2}x^\top A x + a^\top x + \alpha</math>.
}}
{{Resolução
|Sendo uma função quadrática, tem-se <math>f(x)=\frac{1}{2}x^\top A x + a^\top x + \alpha</math>. A matriz <math>A</math> pode ser suposta simétrica, pois caso não seja, toma-se <math>B = \frac{A+A^\top}{2}</math> (simétrica), e segue <math>x^\top B x = x^\top A x</math> (verifique).

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


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

Tem-se:
:<math>\nabla f  :  \mathbb{R}^n \mapsto \mathbb{R}</math>
:<math>\nabla^2 f  :  \mathbb{R}^n \mapsto (\mathbb{R}^n)^2</math>

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

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

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 <math>x</math> seja um ponto de mínimo para a função <math>f</math> é que <math>\nabla f(x) = 0</math>. Para o presente caso, a função <math>f</math> é convexa, então, a condição necessária <math>\nabla f(x) = 0</math> também é suficiente.

{{Exercício
|Prove que se <math>A</math> é uma matriz simétrica definida positiva, então <math>f: \mathbb{R}^n \mapsto \mathbb{R}</math> dada por <math>f(x)=x^\top A x + a^\top x + \alpha</math> possui um único ponto de mínimo.
}}
{{Resolução
|Uma vez que <math>A</math> é simétrica definida positiva, a função <math>f</math> é 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 [[w:Teorema de Weierstrass|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, <math>f</math> possui um único ponto de mínimo.
}}

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

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

:<math>A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 10^{20}+1 \end{bmatrix}</math>

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 <math>Ax = -a</math> tem <math>x = -A^{-1}a</math> 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 <math>A</math> e multiplicar pelo vetor <math>-a</math>?'' 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 <math>\{ x_k, t_k, d_k \}\in \mathbb{N}</math>, com <math>x_{k+1} = x_k + t_kd_k</math> e <math>t_k</math> é um minimizador de <math>{f(x_k + td_k): t\in \mathbb{R}}</math>

:<math>\nabla f(x) = Ax + a</math>

e

:<math>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</math>

Logo, <math>t_k A d_k = -(A x_k + a)</math> e multiplicando por <math>d_k^\top</math> obtem-se <math>t_k d_k^\top A d_k = -(d_k^\top A x_k + d_k^\top a)</math>. Consequentemente, o valor de <math>t_k</math> é dado por

:<math>t_k= \frac{-(d_k^\top A x_k + a^\top d_k)}{ d_k^\top A d_k }</math>

Deste modo, o método consistirá de escolher em cada etapa <math>k</math> uma direção <math>d_k</math>, e calcular o coeficiente <math>t_k</math> pela fórmula anterior, para gerar o próximo ponto <math>x_{k+1}</math>. Mas como escolher a direção <math>d_k</math>?

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

:<math>\theta^\prime(t) = \nabla f(x_k + t d_k)^\top d_k</math>

Então, no ponto de mínimo, <math>x_{k+1}</math>, tem-se

:<math>0 = \nabla f(x_{k+1})^\top d_k</math>

Ou seja, a direção <math>d_k</math> a ser seguida a partir do ponto <math>x_k</math> é ortogonal ao gradiente da função <math>f</math>, no ponto <math>x_{k+1}</math>.

==Esquema do método de descida==

:<math>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</math>

Seja <math>\bar{x}</math> o minimizador da função <math>f</math>. Tem-se

:<math>x_{k+1} - \bar{x} = x_1 - \bar{x} + \sum_{i = 1}^{k}t_i d_i</math>

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

:<math>\nabla f(x) = Ax + a = Ax - A\bar{x} = A(x - \bar{x})</math>

e consequentemente

:<math>A(x_{k+1} - \bar{x}) = A(x_1 - \bar{x}) + \sum_{i = 1}^{k}t_i A d_i</math>

Donde <math>\nabla f(x_{k+1}) = \nabla f(x_{1}) + \sum_{i = 1}^{k}t_i A d_i</math>. Portanto <math>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</math>.

{{Exercício
|Provar que se <math>A</math> é uma matriz simétrica, definida positiva, então existe uma matriz simétrica <math>B</math>, de modo que <math>A = B^2</math>
}}
{{Resolução
|Sendo uma matriz simétrica, tem-se <math>A = P^\top \Lambda P</math>, com <math>P</math> unitária e

:<math>\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</math>

Logo <math>A = P^\top\Lambda_1\Lambda_1 P = (P^\top\Lambda_1 P) (P^\top \Lambda_1 P) = B^2</math>
}}

Usando o resultado desse exercício, tem-se ainda que <math>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)</math>

Fazendo <math>\delta = Bd</math>, o método do gradiente conjugado escolhe as direções de descida tais que <math>\delta_i^\top d_j = 0, \forall i \not = j</math>. Mas quando <math>i \not = j</math>, tem-se na expressão apresentada anteriormente apenas <math>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</math>

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

==Algoritmo de Hestenes-Stiefel==
[[Image:Conjugate gradient illustration.svg|right|thumb|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 <math>x_0 \in \mathbb{R}^n</math>
   Se <math>\nabla f(x_0) = 0</math>, então pare: <math>\bar{x} = x_0</math>
   Senão: <math>d_0 = -\nabla f(x_0) = -A x_0 - a</math>
   Calcular <math>t_0 = \frac{\|\nabla f(x_0) \|^2}{d_0^\top A d_0}</math>
   <math>x_1 = x_0 + t_0 d_0</math>
 
 
 '''Passo iterativo <math>k</math>:''' Dado <math>x_k \in \mathbb{R}^n</math>
   Se <math>\nabla f(x_k) = 0</math>, então pare: <math>\bar{x} = x_k</math>
   Senão: <math>d_k = -\nabla f(x_k) + \frac{\nabla f(x_k)^\top A d_k}{d_k^\top A d_k} d_k</math>
   <math>t_k = \frac{\|\nabla^2 f(x_k) \|^2}{d_k^\top A d_k}</math>
   <math>x_{k+1} = x_k + t_k d_k</math>

Pode-se verificar facilmente que <math>d_{k+1} \perp d_k</math>. De fato, como <math>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</math>, tem-se <math>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</math>. Logo, <math>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</math>.

{{Exercício
|Provar que se <math>y = Bx</math> então <math>\|y_1 - \bar{y}\|^2 = \|y_{k+1} - \bar{y}\|^2 + \sum_{i=1}^{k}t_i^2 \|\delta_i\|^2</math>.
}}
{{Resolução
|Tem-se
:<math>x_{k+1} = x_1 + \sum_{i = 1}^{k}t_i d_i</math>
Multiplicando ambos os membros por <math>B</math>, e trocando <math>x_1</math> de lugar com <math>x_{k+1}</math> resulta:
:<math>-Bx_1 = -Bx_{k+1} + \sum_{i = 1}^{k}t_i Bd_i</math>,
ou seja,
:<math>-y_1 = -y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i</math>,
somando <math>\bar{y}</math> em ambos os lados, segue que
:<math>\bar{y}-y_1 = \bar{y}-y_{k+1} + \sum_{i = 1}^{k}t_i \delta_i</math>,
Então
:<math>\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}</math>

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

}}

----

===Exemplos===
Considere <math>f: \mathbb{R}^2 \mapsto \mathbb{R}</math> definida por <math>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}</math>.
Em outros termos, tomando <math>u = \begin{bmatrix} x \\ y\end{bmatrix}</math>, tem-se <math>f(u)=\frac{1}{2}u^\top A u</math>, onde <math>A = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} = I_{2\times 2}</math>.

Pode-se aplicar o método de direções conjugadas ao seguinte problema
:<math>(P) \left\{\begin{matrix}
min f(u)\\
u \in \mathbb{R}^2
\end{matrix}\right.</math>

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

;Inicio
* Toma-se <math>x_0</math> arbitrário, por exemplo, <math>x_0 = \begin{bmatrix} 2 \\ 1\end{bmatrix}</math>.
* Avalia-se o gradiente da função <math>f</math> neste ponto inicial: <math>\nabla f (x_0) = A x_0 = I_{2\times 2} x_0 = x_0</math>

;Iteração 1
* <math>d_0 = -\nabla f(x_0) = \begin{bmatrix} -2 \\ -1\end{bmatrix}</math>
* <math>t_0 = \frac{\| \nabla f(x_0) \|^2}{d_0^\top A d_0} = \frac{5}{5} = 1</math>
* <math>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}</math>

A seguir, verifica-se se o gradiente se anula no novo ponto <math>x_1</math>:
* <math>\nabla f(x_1) = A \begin{bmatrix} 0 \\ 0\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}</math>

Como o gradiente já é nulo, não é preciso fazer a segunda iteração, e o ponto <math>x_1</math> é o (único) minimizador global de <math>f</math>.

----

Em um caso mais geral, considerando <math>f: \mathbb{R}^2 \mapsto \mathbb{R}</math> definida por <math>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}</math>, tem-se cálculos muito parecidos em cada passo.

O conjunto solução continua sendo <math>S = \{ \begin{bmatrix} 0 \\ 0\end{bmatrix}\}</math>.

;Inicio
* Considere <math>x_0</math> como no primeiro exemplo, ou seja, <math>x_0 = \begin{bmatrix} 2 \\ 1\end{bmatrix}</math>.
* Avalia-se o gradiente da função <math>f</math> neste ponto inicial: <math>\nabla f (x_0) = A x_0 = \lambda x_0</math>

;Iteração 1
* <math>d_0 = -\nabla f(x_0) = \lambda \begin{bmatrix} -2 \\ -1\end{bmatrix}</math>
* <math>t_0 = \frac{\| \nabla f(x_0) \|^2}{d_0^\top A d_0} = \frac{5 \lambda^2}{5 \lambda^3} = \frac{1}{\lambda}</math>
* <math>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}</math>

A seguir, verifica-se se o gradiente se anula no novo ponto <math>x_1</math>:
* <math>\nabla f(x_1) = A \begin{bmatrix} 0 \\ 0\end{bmatrix} = \lambda\begin{bmatrix} 0 \\ 0\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix}</math>

Novamente, o gradiente se anula já na primeira iteração, de modo que <math>x_1</math> é o minimizador global de <math>f</math>.

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

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

Logo, os autovalores de <math>A</math> são <math>\lambda = 1</math> e <math>\lambda = 3</math>. 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 <math>x_0 = \begin{bmatrix} 10 \\ 20\end{bmatrix}</math>;
* Verifica-se se tal ponto é o minimizador global, avaliando nele o gradiente da função:
:<math>\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}</math>.
* 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
:<math>d_0 = - \nabla f(x_0) = \begin{bmatrix} 0 \\ -30\end{bmatrix}</math>
:<math>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}</math>
Feitos esses cálculos, o próximo ponto é dado por
:<math>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}</math>
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:
:<math>\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}</math>.
Novamente, será preciso calcular uma nova direção e um novo comprimento de passo:
;Iteração 2
:<math>d_0 = \begin{bmatrix} -15 \\ 0\end{bmatrix} + \beta \begin{bmatrix} 0 \\ -30\end{bmatrix} = \begin{bmatrix} -15 \\ -30\beta\end{bmatrix}</math>
onde <math>\beta</math>, no algoritmo de Hestenes é dado por:
:<math>\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}</math>
Portanto
:<math>d_0 = -15 \begin{bmatrix} 1 \\ 1/2\end{bmatrix}</math>
Além disso, o tamanho do passo é dado por
:<math>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}</math>

Portanto
:<math>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} </math>

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

{{Exercício
|Implementar o algoritmo de Hestenes-Stiefel em alguma linguagem de programação, por exemplo em [[w:Scilab|Scilab]], ou [[w:Matlab|Matlab]].
}}

{{Exercício
|Seja <math> f</math> um função quadrática fortemente convexa. Verifique as seguintes igualdades:
* <math> \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}</math>
* <math> \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}</math>
* <math> \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}</math>
}}

=== Implementação em Scilab ===
Abaixo é apresentada uma implementação deste algoritmo na linguagem de programação utilizada pelo [[w:Scilab|Scilab]].

<syntaxhighlight lang="matlab">
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
</syntaxhighlight>

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:

<syntaxhighlight lang="matlab">
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

</syntaxhighlight>

==Algoritmo de Fletcher-Reeves==

{{Tarefa||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 <math> x_0 \in \mathbb{R}^n</math>
  Se <math> \nabla f(x_0) = 0</math>, então pare: <math> \bar{x} = x_0</math>
  Senão: <math> d_0 = -\nabla f(x_0)</math> (como em todo método de descida)
  Calcular <math> t_0</math>, através de uma busca linear
  <math> x_1 = x_0 + t_0 d_0</math>
 '''Passo iterativo:'''
  Se <math> \nabla f(x_k) = 0</math>, então pare: <math> \bar{x} = x_k</math>
  Senão: <math> d_k = -\nabla f(x_k) + \frac{\|\nabla f(x_k)\|^2}{\|\nabla f(x_{k-1})\|^2}d_{k-1}</math>
  Calcular <math> t_k</math>, através de uma busca linear
  <math> x_{k+1} = x_k + t_k d_k</math>
  <math> k = k+1</math>

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

==Algoritmo de Polak-Ribière==

{{Tarefa||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 <math> x_0 \in \mathbb{R}^n</math>
  Se <math> \nabla f(x_0) = 0</math>, então pare: <math> \bar{x} = x_0</math>
  Senão: <math> d_0 = -\nabla f(x_0)</math> (como em todo método de descida)
  Calcular <math> t_0</math>, através de uma busca linear
  <math> x_1 = x_0 + t_0 d_0</math>
  <math> k = 1</math>
 '''Passo iterativo:'''
  Se <math> \nabla f(x_k) = 0</math>, então pare: <math> \bar{x} = x_k</math>
  Senão: <math> d_k = -\nabla f(x_k) + \frac{\nabla f(x_k)^\top (\nabla f(x_k) - \nabla f(x_{k-1}))}{\|\nabla f(x_{k-1})\|^2}d_{k-1}</math>
  Calcular <math> t_k</math>, através de uma busca linear
  <math> x_{k+1} = x_k + t_k d_k</math>
  <math> k = k+1</math>

{{Tarefa||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 <math> f</math> 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 <math>f(x) = e^{-x} + e^x</math>. Implemente o método de gradientes conjugados, e utilize o algoritmo para determinar o ponto de mínimo da função <math>f</math>. 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 <math>-\nabla f(x_0)</math>. Observe ainda que <math>f</math> é uma função coerciva fortemente convexa.
}}

== 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 [[#Algoritmo de Fletcher-Reeves|Fletcher-Reeves]] ou a de [[#Algoritmo de Polak-Ribière|Polak-Ribière]]. Uma possibilidade é a busca de linear de Armijo (ver [[Otimização/Bibliografia#Izmailov & Solodov (2007)|Izmailov & Solodov (2007), vol 2, pag. 65]]), cujo algoritmo é esboçado a seguir:

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

</syntaxhighlight>
com:
* <math>pred = -t\theta</math>
* <math>\theta(t) = f(x+td)</math>
* <math>\theta'(t) = \nabla f(x+td)^\top d</math>

{{Tarefa||Implementar a regra de Armijo e corrigir o esboço acima.}}

{{AutoCat}}