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) \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 $\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 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 $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_kd_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 $\theta (t) = f(x_k + t d_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 + 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}$.

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 $\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 (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.

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) = -A x_0 - a$
  Calcular $t_0 = \frac{\|\nabla f(x_0) \|^2}{d_0^\top A d_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 A d_k}{d_k^\top A d_k} d_k$
  $t_k = \frac{\|\nabla^2 f(x_k) \|^2}{d_k^\top A d_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 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 = 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_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$.

---

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

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.