Otimização/Imprimir

Capa

Otimização

Introdução

Este material está sendo elaborado com base nas notas de aula da disciplina Otimização II, do curso de Matemática Industrial oferecido pela UFPR, ministrada pelo professor Wilfredo, no segundo semestre letivo do ano de 2008.

O conteúdo do livro não precisa (nem deve) se limitar àquele que consta atualmente no índice. Sendo assim, a qualquer momento o livro pode ser revisto e ampliado.

Sinta-se a vontade para ler este ou quaisquer outros livros do projeto, melhorando-os conforme lhe for possível. Com isso estará ajudando a aumentar a quantidade e a qualidade dos textos didáticos disponíveis em língua portuguesa, ao mesmo tempo em que colaborará com o crescimento projeto Wikilivros como um todo.

Se tiver dúvidas ou sugestões sobre páginas específicas, utilize as páginas de discussão correspondentes para deixar um comentário a respeito.

Ainda há muito por fazer, mas cada um daqueles que contribuem acredita estar fazendo o possível para oferecer o melhor a todos.

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.

Métodos de penalidades

Os métodos que recebem este nome fazem parte de uma família de métodos baseados em:

• Simplicidade conceitual;
• Eficiência prática;

Algumas das primeiras funções de penalidade foram desenvolvidas por:

• Courant (1943)
• Ablate Brighham (1955) (qual o artigo?)
• AV Fiacco, GP McCormick (1968) (qual o artigo?)

Para compreender este tipo de método, será considerado o seguinte problema:

