Otimização/Método de gradientes conjugados

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.

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$

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)\leq 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 }Ax+a^{\top }x+\alpha$.

Resolução

Sendo uma função quadrática, tem-se $f(x)={\frac {1}{2}}x^{\top }Ax+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 }Bx=x^{\top }Ax$ (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 }Ax+a^{\top }x+\alpha$, segue em particular que $\nabla f=Ax+a$ e $\nabla ^{2}f=A=P^{\top }\Lambda P$, onde $\Lambda$ é 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 }Ax+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 Ax+a=0$, ou seja, $x$ é solução do sistema linear $Ax=-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 $Ax=-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})=Ax_{k+1}+a=A(x_{k}+t_{k}d_{k})+a=Ax_{k}+t_{k}Ad_{k}+a$

Logo, $t_{k}Ad_{k}=-(Ax_{k}+a)$ e multiplicando por $d_{k}^{\top }$ obtem-se $t_{k}d_{k}^{\top }Ad_{k}=-(d_{k}^{\top }Ax_{k}+d_{k}^{\top }a)$. Consequentemente, o valor de $t_{k}$ é dado por

$t_{k}={\frac {-(d_{k}^{\top }Ax_{k}+a^{\top }d_{k})}{d_{k}^{\top }Ad_{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 $\theta (t)=f(x_{k}+td_{k})$, ou seja, $\theta$ é 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 $\theta$ em relação a $t$, obtem-se

$\theta ^{\prime }(t)=\nabla f(x_{k}+td_{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}$.

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}Ad_{i}$

Donde $\nabla f(x_{k+1})=\nabla f(x_{1})+\sum _{i=1}^{k}t_{i}Ad_{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 }Ad_{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}(Bd_{i})^{\top }(Bd_{k})$

Fazendo $\delta =Bd$, 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}(Bd_{k})^{\top }(Bd_{k})=\nabla f(x_{1})^{\top }d_{k}+t_{k}d_{k}Ad_{k}$

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

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 $x_{0}\in \mathbb {R} ^{n}$
  Se $\nabla f(x_{0})=0$, então pare: ${\bar {x}}=x_{0}$
  Senão: $d_{0}=-\nabla f(x_{0})=-Ax_{0}-a$
  Calcular $t_{0}={\frac {\|\nabla f(x_{0})\|^{2}}{d_{0}^{\top }Ad_{0}}}$
  $x_{1}=x_{0}+t_{0}d_{0}$


Passo iterativo $k$: Dado $x_{k}\in \mathbb {R} ^{n}$
  Se $\nabla f(x_{k})=0$, então pare: ${\bar {x}}=x_{k}$
  Senão: $d_{k}=-\nabla f(x_{k})+{\frac {\nabla f(x_{k})^{\top }Ad_{k}}{d_{k}^{\top }Ad_{k}}}d_{k}$
  $t_{k}={\frac {\|\nabla ^{2}f(x_{k})\|^{2}}{d_{k}^{\top }Ad_{k}}}$
  $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 }Ad_{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 }Ad_{k}}{d_{k}^{\top }Ad_{k}}}Ad_{k}$. Logo, $d_{k}^{\top }Ad_{k+1}=-\nabla f(x_{k+1})^{\top }Ad_{k}+{\frac {\nabla f(x_{k+1})^{\top }Ad_{k}}{d_{k}^{\top }Ad_{k}}}d_{k}^{\top }Ad_{k}=-\nabla f(x_{k+1})^{\top }Ad_{k}+\nabla f(x_{k+1})^{\top }Ad_{k}=0$.

Exercício

Provar que se $y=Bx$ 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_{i}t_{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$.

---

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 }Au$, 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}minf(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})=Ax_{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 }Ad_{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})=Ax_{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 }Ad_{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 }Au$. Observe que tal matriz é simétrica e definida positiva:

$\det(A-\lambda I)=(2-\lambda )(3-\lambda )-1=\lambda ^{2}-4\lambda -3=(\lambda -3)(\lambda -1)$

Logo, os autovalores de $A$ são $\lambda =1$ e $\lambda =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 $\beta$, 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 }Ad_{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 }Ad_{k}}}={\frac {\|\nabla f(x_{k})\|^{2}}{d_{k}^{\top }Ad_{k}}}$
• ${\frac {\nabla f(x_{k})^{\top }Ad_{k-1}}{d_{k-1}^{\top }Ad_{k-1}}}={\frac {\|\nabla f(x_{k})\|^{2}}{\|\nabla f(x_{k-1})\|^{2}}}$
• ${\frac {\nabla f(x_{k})^{\top }Ad_{k}}{d_{k-1}^{\top }Ad_{k-1}}}={\frac {\nabla f(x_{k})^{\top }(\nabla f(x_{k})-\nabla f(x_{k-1}))}{\|\nabla f(x_{k-1})\|^{2}}}$

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

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 $x_{0}\in \mathbb {R} ^{n}$
 Se $\nabla f(x_{0})=0$, então pare: ${\bar {x}}=x_{0}$
 Senão: $d_{0}=-\nabla f(x_{0})$ (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 $\nabla f(x_{k})=0$, então pare: ${\bar {x}}=x_{k}$
 Senão: $d_{k}=-\nabla f(x_{k})+{\frac {\|\nabla f(x_{k})\|^{2}}{\|\nabla f(x_{k-1})\|^{2}}}d_{k-1}$
 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.

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 $x_{0}\in \mathbb {R} ^{n}$
 Se $\nabla f(x_{0})=0$, então pare: ${\bar {x}}=x_{0}$
 Senão: $d_{0}=-\nabla f(x_{0})$ (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 $\nabla f(x_{k})=0$, então pare: ${\bar {x}}=x_{k}$
 Senão: $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}$
 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.

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:

• $pred=-t\theta$
• $\theta (t)=f(x+td)$
• $\theta '(t)=\nabla f(x+td)^{\top }d$

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