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,\dots ,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:

$\{\begin{array}{ccc}2002m+n& =& 17734\\ 2003m+n& =& 19669\\ 2004m+n& =& 20943\\ 2005m+n& =& 24083\end{array}$

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,\dots ,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:

$\underset{x\in {ℝ}^n}{\min }\sum _{i=1}^p{[h_i(x)]}^2$

O caso linear

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

$h_i:{ℝ}^n\mapsto ℝ$

$h_i(x)=a_i^{\mathrm{\top }}x+b_i$

onde $a_i\in {ℝ}^n$ e $b_i\in ℝ$ para cada $i=1,\dots ,p$. Deste modo, pode-se definir uma função $H$ como

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

$H(x)=\left[\begin{array}{c}{h}_1(x)\\ ⋮\\ h_p(x)\end{array}\right]=\left[\begin{array}{c}{a}_1^{\mathrm{\top }}x+b_1\\ ⋮\\ a_p^{\mathrm{\top }}x+b_p\end{array}\right]=\left[\begin{array}{c}{a}_1^{\mathrm{\top }}\\ ⋮\\ a_p^{\mathrm{\top }}\end{array}\right]x+\left[\begin{array}{c}{b}_1\\ ⋮\\ b_p\end{array}\right]$

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

$A=\left[\begin{array}{c}{a}_1^{\mathrm{\top }}\\ ⋮\\ a_p^{\mathrm{\top }}\end{array}\right],b=\left[\begin{array}{c}{b}_1\\ ⋮\\ b_p\end{array}\right]$

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:

$\underset{x\in {ℝ}^n}{\min }\frac{1}{2}\Vert Ax+b{\Vert }^2$

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

$f(x)=\frac{1}{2}\Vert Ax+b{\Vert }^2=\frac{1}{2}⟨Ax+b,Ax+b⟩=\frac{1}{2}{x}^{\mathrm{\top }}{A}^{\mathrm{\top }}Ax+b^{\mathrm{\top }}Ax+\frac{1}{2}{b}^{\mathrm{\top }}b$

Logo, como $B=A^{\mathrm{\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 {ℝ}^n$ tal que $\mathrm{\nabla }f(x)=0$ é solução do problema aproximado.

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

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

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

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

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

Exemplo de aplicação

Considere dado um conjunto de pontos do plano ${ℝ}^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?

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

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^{\mathrm{\top }}A$ sejam cada vez mais dispersos. Com isso, $A^{\mathrm{\top }}A$ torna-se mal condicionada.

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)\{\begin{array}{c}\min f(x)\\ x\in {ℝ}^n\end{array}$

sendo $h_j:{ℝ}^n\to ℝ$, e portanto $f=\frac{1}{2}\sum _{j=1}^p[h_j{]}^2$ de classe ${𝒞}^2\left({ℝ}^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 é $\overline{x}$, então:

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

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

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

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

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

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

donde

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

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

$x^{*}=\overline{x}-{[{\mathrm{\nabla }}^{2}f(\overline{x})]}^{-1}\mathrm{\nabla }f(\overline{x})$

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

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

Algoritmo de Newton (puro)

Início: Tome $x_0\in {ℝ}^n$
  Se $\mathrm{\nabla }f(x_0)=0$, pare: $x_0$ é ponto crítico.
  Senão, Calcule $d_0$, solução de
    ${\mathrm{\nabla }}^{2}f(x_0)d_0=-\mathrm{\nabla }f(x_0)$
    Faça $x_1=x_0+d_0$ e $k=1$
Iteração: Se $\mathrm{\nabla }f(x_k)=0$, pare: $x_k$ é ponto crítico.
  Senão, calcular $d_k$, solução de
    ${\mathrm{\nabla }}^{2}f(x_k)d_k=-\mathrm{\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{[h_i(x)]}^2$

Calculando o gradiente desta função, resulta:

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

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

$H(x)=\left[\begin{array}{c}{h}_1(x)\\ ⋮\\ h_p(x)\end{array}\right]$

Deste modo, a Jacobiana de $H$ verifica:

$\mathrm{\nabla }f(x)=[J_H(x)]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

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

Seja

$B(x)=\sum _{i=1}^p\mathrm{\nabla }{h}_i(x)\mathrm{\nabla }{h}_i(x{)}^{\mathrm{\top }}=J_H(x)J_H(x{)}^{\mathrm{\top }}$

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

$x^{\mathrm{\top }}\mathrm{\nabla }{h}_i(x)\mathrm{\nabla }{h}_i(x{)}^{\mathrm{\top }}x={[\mathrm{\nabla }{h}_i(x{)}^{\mathrm{\top }}x]}^2\ge 0$

Para que $D$ seja definida positiva, é necessário que $S_P(\mathrm{\nabla }{h}_i(x))={ℝ}^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 {ℝ}^n$
  Se $\mathrm{\nabla }f(x_0)=0$, pare: $x_0$ é ponto crítico.
  Senão, Calcule $d_0$, solução de
    $B(x_0)d_0=-\mathrm{\nabla }f(x_0)$, onde $B(x)=\sum _{i=1}^p\mathrm{\nabla }{h}_i(x)\mathrm{\nabla }{h}_i(x{)}^{\mathrm{\top }}$
    Faça $x_1=x_0+d_0$ e $k=1$
Iteração: Se $\mathrm{\nabla }f(x_k)=0$, pare: $x_k$ é ponto crítico.
  Senão, calcular $d_k$, solução de
    $B(x_k)d_k=-\mathrm{\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 {ℝ}^n$
  Se $\mathrm{\nabla }f(x_0)=0$, pare: $x_0$ é ponto crítico.
  Senão, Calcule $d_0$, solução de
    $(B(x_0)+\rho I)d_0=-\mathrm{\nabla }f(x_0)$, onde $B(x)=\sum _{i=1}^p\mathrm{\nabla }{h}_i(x)\mathrm{\nabla }{h}_i(x{)}^{\mathrm{\top }}$
    Faça $x_1=x_0+d_0$ e $k=1$
Iteração: Se $\mathrm{\nabla }f(x_k)=0$, pare: $x_k$ é ponto crítico.
  Senão, calcular $d_k$, solução de
    $(B(x_0)+\rho I)d_k=-\mathrm{\nabla }f(x_k)$
    Faça $x_{k+1}=x_k+d_k$ e $k=k+1$

Predefinição:AutoCat