$(P)\left\{{\begin{matrix}minf(x)\\g_{i}(x)\leq 0\forall i=1,\ldots ,p\\h_{j}(x)=0\forall j=1,\ldots ,q\end{matrix}}\right.$

Adicionalmente, se for considerado o conjunto de pontos $C=\{x\in \mathbb {R} ^{n}:g_{i}(x)\leq 0,\forall i=1,\ldots ,p;h_{j}(x)=0,\forall j=1,\ldots ,q\}$, o problema (P) se escreve ainda como

$(P)\left\{{\begin{matrix}minf(x)\\x\in C\end{matrix}}\right.$

Uma primeira idéia para a resolução desse tipo de problema é fazer uso de funções indicadoras, conforme é definido a seguir:

Definição

Dado um conjunto não vazio $C\subset \mathbb {R} ^{n}$, define-se a função indicadora $I_{C}:\mathbb {R} ^{n}\mapsto \mathbb {R} \cup \{\infty \}$ por:

$I_{C}(x)={\begin{cases}0,&{\mbox{se }}x\in C\\\infty ,&{\mbox{se }}x\not \in C\end{cases}}$

Notas: A função indicadora $I_{C}$ é às vezes denotada por $\chi _{C}$.

Para transformar um conjunto (P) em um problema sem restrições, pode-se proceder da seguinte maneira: Considerar o problema (PD), dado por:

$(PD)\left\{{\begin{matrix}minf(x)+I_{C}(x)\\x\in \mathbb {R} ^{n}\end{matrix}}\right.$

Considerando $i_{C}:\mathbb {R} \mapsto \mathbb {R} \cup \{\infty \}$ definida por:

$i_{C}(t)={\begin{cases}0,&{\mbox{se }}t\leq 0\\\infty ,&{\mbox{se }}t>0\end{cases}}$

Tem-se ainda o problema (PP), dado por:

$(PP)\left\{{\begin{matrix}minf(x)+\sum _{i=1}^{p}i(g_{i}(x))+\sum _{j=1}^{p}i(h_{j}(x))\\x\in \mathbb {R} ^{n}\end{matrix}}\right.$

Comentários:

• A grande vantagem desta idéia é que ela transforma um problema com restrições em um problema irrestrito.
• A principal desvantagem é que a função $\phi$ definida a seguir é descontínua em cada ponto $x\in \partial C$:

$\phi :\mathbb {R} ^{n}\mapsto \mathbb {R} \cup \{\infty \}$

$\phi (x)=f(x)+\sum _{i=1}^{p}i(g_{i}(x))+\sum _{j=1}^{p}i(h_{j}(x))$

Método de penalidade exterior

Para contornar a desvantagem da descontinuidade da função apresentada anteriormente, surgem outras funções, como por exemplo:

Gráfico de uma função de penalidade

$H:\mathbb {R} \mapsto \mathbb {R}$

$H(t)={\begin{cases}0,&{\mbox{se }}t\leq 0\\t^{2},&{\mbox{se }}t>0\end{cases}}$

Definição

Dada uma função $g:A\mapsto \mathbb {R}$, definida em um conjunto arbitrário $A$, define-se a parte positiva de $g$, como:

$g^{+}(x)=max\{0,g(x)\}$

Um dos autores deste material sugeriu a adição de uma imagem para comparar uma função $g(x)$ e sua correspondente $g^{+}(x)=max\{0,g(x)\}$.

Exercício

Verifique que para cada $x\in \mathbb {R} ^{n}$ tem-se, para cada $i$, a igualdade $H(g_{i}(x))=\left(g_{i}^{+}(x)\right)^{2}$, onde $g_{i}^{+}$ é a parte positiva de $g_{i}$ e $H$ é a função definida no exemplo anterior.

Resolução

De fato, dada qualquer função $g:A\mapsto \mathbb {R}$, segue da própria definição de $H$ que $H(g(x))={\begin{cases}0,&{\mbox{se }}g(x)\leq 0\\g(x)^{2},&{\mbox{se }}g(x)>0\end{cases}}$ Mas $g^{+}(x)=max\{0,g(x)\}={\begin{cases}0,&{\mbox{se }}g(x)\leq 0\\g(x),&{\mbox{se }}g(x)>0\end{cases}}$ implica que $(g^{+}(x))^{2}=max\{0,g(x)\}={\begin{cases}0,&{\mbox{se }}g(x)\leq 0\\(g(x))^{2},&{\mbox{se }}g(x)>0\end{cases}}$ Portanto, $H(g(x))=\left(g^{+}(x)\right)^{2}$. Como $g$ pode ser tomada arbitrariamente, tem-se em particular que $H(g_{i}(x))=\left(g_{i}^{+}(x)\right)^{2}$.

Exercício

Verifique que para cada $x\in \mathbb {R} ^{n}$ tem-se, para cada $j$, a igualdade $H\left(h_{j}(x)^{2}\right)=\left(h_{j}(x)\right)^{2}$.

Um dos autores deste material sugeriu conferir o enunciado do exercício anterior. A afirmação parece ser falsa nos casos em que $h_{j}(x)<0$.

A idéia de aplicar penalizações aos pontos que não pertencem ao conjunto viável é formalizada na seguinte definição:

Definição

Seja $H:\mathbb {R} ^{n}\mapsto \mathbb {R}$. A função $H$ é chamada de função de penalidade exterior se possui as seguintes propriedades:

• $H(x)\geq 0,\forall x\in \mathbb {R} ^{n}$
• $H(x)=0,\Leftrightarrow x\in C$
• $H$ é contínua

Nota: Lembre-se que $C=\{x\in \mathbb {R} ^{n}:g_{i}(x)\leq 0,\forall i=1,\ldots ,p;h_{j}(x)=0,\forall j=1,\ldots ,q\}$ é o conjunto viável do problema (P).

Em particular, as funções $H\circ g_{i}$ são funções de penalidade exterior.

Definição

Seja $\theta :\mathbb {R} ^{n}\mapsto \mathbb {R}$. A função $\theta$ é dita coerciva se $\lim _{\|x\|\to \infty }\theta (x)=+\infty$

Um dos autores deste material sugeriu a adição de exemplos de funções de penalidade, juntamente com algumas imagens ilustrando os seus gráficos.

Nota: Conforme o Wikcionário, o termo coercivo significa: que coage; que reprime; que impõe pena; coercitivo. Nesse sentido, esse é um termo adequado ao tratar do conceito anterior, no contexto dos métodos de penalidade.

Exercício

Verifique que $\theta$ é uma função coerciva se, e somente se, $L_{\theta }(\lambda )=\{x\in \mathbb {R} ^{n}:\theta (x)\leq \lambda \}$ é limitado para todo $\lambda \in \mathbb {R}$.

Resolução

Pela definição de limite, a afirmação $\lim _{\|x\|\to \infty }\theta (x)=+\infty$ é equivalente a dizer que $\forall K>0,\exists M>0$ tal que $\forall x$ tem-se $\|x\|>M\Rightarrow \theta (x)>K$. Esta última implicação, é equivalente à $\theta (x)\leq K\Rightarrow \|x\|\leq M$ (sua contrapositiva). Pela definição de conjunto de nível, isso equivale à $x\in L_{\theta }(K)\Rightarrow \|x\|\leq M$. A existência de um número $M$ com tal propriedade significa que $L_{\theta }(K)$ é um conjunto limitado, donde conclui-se a equivalência entre coercividade e limitação dos conjuntos de níveis

Um dos autores deste material sugeriu a adição de uma imagem que ilustre geometricamente a relação entre coercividade e conjuntos de nível.

Exercício

Verifique que $H:\mathbb {R} ^{n}\mapsto \mathbb {R}$ definida por $H(x)=\sum _{i=1}^{p}\left(g_{i}^{+}(x)\right)^{2}+\sum _{j=1}^{q}\left(h_{j}(x)\right)^{2}$ é uma função de penalidade exterior, onde conforme anteriormente, $g_{i}^{+}$ é a parte positiva de $g_{i}$.

Resolução

Primeiramente, $H(x)$ é uma função não negativa, pois o quadrado de qualquer número real é não negativo, assim como a soma de números reais não negativos. Em segundo lugar, tem-se $H(x)=0$ se, e somente se, para cada índice $i$ e cada índice $j$ vale $g_{i}^{+}(x)=0$ e $h_{j}(x)=0$. Como $g_{i}^{+}(x)=0\Leftrightarrow g_{i}\leq 0$, tem-se $H(x)=0\Leftrightarrow x\in C$ Finalmente, como a soma e o produto de funções contínuas resulta em uma função contínua, segue que $H(x)$ é contínua, pois $g_{i},h_{j}$ são contínuas.

Exercício

Verifique que se $g:\mathbb {R} ^{n}\mapsto \mathbb {R}$ é uma função contínua coerciva, então existe ${\bar {x}}\in \mathbb {R} ^{n}$ tal que $g({\bar {x}})\leq g(x),\forall x\in \mathbb {R} ^{n}$.

Resolução

Considere um ponto arbitrário ${\hat {x}}\in \mathbb {R} ^{n}$. Tome $\lambda =g({\hat {x}})$. Nestas condições, $L_{g}(\lambda )$ é limitado (conforme um dos exercícios anteriores) e fechado (pois é pré-imagem de um conjunto fechado por uma função contínua), portanto compacto. Neste caso, conforme o teorema de Weierstrass, a função $g$ possui algum ponto de mínimo no conjunto $L_{g}(\lambda )$, ou seja, existe ${\bar {x}}\in L_{g}(\lambda )$ tal que $g({\bar {x}})\leq g(x),\forall x\in L_{g}(\lambda )$. Ne verdade, tal ponto ${\bar {x}}$ é também um minimizador global da função $g$, pois se $x\not \in L_{g}(\lambda )$ então $g(x)>\lambda$ (pela definição do conjunto $L_{g}(\lambda )$).

Um dos autores deste material sugeriu a adição de uma imagem ilustrando a existência de minimizadores globais para funções coercivas.

Alguns conceitos da topologia

Em algumas situações, é interessante ter em mente que certos conceitos definidos no contexto da Otimização são, na verdade, instanciações de conceitos mais gerais, muitos deles provenientes da topologia. Alguns exemplos são apresentados a seguir.

Definição

Dado um conjunto $X$, uma coleção $\Gamma \subset {\mathcal {P}}$ de subconjuntos de $X$ é chamada de topologia se:

• $\emptyset ,X\in \Gamma$
• $\left\{A_{i}\right\}_{i\in I}\subset \Gamma \Rightarrow \cup _{i\in I}A_{i}\in \Gamma$
• $\left\{B_{i}\right\}_{i=1}^{p}\subset \Gamma \Rightarrow \cap _{i=1}^{p}B_{i}\in \Gamma$

Exemplos

• Topologia euclidiana:

$X=R$

$\Gamma =\Gamma _{1}=\{\emptyset \}\cup \{X\}\cup \left\{A\subset \mathbb {R} :\forall x\in A,\exists \epsilon _{x}>0,\left(x-\epsilon _{x},x+\epsilon _{x}\right)\subset A\right\}$

Em outras palavras, a topologia euclideana é a coleção de todos os conjuntos abertos contidos em $\mathbb {R}$. Pode-se verificar com facilidade que de fato são satizfeitas as três propriedades que definem uma topologia.

Outro exemplo muito comum é o seguinte:

• Topologia euclideana estendida:

$X=R\cup \{\infty \}$

$\Gamma =\Gamma _{1}\cup \left\{A\subset \mathbb {R} \cup \{\infty \}:\exists a\in \mathbb {R} ,A=\left]a,\infty \right]\right\}$

Em geral, a noção de limite seria caracterizada topologicamente da seguinte forma:

Definição

$\lim _{n\to \infty }x_{n}=a\Leftrightarrow \forall I\ni x,\exists n_{0}\in \mathbb {N} ,\forall n\geq n_{0}\Rightarrow x_{n}\in I$

De volta ao método

Uma vez apresentados os conceitos iniciais, pode-se provar o seguinte teorema:

Teorema

Considere:

• $H:\mathbb {R} ^{n}\mapsto \mathbb {R}$ uma função de penalidade exterior;
• $f:\mathbb {R} \mapsto \mathbb {R}$ uma função contínua;
• $C=\{x\in \mathbb {R} ^{n}:g_{i}(x)\leq 0,\forall i=1,\ldots ,p;h_{j}(x)=0,\forall j=1,\ldots ,q\}$ um conjunto fechado;
• $\left\{r_{k}\right\}$ uma sequência de termos positivos tal que $\lim _{k\to \infty }r_{k}=0$.

Suponha que é válida uma das seguintes propriedades:

1. $f$ é coerciva.
2. $C$ é limitado e $H$ é coerciva.

Se, para cada $k$, for escolhido ${\bar {x}}(r_{k})\in \arg \min\{f(x)+{\frac {1}{r_{k}}}H(x)\}$,então:

1. $\left\{{\bar {x}}(r_{k})\right\}$ possui algum ponto de acumulação;
2. Todo ponto de acumulação de $\left\{{\bar {x}}(r_{k})\right\}$ é solução do problema (P);
3. $\lim _{k\to \infty }H({\bar {x}}(r_{k}))=0$.

Um dos autores deste material sugeriu a adição de uma imagem que ilustre geometricamente o significado do teorema acima.

Demonstração

Se (1) acontece, então $f+{\frac {1}{r}}H$ é coerciva, pois $\lim _{\|x\|\to \infty }f(x)+{\frac {1}{r}}H(x)\geq \lim _{\|x\|\to \infty }f(x)=\infty$ A desigualdade é válida pois $H$ é uma função de penalidade, portanto não negativa, e a igualdade se deve à hipotese sobre $f$. Além de coerciva, tal função é contínua (pois é combinação linear de funções contínuas). Logo, a função possui um ponto de mínimo global, ou seja, existe ${\bar {x}}\in \mathbb {R} ^{n}$ tal que $f({\bar {x}})+{\frac {1}{r}}H({\bar {x}})\leq f(x)+{\frac {1}{r}}H(x)$ , para qualquer $x\in \mathbb {R} ^{n}$. Mas $f$ é contínua e coerciva, então também existe $x^{*}\in \mathbb {R} ^{n}$ tal que $f(x^{*})\leq f(x),\forall x\in C$ Por outro lado, se vale (2), então $f+{\frac {1}{r}}H$ é contínua e C compacto. Analogamente, $f$ é contínua e $C$ compacto, donde tem-se algum $x^{*}\in C$ tal que $f(x^{*})\leq f(x),\forall x\in C$ Em ambos os casos, dada uma sequência $\left\{r_{k}\right\}$ tal que $0<r_{k+1}<r_{k}$ e $\lim _{k\to \infty }r_{k}=\infty$, defina-se $\varphi :\mathbb {R} ^{n}\times \mathbb {R} \mapsto \mathbb {R}$ por $\varphi (x,r)=f(x)+{\frac {1}{r}}H(x)$. Logo, ${\begin{array}{lcl}\varphi ({\bar {x}}(r_{k}),r_{k})&=&f({\bar {x}}(r_{k}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))\\\ &\leq &f({\bar {x}}(r_{k+1}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k+1})\\\ &\leq &f({\bar {x}}(r_{k+1}))+{\frac {1}{r_{k+1}}}H({\bar {x}}(r_{k+1})\\\ &=&\varphi ({\bar {x}}(r_{k+1}),r_{k+1})\end{array}}$ Sendo que a primeira desigualdade se deve ao fato de ${\bar {x}}(r_{k})$ ser um minimizador, por construção, e a segunda segue por que $\{r_{k}\}$ é decrescente. Portanto, $\varphi ({\bar {x}}(r_{k}),r_{k})\leq \varphi ({\bar {x}}(r_{k+1}),r_{k+1}),\forall k$. Além disso, tem-se as seguintes desigualdades: $f({\bar {x}}(r_{k}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))\leq f({\bar {x}}(r_{k+1}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k+1}))$ $f({\bar {x}}(r_{k+1}))+{\frac {1}{r_{k+1}}}H({\bar {x}}(r_{k+1}))\leq f({\bar {x}}(r_{k}))+{\frac {1}{r_{k+1}}}H({\bar {x}}(r_{k}))$ Logo, somando os membros correspondentes, obtem-se: ${\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))+{\frac {1}{r_{k+1}}}H({\bar {x}}(r_{k+1}))\leq {\frac {1}{r_{k}}}H({\bar {x}}(r_{k+1}))+{\frac {1}{r_{k+1}}}H({\bar {x}}(r_{k}))$ ou seja, $\left({\frac {1}{r_{k+1}}}-{\frac {1}{r_{k}}}\right)H({\bar {x}}(r_{k+1}))\leq \left({\frac {1}{r_{k+1}}}-{\frac {1}{r_{k}}}\right)H({\bar {x}}(r_{k}))$ Portanto, $H({\bar {x}}(r_{k+1}))\leq H({\bar {x}}(r_{k})),\forall k$ Por outro lado, $f({\bar {x}}(r_{k}))\leq f({\bar {x}}(r_{k}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))\leq f(x^{*})+{\frac {1}{r_{k}}}H(x^{*})=f(x^{*}),\forall k$, sendo primeira desigualdade válida por $H$ ser não negativa e $r_{k}$ positivo, e a segunda devida à própria definição de ${\bar {x}}(r_{k})$. Logo, $f({\bar {x}}(r_{k}))\leq f(x^{*})$. Se a primeira das hipóteses acontece, segue da coercividade e dessa última desigualdade que $\{{\bar {x}}(r_{k})\}\subset L_{f}(f(x^{*}))$. Então $\{{\bar {x}}(r_{k})\}_{i=1}^{\infty }$ é uma sequência limitada. Se ocorre a segunda, então $H$ é coerciva, mas foi mostrado que $H({\bar {x}}(r_{k+1}))\leq H({\bar {x}}(r_{k}))$, consequentemente $\{{\bar {x}}(r_{k+1})\}\subset L_{H}(H({\bar {x}}(r_{1})))$. Sendo $H$ coerciva, conclui-se novamente que $\{{\bar {x}}(r_{k})\}_{i=1}^{\infty }$ é uma sequência limitada. Portanto, em qualquer caso, $\{{\bar {x}}(r_{k})\}_{i=1}^{\infty }$ possui algum ponto de acumulação. Seja ${\bar {x}}$ um ponto de acumulação de $\{{\bar {x}}(r_{k})\}_{i=1}^{\infty }$. Então existe $\{{\bar {x}}(r_{k_{n}})\}\subset \{{\bar {x}}(r_{k})\}$ tal que $\lim _{n\to \infty }{\bar {x}}(r_{k_{n}})={\bar {x}}$. Logicamente, $\lim _{n\to \infty }r_{k_{n}}=0$. Pela continuidade de $f$ e sabendo que $f({\bar {x}}(r_{k}))\leq f(x^{*})$, se deduz que $f({\bar {x}})=\lim _{n\to \infty }f({\bar {x}}(r_{k_{n}}))\leq f(x^{*})$. Mas já havia sido verificado que $f(x^{*})\leq f(x),\forall x\in C$, então segue a igualdade $f({\bar {x}})=f(x^{*})$. A sequência $\{\varphi ({\bar {x}}(r_{k}),r_{k})\}$ é crescente. Seja ${\bar {\varphi }}=\lim _{n\to \infty }\varphi ({\bar {x}}(r_{k}),r_{k})$ (por que razão ele existe?). Como $\varphi ({\bar {x}}(r_{k}),r_{k})=f({\bar {x}}(r_{k}))+{\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))$, tem-se ${\bar {\varphi }}=\lim _{n\to \infty }{\frac {1}{r_{k}}}H({\bar {x}}(r_{k}))=\lim _{n\to \infty }\varphi ({\bar {x}}(r_{k}),r_{k})-\lim _{n\to \infty }f({\bar {x}}(r_{k}))={\bar {\varphi }}-f({\bar {x}})$. Portanto, $\lim _{n\to \infty }H({\bar {x}}(r_{k_{n}}))=0$. Logo $\lim _{k\to \infty }H({\bar {x}}(r_{k}))$, ou seja, $\lim _{r\to 0}H({\bar {x}}(r))$. Como $H$ é contínua, $\lim _{n\to \infty }H({\bar {x}}(r_{k_{n}}))=H({\bar {x}})$, e este valor é nulo se, e somente se, ${\bar {x}}\in C$. Logo, $f({\bar {x}})\leq f(x^{*})\leq f({\bar {x}})$, donde $f({\bar {x}})=f(x^{*})$. Assim, tem-se $f({\bar {x}})=f(x^{*})\leq f(x),\forall x\in C$, ou seja, ${\bar {x}}$ é solução de (P).

