# 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$.
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.
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}$
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

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

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

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

## Algoritmo de Polak-Ribière

Uma outra versão é a seguinte:

### 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:
com:
* $pred=-t\theta$
* $\theta (t)=f(x+td)$
* $\theta '(t)=\nabla f(x+td)^{\top }d$