Otimização/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:

$\left\{{\begin{matrix}2002m+n&=&17734\\2003m+n&=&19669\\2004m+n&=&20943\\2005m+n&=&24083\end{matrix}}\right.$

Note que a partir dessas equações poderiam ser definidas as funções $h_{i}$ como:

$h_{1}(m,n)=2002m+n-17734$

$h_{2}(m,n)=2003m+n-19669$

$h_{3}(m,n)=2004m+n-20943$

$h_{4}(m,n)=2005m+n-24083$

É de se esperar que o sistema de equações obtido a pouco não admitirá uma solução exata, pois se tem mais equações do que variáveis.

Isso geralmente acontece, pois é comum haver uma quantidade $p$ de equações bem maior que o número $n$ de parâmetros a identificar. Em particular, quando todas as equações $h_{i}$ são lineares em suas $n$ variáveis, dificilmente existirá uma solução exata para o sistema linear resultante, pois este terá mais equações do que incógnitas (como no exemplo). Em geral, não é possível encontrar parâmetros que satisfaçam exatamente todas as equações. Por isso, costuma-se tentar identificar os parâmetros que "melhor se aproximam" de uma solução exata, em algum sentido.

Uma forma de obter uma solução aproximada (uma "quase-solução") resulta da seguinte observação: o valor de cada função $h_{i}$ em uma solução exata $x$ deveria ser zero. Se tal exigência é restritiva demais, e com ela não é possível encontrar qualquer solução, uma possibilidade seria exigir um pouco menos. Por exemplo, poderia ser exigido apenas que o valor de $h_{i}$, para $i=1,\ldots ,p$ seja, em geral, pequeno. Uma das formas de capturar essa idéia em termos mais precisos é dizer que se pretende minimizar a soma dos quadrados dos valores de cada $h_{i}$. Em símbolos, o problema passaria a ser:

$\min _{x\in \mathbb {R} ^{n}}\sum _{i=1}^{p}\left[h_{i}(x)\right]^{2}$

O caso linear

Neste caso, para cada índice $i=1,\ldots ,p$, a função $h_{i}$ é afim linear, ou seja:

$h_{i}:\mathbb {R} ^{n}\mapsto \mathbb {R}$

$h_{i}(x)=a_{i}^{\top }x+b_{i}$

onde $a_{i}\in \mathbb {R} ^{n}$ e $b_{i}\in \mathbb {R}$ para cada $i=1,\ldots ,p$. Deste modo, pode-se definir uma função $H$ como

$H:\mathbb {R} ^{n}\mapsto \mathbb {R} ^{p}$ de modo que

$H(x)={\begin{bmatrix}h_{1}(x)\\\vdots \\h_{p}(x)\end{bmatrix}}={\begin{bmatrix}a_{1}^{\top }x+b_{1}\\\vdots \\a_{p}^{\top }x+b_{p}\end{bmatrix}}={\begin{bmatrix}a_{1}^{\top }\\\vdots \\a_{p}^{\top }\end{bmatrix}}x+{\begin{bmatrix}b_{1}\\\vdots \\b_{p}\end{bmatrix}}$

Motivando a introdução da seguinte notação:

$A={\begin{bmatrix}a_{1}^{\top }\\\vdots \\a_{p}^{\top }\end{bmatrix}},\quad b={\begin{bmatrix}b_{1}\\\vdots \\b_{p}\end{bmatrix}}$

Assim,

$H(x)=Ax+b$

Logo, buscar uma solução exata é o mesmo que procurar uma solução para o seguinte sistema linear:

$Ax=-b$

E uma solução aproximada poderia ser buscada a partir do seguinte problema de minimização:

$\min _{x\in \mathbb {R} ^{n}}{\frac {1}{2}}\|Ax+b\|^{2}$

Exercício

Verifique que se houver uma solução exata, ela também será um minimizador de $\sum _{i=1}^{p}\left[h_{i}(x)\right]^{2}$. A recíproca é verdadeira?

Resolução

Primeiramente, observe que para qualquer $x$ tem-se $\sum _{i=1}^{p}\left[h_{i}({\bar {x}})\right]^{2}\geq 0$ Isso significa que $0$ é uma cota inferior para os valores assumidos por esta soma de quadrados. Além disso, se ${\bar {x}}$ é uma solução exata, então para cada índice $i$ tem-se: $h_{i}({\bar {x}})=0$ e consequentemente $\sum _{i=1}^{p}\left[h_{i}({\bar {x}})\right]^{2}=\sum _{i=1}^{p}0^{2}=0$ Então, em ${\bar {x}}$ a soma dos quadrados assume seu valor mínimo. Não vale a recíproca, pois se para um certo ${\tilde {x}}$ a soma dos quadrados assumir um valor positivo $\epsilon$, então $\sum _{i=1}^{p}\left[h_{i}({\bar {x}})\right]^{2}=\epsilon >0$ implica que existe algum $i$ tal que $h_{i}({\tilde {x}})>0$ Isto significa que ${\tilde {x}}$ não satisfaz a equação de índice $i$, e portanto não pode ser uma solução exata.

Analisando a função objetivo do problema de minimização anterior, tem-se:

$f(x)={\frac {1}{2}}\|Ax+b\|^{2}={\frac {1}{2}}\langle Ax+b,Ax+b\rangle ={\frac {1}{2}}x^{\top }A^{\top }Ax+b^{\top }Ax+{\frac {1}{2}}b^{\top }b$

Logo, como $B=A^{\top }A$ é simétrica e semi-definida positiva, tem-se $f$ convexa. Isso implica que a condição necessária de primeira ordem é também suficiente. Assim, qualquer ponto $x\in \mathbb {R} ^{n}$ tal que $\nabla f(x)=0$ é solução do problema aproximado.

Calculando o gradiente da função objetivo tem-se:

$0=\nabla f(x)=A^{\top }Ax+A^{\top }b$

Deste modo, a solução do problema é obtida resolvendo o sistema

$A^{\top }Ax=-A^{\top }b$

Observe também que isso implica $A^{\top }(Ax+b)=0$, ou seja, $Ax+b\in \ker(A^{\top })$.

Exemplo de aplicação

Considere dado um conjunto de pontos do plano $\mathbb {R} ^{2}$, por exemplo $\{(x_{i},y_{i})\}_{i=1}^{p}$, representando dados obtidos experimentalmente.

Perguntas:

1. Qual é a função afim linear que melhor se aproxima dos dados experimentais?

Resolução

As funções afim lineares no plano são aquelas que possuem a forma $mx+b$. Então os parâmetros a identificar são $m$ e $b$. Como são as funções $h_{i}$? Lembrando que tais funções teriam imagem igual a zero em uma solução exata, pode-se tomar $h_{i}(m,b)=mx_{i}+b-y_{i}$ Assim, $H(x)={\begin{bmatrix}h_{1}(x)\\\vdots \\h_{p}(x)\end{bmatrix}}={\begin{bmatrix}mx_{1}+b-y_{1}\\\vdots \\mx_{p}+b-y_{p}\end{bmatrix}}={\begin{bmatrix}mx_{1}+b\\\vdots \\mx_{p}+b\end{bmatrix}}+{\begin{bmatrix}-y_{1}\\\vdots \\-y_{p}\end{bmatrix}}={\begin{bmatrix}x_{1}&1\\\vdots &\vdots \\x_{p}&1\end{bmatrix}}{\begin{bmatrix}m\\b\end{bmatrix}}+{\begin{bmatrix}-y_{1}\\\vdots \\-y_{p}\end{bmatrix}}$ A expressão final pode ser simplificada adotando a seguinte notação $A={\begin{bmatrix}x_{1}&1\\\vdots &\vdots \\x_{p}&1\end{bmatrix}},\quad x={\begin{bmatrix}m\\c\end{bmatrix}},\quad c={\begin{bmatrix}-y_{1}\\\vdots \\-y_{p}\end{bmatrix}}$ De fato, nesses termos tem-se $H(x)=Ax+c$

2. Qual é a função quadrática que melhor se aproxima dos dados experimentais?

Resolução

As funções quadráticas possuem a forma $ax^{2}+bx+c$. Então os parâmetros a identificar são $a$, $b$ e $c$. Pode-se tomar $h_{i}(a,b,c)=ax_{i}^{2}+bx_{i}+c-y_{i}$. Procedendo como antes, tem-se $H(x)=Ax+c$ onde $A={\begin{bmatrix}x_{1}^{2}&x_{1}&1\\\vdots &\vdots &\vdots \\x_{p}^{2}&x_{p}&1\end{bmatrix}},\quad x={\begin{bmatrix}a\\b\\c\end{bmatrix}},\quad c={\begin{bmatrix}-y_{1}\\\vdots \\-y_{p}\end{bmatrix}}$

Exercício

Explorar o exemplo anterior com dados concretos, implementando um dos três métodos de gradientes conjugados.

Resolução
A resolução deste exercício é deixada a cargo do leitor. Sinta-se livre para melhorar a qualidade deste texto, incluindo-a neste módulo.

Observações

Conforme se aumenta o grau do polinômio que faz a aproximação dos dados, as colunas de $A$ têm elementos elevados a potências cada vez maiores, fazendo com que os autovalores de $A^{\top }A$ sejam cada vez mais dispersos. Com isso, $A^{\top }A$ torna-se mal condicionada.

Exercício

É possível, usando o problema de mínimos quadrados, verificar o postulado de Euclides que diz que "por dois pontos distintos passa uma única reta"?

Resolução
A resolução deste exercício é deixada a cargo do leitor. Sinta-se livre para melhorar a qualidade deste texto, incluindo-a neste módulo.

Exercício

Usando o problema de mínimos quadrados, é possível inferir que "por três pontos distintos passa uma única curva quadrática"?

Resolução
A resolução deste exercício é deixada a cargo do leitor. Sinta-se livre para melhorar a qualidade deste texto, incluindo-a neste módulo.

O caso não linear

Para esse tipo de problemas, há dois métodos:

1. Gauss-Newton
2. Levemberg-Marquardt

Ambos são métodos do tipo Newton. Então, para entender cada um deles é preciso entender o Método de Newton.

O método de Newton

Para entender a essência do método de Newton, primeiro considere que o problema a ser resolvido é

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

sendo $h_{j}:\mathbb {R} ^{n}\rightarrow \mathbb {R}$, e portanto $f={\frac {1}{2}}\sum _{j=1}^{p}[h_{j}]^{2}$ de classe ${\mathcal {C}}^{2}\left(\mathbb {R} ^{n}\right)$. A idéia de Newton é usar o desenvolvimento até segunda ordem da série de Taylor da função $f$ em cada ponto iterado. Isto é, se o iterado é ${\bar {x}}$, então:

$Q(x)=f({\bar {x}})+\nabla f({\bar {x}})^{\top }(x-{\bar {x}})+{\frac {1}{2}}(x-{\bar {x}})\nabla ^{2}f({\bar {x}})(x-{\bar {x}})$

Então a condição de Newton é que em cada iteração a Hessiana $\nabla ^{2}f$ deve ser definida positiva.

Calculando o gradiente da função $Q$, segue:

$\nabla Q(x)=\nabla f({\bar {x}})+\nabla ^{2}f({\bar {x}})(x-{\bar {x}})$

Se $x^{*}$ é o (único) minimizador de $Q$, então

$0=\nabla Q(x^{*})=\nabla f({\bar {x}})+\nabla ^{2}f({\bar {x}})(x^{*}-{\bar {x}})$

donde

$\nabla ^{2}f({\bar {x}})(x^{*}-{\bar {x}})=-\nabla f({\bar {x}})$

Sendo $\nabla ^{2}f$ definida positiva, tal matriz é também inversível. Portanto:

$x^{*}={\bar {x}}-\left[\nabla ^{2}f({\bar {x}})\right]^{-1}\nabla f({\bar {x}})$

Assim, pode-se usar a seguinte iteração:

$x_{k+1}=x_{k}-\left[\nabla ^{2}f(x_{k})\right]^{-1}\nabla f(x_{k})$

Algoritmo de Newton (puro)

Início: Tome $x_{0}\in \mathbb {R} ^{n}$
  Se $\nabla f(x_{0})=0$, pare: $x_{0}$ é ponto crítico.
  Senão, Calcule $d_{0}$, solução de
    $\nabla ^{2}f(x_{0})d_{0}=-\nabla f(x_{0})$
    Faça $x_{1}=x_{0}+d_{0}$ e $k=1$
Iteração: Se $\nabla f(x_{k})=0$, pare: $x_{k}$ é ponto crítico.
  Senão, calcular $d_{k}$, solução de
    $\nabla ^{2}f(x_{k})d_{k}=-\nabla f(x_{k})$
    Faça $x_{k+1}=x_{k}+d_{k}$ e $k=k+1$

Voltando ao problema original, de mínimos quadrados, se tinha:

$f(x)={\frac {1}{2}}\sum _{i=1}^{p}\left[h_{i}(x)\right]^{2}$

Calculando o gradiente desta função, resulta:

$\nabla f(x)=\sum _{i=1}^{p}h_{i}(x)\nabla h_{i}(x)$

Considera-se $H:\mathbb {R} ^{n}\mapsto \mathbb {R} ^{p}$ definida por

$H(x)={\begin{bmatrix}h_{1}(x)\\\vdots \\h_{p}(x)\end{bmatrix}}$

Deste modo, a Jacobiana de $H$ verifica:

$\nabla f(x)=\left[J_{H}(x)\right]H(x)$

pois o produto de uma matriz por um vetor tem como resultado um vetor que é a combinação linear das colunas da matriz, com coeficientes que são as coordenadas do vetor.

Além disso, tem-se

$\nabla ^{2}f(x)=\sum _{i=1}^{p}h_{i}(x)\nabla ^{2}h_{i}(x)+\sum _{i=1}^{p}\nabla h_{i}(x)\nabla h_{i}(x)^{\top }$

Seja

$B(x)=\sum _{i=1}^{p}\nabla h_{i}(x)\nabla h_{i}(x)^{\top }=J_{H}(x)J_{H}(x)^{\top }$

Sabe-se que uma matriz $D$ é definida positiva se $x^{t}Dx>0$ para qualquer $x\not =0$. Fazendo $D=\nabla h_{i}(x)\nabla h_{i}(x)^{\top }$, tem-se:

$x^{\top }\nabla h_{i}(x)\nabla h_{i}(x)^{\top }x=\left[\nabla h_{i}(x)^{\top }x\right]^{2}\geq 0$

Para que $D$ seja definida positiva, é necessário que $S_{P}\left(\nabla h_{i}(x)\right)=\mathbb {R} ^{n}$ ( deve gerar todo o espaço), neste caso, se diz que $J_{H}(x)$ é de posto máximo.

Algoritmo de Gauss-Newton

Início: Tome $x_{0}\in \mathbb {R} ^{n}$
  Se $\nabla f(x_{0})=0$, pare: $x_{0}$ é ponto crítico.
  Senão, Calcule $d_{0}$, solução de
    $B(x_{0})d_{0}=-\nabla f(x_{0})$, onde $B(x)=\sum _{i=1}^{p}\nabla h_{i}(x)\nabla h_{i}(x)^{\top }$
    Faça $x_{1}=x_{0}+d_{0}$ e $k=1$
Iteração: Se $\nabla f(x_{k})=0$, pare: $x_{k}$ é ponto crítico.
  Senão, calcular $d_{k}$, solução de
    $B(x_{k})d_{k}=-\nabla f(x_{k})$
    Faça $x_{k+1}=x_{k}+d_{k}$ e $k=k+1$

Algoritmo de Levemberg-Marquardt

A idéia de Levemberg-Marquardt foi perturbar a matriz $B(x)$, considerando $B(x)+\rho I$, para algum $\rho >0$ pequeno.

Início: Tome $x_{0}\in \mathbb {R} ^{n}$
  Se $\nabla f(x_{0})=0$, pare: $x_{0}$ é ponto crítico.
  Senão, Calcule $d_{0}$, solução de
    $\left(B(x_{0})+\rho I\right)d_{0}=-\nabla f(x_{0})$, onde $B(x)=\sum _{i=1}^{p}\nabla h_{i}(x)\nabla h_{i}(x)^{\top }$
    Faça $x_{1}=x_{0}+d_{0}$ e $k=1$
Iteração: Se $\nabla f(x_{k})=0$, pare: $x_{k}$ é ponto crítico.
  Senão, calcular $d_{k}$, solução de
    $\left(B(x_{0})+\rho I\right)d_{k}=-\nabla f(x_{k})$
    Faça $x_{k+1}=x_{k}+d_{k}$ e $k=k+1$

Exercício

Enumere as diferenças que existem nos seguintes algoritmos:

• Newton
• Gauss-Newton
• Levenberg-Marquardt

Resolução

Newton Desenvolvendo $f$ em série de Taylor e exigindo que a matriz hessiana de $f$ seja definida positiva o método de Newton usa $d_{k}=-[\nabla ^{2}f(x_{k})]^{-1}\nabla f(x_{k})$. Gauss-Newton Em alguns casos a matriz hessiana de $f$ pode ser bem complicada. Nesses casos então a idéia do método de Gauss-Newton é usar uma boa aproximação para essa matriz. Podemos escrever $f(x)={\frac {1}{2}}\sum _{i=1}^{p}[h_{i}(x)]^{2}={\frac {1}{2}}\|h(x)\|^{2}={\frac {1}{2}}h(x)^{\top }h(x)$. Usando Taylor em $h$ temos: $h({\bar {x}})=h(x)+J(x)({\bar {x}}-x)$ onde $J(x)$ é a jacobiana de $h$. Substituindo $h({\bar {x}})$ na $f$ temos ${\begin{aligned}f({\bar {x}})&={\frac {1}{2}}h({\bar {x}})^{\top }h({\bar {x}})\\&={\frac {1}{2}}(h(x)+J(x)({\bar {x}}-x))^{\top }(h(x)+J(x)({\bar {x}}-x))\\&={\frac {1}{2}}[(h(x)^{\top }h(x)+2h(x)^{\top }(J(x)({\bar {x}}-x))+(J(x)({\bar {x}}-x))^{\top }(J(x)({\bar {x}}-x))]\\&=f(x)+h(x)^{\top }(J(x)({\bar {x}}-x))+(J(x)({\bar {x}}-x))^{\top }(J(x)({\bar {x}}-x))\end{aligned}}$ Calculando o gradiente e depois a segunda derivada temos: $\nabla f({\bar {x}})=J(x)^{\top }h(x)+J(x)^{\top }(J(x)({\bar {x}}-x))$ e $\nabla ^{2}f({\bar {x}})=J(x)^{\top }J(x)$ Denotando $B(x)=J(x)^{\top }J(x)$, então o método de Gauss-Newton escolhe $d_{k}$ tal que $B(x_{k})d_{k}=-\nabla f(x_{k})$. Levenberg-Marquardt A idéia do método de Levenberg-Marquardt é pertubar a matriz $B(x)$ um pouco, considerando $B(x)+\rho I$ com $\rho >0$ pequeno. Portanto o método de Levenberg-Marquardt escolhe $d_{k}$ tal que $(B(x_{k})+\rho I)d_{k}=-\nabla f(x_{k})$.