Exercício

Dado o problema (P), considere $m=\inf \left\{f(x):g_{i}(x)\leq 0,\forall i\ {\text{e}}\ h_{j}=0,\forall j\right\}$ (isso não quer dizer que o problema tenha solução). suponha-se que $f,g$ e $h$ são funções contínuas e que $C=\left\{x:g_{i}(x)\leq 0,\forall i;h_{j}(x)=0,\forall j\right\}$ seja não vazio, ou seja, que $C$ é factível. Tome $H:\mathbb {R} ^{n}\mapsto \mathbb {R}$ como $H(x)=\sum _{i=1}^{p}[g_{i}^{+}(x)]^{2}+\sum _{j=1}^{q}[h_{j}(x)]^{2}$, onde $g_{i}^{+}$ denota a parte positiva de $g_{i}$, como de costume. Considere ainda $\varphi :\mathbb {R} ^{n}\times \mathbb {R} ^{+}\mapsto \mathbb {R}$ dada por $\varphi (x,r)=f(x)+{\frac {1}{r}}H(x)$ e, para cada $r>0$, seja $m(r)=\inf _{x\in \mathbb {R} ^{n}}\varphi (x,r)$ Nessas condições, provar que:

1. $H$ é uma função de penalidade exterior
2. $m>-\infty$
3. $\varphi (x,r)=f(x),\forall x\in C,\forall r>0$
4. Se $0<r<s$ então $m(s)<m(r)$
5. Se $f,g_{i},h_{j}$ são covexas, então $\varphi$ é convexa.
6. Se $f,g_{i},h_{j}$ são diferenciáveis, então $\varphi$ é diferenciável em $x$, e $\nabla _{x}\varphi (x,r)=\nabla f(x)+\sum _{i=1}^{p}g_{i}^{+}(x)\nabla g_{i}(x)+\sum _{j=1}^{q}h_{j}(x)\nabla h_{j}(x)$
7. Se $0<r_{k+1}<r_{k}$ e $\lim _{k\to \infty }r_{k}=0$ e se $\{x_{k}\}$ é uma sequência tal que $\varphi (x_{k},r_{k})=m(r_{k})$ e $\lim _{k\to \infty }x_{k}={\bar {x}}$ então ${\bar {x}}$ é solução de (P).

