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 \ge 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 A x + b^\top A x + \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 A x + A^\top b$

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

$A^\top A x = -A^\top b$

Observe também que isso implica $A^\top (A x + b) = 0$, ou seja, $A x + 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^tDx > 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 \ge 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^2f(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{align} 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{align}$ 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^2f(\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)$.