Algoritmo de penalidade exterior

Primeiro passo: Escolha $x_{0}\in \mathbb {R} ^{n}$, $r>0$ e $k=1$

Passo iterativo $k$: Calcular $x_{k}={\bar {x}}(r_{k})$, solução global de 
$(P)\left\{{\begin{matrix}\min \varphi (x,r_{k})&=&f(x)+{\frac {1}{2r}}H(x)\\x\in \mathbb {R} ^{n}&&\end{matrix}}\right.$
Escolha $r_{k}$ tal que $0<r_{k}<r_{k-1}$

Nota: Neste algoritmo, $H$ é uma função de penalidade exterior.

Exercício

Sejam $f(x,y)=x^{2}+y^{2}$ e $C=\{(x,y)\in \mathbb {R} ^{2}:-x\leq 0;-y\leq 0;x+y-1\leq 0\}$. Considere o problema: $(P)\left\{{\begin{matrix}\min f(x,y)\\(x,y)\in \mathbb {R} ^{n}\end{matrix}}\right.$ Utilize o método de penalidade exterior para determinar o ponto de mínimo da função $f$ sobre o conjunto $C$.

Resolução
Esboço: Pode-se tomar $H(x,y)=\sum _{i=1}^{3}(g_{i}^{+}(x,y))^{2}$, onde: $g_{1}(x,y)=-x$ $g_{2}(x,y)=-y$ $g_{3}(x,y)=x+y-1$ Tem-se então o problema irrestrito: $(P_{r})\left\{{\begin{matrix}\min f(x,y)+{\frac {1}{r}}H(x,y)\\(x,y)\in \mathbb {R} ^{2}\end{matrix}}\right.$, onde pode ser aplicado o método de gradientes conjugados.

Implementação do algoritmo de penalidade exterior em SciLab

Considerando o problema

$(P)\left\{{\begin{matrix}minf(x)\\g_{i}(x)\leq 0\forall i=1,\ldots ,p\\h_{j}(x)=0\forall j=1,\ldots ,q\end{matrix}}\right.$

com $f$ uma função quadrática, $A$ simétrica positiva definida, $g_{i},h_{j}$ lineares, e a função de penalidade, $H$, dada por:

$H(x)=\sum _{i=1}^{p}(g_{i}^{+}(x,y))^{2}+\sum _{j=1}^{q}(h_{j}(x,y))^{2}$

Pode-se usar o método dos gradientes conjugados para resolver o seguinte problema:

$(P)\left\{{\begin{matrix}minf(x)+H(x)\\x\in \mathbb {R} ^{n}\end{matrix}}\right.$

onde $f+H$ terá sua matriz Hessiana definida positiva.

Incluir o código para uma implementação do método em SciLab.

Método de penalidade interior

Este método também é conhecido como método de barreira. Ele consiste em trabalhar com funções de penalidade tais que $\forall x\in \partial C$ e qualquer que seja a sequência $\{x_{n}\}\subset C=\{x:g_{i}(x)<0\}$ para a qual $x_{n}\to x$, se tem que a função de penalidade tende a $+\infty$.

Um dos autores deste material sugeriu a adição de uma imagem para ilustrar uma sequência de pontos que tende a um ponto da fronteira de um conjunto C, de preferência junto com o gráfico de uma função de penalidade deste novo tipo.

Mas como conseguir esse tipo de função?

Exemplos

Considere o caso em que $g_{i}(x)=a_{i}^{\top }x-b_{i}$. Lembrando que a função logarítmica tem a propriedade:

$\lim _{x\to 0^{+}}-\ln(x)=+\infty$

pode-se tomar $H_{1}(x)=-\sum _{i=1}^{p}\ln(b_{i}-a_{i}^{\top }x)$, pois

$\lim _{a_{i}^{\top }x\to b_{i}}-\ln(b_{i}-a_{i}^{\top }x)=+\infty$

implica que

$\lim _{x\to \partial C}H_{1}(x)=+\infty$

Analogamente, tem-se

$\lim _{x\to 0^{+}}{\frac {1}{x}}=+\infty$

então, uma outra função com a propriedade desejada é $H_{2}(x)=\sum _{i=1}^{p}{\frac {1}{b_{i}-a_{i}^{\top }x}}$.

O problema a ser resolvido quando se quer aplicar o método de barreira é o seguinte:

$(P)\left\{{\begin{matrix}\min f(x)\\g_{i}(x)\leq 0;i=1,\ldots ,p\end{matrix}}\right.$

onde $C=\{x:g_{i}(x)\leq 0,\forall i\}$ é tal que $C^{0}=\{x:g_{i}(x)<0,\forall i\}\not =\emptyset$

Observação: $C^{0}$ não é necessariamente igual ao interior do conjunto $C$. Por exemplo:

Um dos autores deste material sugeriu a adição de uma imagem para ilustrar um conjunto C cujo correspondente C0 seja diferente do interior de C.

Definição

Se diz que uma função $B:C^{0}\mapsto \mathbb {R}$ é uma função de barreira se ela possui as seguintes propriedades:

1. $B$ é limitada inferiormente em $C$.
2. Para qualquer sequência $\{x_{n}\}\subset C^{0}$ tal que $\lim _{n\to \infty }x_{n}=x\in \partial C$ vale $\lim _{n\to \infty }B(x_{n})=+\infty$
3. $B$ é contínua

Algoritmo de barreira

Primeiro passo: Escolha $x_{0}\in C^{0}$, $r>0$ e $t_{0}>0$

Passo iterativo $k$: Calcular $x_{k+1}={\bar {x}}(t_{k})$, solução global de 
$(P)\left\{{\begin{matrix}\min f(x)+t_{k}B(x)\\x\in C^{0}\end{matrix}}\right.$
Escolha $t_{k+1}$ tal que $0<t_{k+1}<t_{k}$

Observação: Neste método os termos da sequência $\{x_{n}\}$ estão sempre em $C^{0}$.

Implementação do algoritmo de barreira em SciLab

Considerando o problema

$\left\{{\begin{matrix}minf(x)\\g_{i}(x)\leq 0\forall i=1,\ldots ,p\\\end{matrix}}\right.$

pode-se usar o método de barreira no conjunto:

$C=\{x\in \mathbb {R} ^{n}:g_{i}(x)\leq 0\}$

Pode-se usar qualquer das funções de penalidade interior apresentadas, por exemplo:

$B(x)=\sum _{i=1}^{p}{\frac {1}{g_{i}(x)}}$

ou ainda

$B(x)=\sum _{i=1}^{p}-\ln(-g_{i}(x))$

Como $C^{0}=\{x\in \mathbb {R} ^{n}:g_{i}(x)<0\}$, o problema passa a ser:

$\left\{{\begin{matrix}minf(x)+B(x)\\x\in C^{0}\end{matrix}}\right.$

Observação: Note que é necessário conhecer um ponto inicial $x_{0}\in C^{0}$, para servir de primeira aproximação, ou ponto de partida para o algoritmo.

Métodos de região de confiança

Considere o seguinte problema de programação diferenciável não linear sem restrições:

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

onde $f:\mathbb {R} ^{n}\mapsto \mathbb {R}$ é de classe $C^{2}(\mathbb {R} ^{n})$.

Observação: Como $f$ não é necessariamente convexa, a matriz $\nabla ^{2}f(x)$ pode não ser definida positiva, apesar de ser simétrica. Neste caso, o método de Newton ou suas variantes (direções conjugadas, quase Newton, etc) não servem.

O primeiro método de região de confiança (em inglês, Trust region method), foi introduzido por Powel em 1970 (qual artigo?) mas oficialmente introduzido por Dennis em 1978 (artigo?). Ele consiste no seguinte:

A cada iteração, se constrói um modelo quadrático e uma região de confiança

este princípio pode ser considerado como uma extensão da busca de Armijo unidimensional.

Breve revisão da busca de Armijo

Para entender a geometria do método de região de confiança, é bom lembrar a geometria da busca de Armijo unidimensional.

Seja $d$ uma direção de descida de uma função $f$ a partir do ponto $x$. Então $\nabla f(x)^{\top }d<0$. Agora, considerando $\theta :\mathbb {R} \mapsto \mathbb {R}$, definida por $\theta (t)=f(x+td)$, tem-se:

$\theta '(t)=\nabla f(x+td)^{\top }d$. Assim, $\theta '(0)=\nabla f(x)^{\top }d<0$. Se $\alpha \in (0,1)$, então $\theta '(0)<\alpha \theta '(0)<0$.

Exercício

Mostrar que sempre existe $t>0$ tal que ${\frac {f(x+td)-f(x)}{t}}={\frac {\theta (t)-\theta (0)}{t-0}}\leq \alpha \theta '(0)$.

Resolução
A resolução deste exercício é deixada a cargo do leitor. Sinta-se livre para melhorar a qualidade deste texto, incluindo-a na versão online deste material.

A busca de Armijo consiste em tomar $t>0$.

Se ${\frac {f(x+td)-f(x)}{t}}>\alpha \theta '(0)>\theta '(0)$

Mas como $\theta '(0)=\lim _{t\to 0}{\frac {f(x+td)-f(x)}{t}}$, então existe algum $t>0$ tal que vale ${\frac {f(x+td)-f(x)}{t}}\leq \alpha \theta '(0)$. A tal ponto, chama-se ponto de Armijo.

Introduzindo as seguinte notação

$m(x+h)=f(x)+\nabla f(x)^{\top }h$	(modelo linear para a função $f$)
$pred=m(x)-m(x+td)=-t\theta '(0)$	(redução predita por este modelo)
$ared=f(x)-f(x+td)$	(redução real no valor da função)

a pergunta é:

Quando um ponto $t$ vai ser ponto de Armijo com essa notação?

Tem-se:

$\alpha \cdot pred=-\alpha t\theta '(0)\leq f(x)-f(x+td)=ared$

Logo, se $\alpha \cdot pred\leq ared$ segue que $t$ é um ponto de Armijo.

Observações: Note que a essência da busca linear de Armijo é construir um modelo linear e um intervalo compacto $[0,t]$, sendo $t>0$ e o ponto inicial da busca e logo procurar o ponto de Armijo em $[0,t]$.

De volta ao método

O método de região de confiança será uma generalização da busca de Armijo, consistindo da construção de um modelo quadrático e uma região $R$, chamada de região de confiança, e nessa região calcular o novo iterando.

Algoritmo da região de confiança

Primeiro passo: Escolha $x_{0}\in \mathbb {R} ^{n}$, $r>0$, $\alpha \in (0,1)$, $\beta \in (0,1)$ e $k=0$.

Passo iterativo $k$: Enquanto $\nabla f(x_{k})\not =0$, construa o modelo quadrático:
 $m_{k}(h)=\nabla f(x_{k})^{\top }h+{\frac {1}{2}}h^{\top }H_{k}h$
Calcule $d$, solução de
 $(P)\left\{{\begin{matrix}\min m_{k}(h)\\\|h\|\leq r\end{matrix}}\right.$
Tome $ared=f(x_{k})-f(x_{k}+d)$ e $pred=-m_{k}(d)$
Se ${\frac {ared}{pred}}>\alpha$, fazer $x_{k+1}=x_{k}+d$
Senão $x_{k+1}=x_{k}$, $r=\beta r$ e $k=k+1$

Comentários: No algoritmo anterior, quando se tem um passo falho, a região de confiança sempre diminui. Seria bom incluir casos bons, onde a região deve crescer.

Algoritmo da região de confiança melhorado

Primeiro passo: Escolha $x_{0}\in \mathbb {R} ^{n}$, $r>0$, $0<\alpha <\gamma <1$, $\beta \in (0,1)$, $\beta _{2}>1$, $\epsilon \in (0,1)$ e $k=0$.

Passo iterativo $k$: Enquanto $\|\nabla f(x_{k})\|>\epsilon$, construa o modelo quadrático:
 $m_{k}(h)=\nabla f(x_{k})^{\top }h+{\frac {1}{2}}h^{\top }H_{k}h$
Continue = 1
Enquanto (Continue = 1)
  Calcule $d$, solução de
   $(P)\left\{{\begin{matrix}\min m_{k}(h)\\\|h\|\leq r\end{matrix}}\right.$
  Tome $ared=f(x_{k})-f(x_{k}+d)$ e $pred=-m_{k}(d)$
  Se ${\frac {ared}{pred}}>\alpha$, fazer $x_{k+1}=x_{k}+d$
    Se $r>\gamma$, $r=\beta _{2}r$; Continue = 0
    k = k+1
  Senão $r=\beta _{1}r$

O subproblema quadrático

Conforme foi explicado, o método das regiões de confiança constrói um modelo quadrático da forma:

$m(h)=g^{\top }h+{\frac {1}{2}}h^{\top }Hh$

onde, no método, $g=\nabla f(x_{k})$ e $H=H_{k}$. Em tal modelo, tem-se

• $g$ um vetor não nulo;
• $h$ uma matriz simétrica (que pode não ser definida positiva)

O problema quadrático é

$(PQ)\left\{{\begin{matrix}\min m(h)\\\|h\|\leq r\end{matrix}}\right.$

com $r>0$.

Este é o problema que será tratado a seguir.

Exercício

Provar que $(PQ)$ sempre tem solução.

Resolução
Esboço: Como $h$ é uma função contínua e o conjunto onde se quer minimizar tal função é uma bola, em dimensão finita, trata-se de minimizar uma função contínua em um compacto. Esse é um problema que tem solução, pelo teorema de Weierstrass.

O exercício anterior garante que o método está bem definido, quer dizer, todas as etapas podem ser realizadas.

O problema de mínimos quadrados

Este tipo de problema é muito frequente em ciências experimentais. Para ter um exemplo em mente durante a discussão que será feita mais adiante, considere as seguintes informações:

Segundo dados disponibilizados pelo IBGE, o Produto Interno Bruto per capita na cidade de São Paulo, no período de 2002 a 2005, foi (em reais): 17734, 19669, 20943 e 24083, respectivamente.

Com base nessas informações, como poderia ser feita uma previsão do valor correspondente ao ano seguinte (2006)?

Uma escolha possível seria supor que a cada ano o PIB aumenta uma aproximadamente constante, ou seja, usar um modelo linear para obter tal estimativa (que possivelmente será bem grosseira). Intuitivamente, bastaria analisar os dados disponíveis e a partir deles deduzir qual é o aumento que ocorre a cada ano. Depois, a previsão para 2006 seria aproximadamente igual à de 2005 somada com aquele aumento anual.

Esta idéia poderia até funcionar para o caso deste exemplo, mas o que fazer se a quantidade de dados disponíveis sobre algum fenômeno (ou alguma situação) for significativamente maior?

A melhor escolha, sem dúvida, é fazer uso de um computador para obter o modelo que melhor descreve o "comportamento" dos dados experimentais.

Em geral, os problemas de mínimos quadrados consistem em identificar os valores de determinados parâmetros, de modo que se satisfaçam certas equações $h_{i}(x)=0,\,i=1,\ldots ,p$. No contexto do exemplo anterior, se procura um modelo linear para os dados, ou seja, uma função $y=mx+n$ que os descreva da melhor forma possível. Assim, os parâmetros a considerar são $m$ e $n$. Os valores ideais para essas variáveis seriam aqueles que verificassem as seguintes equações: