Resumo de Otimização: Mínimos Quadrados

Atualização: veja a parte 2 deste texto, onde apresento um método de solução direta do problema e uma formulação regularizada chamada Ridge Regression que evita overfitting. O método de Mínimos Quadrados (Least Squares) é utilizado para resolver sistemas "super"-determinados (ou impossíveis) , i.e. em que existem mais restrições que variáveis desconhecidas. Apesar de não ser possível encontrar uma solução que respeite todas as restrições simultaneamente, é possível calcular uma aproximação. O método de Mínimos Quadrados busca uma solução que minimize a a soma dos quadrados das diferenças entre o valor predito pelo modelo e o valor desejado. 

Dados pares de observações $(x_j, y_j), 1 \leq j \leq m, x_j \in \mathbb{R}^n, m > n$, queremos encontrar uma função $g(x; \theta), \theta \in \mathbb{R}^k$ tal que $g(x_j; \theta) = y_j$ para uma escolha adequada dos parâmetros $\theta$. A função $g(x; \theta)$ é o modelo utilizado e precisa ser definido previamente.

Idealmente, o modelo $g(x; \theta)$ seria flexível o suficiente para representar a transformação desejada. Na prática, porém, o modelo $g(x; \theta)$ pode somente aproximar o valor correto $y_j$ a partir de $x_j$. Desta forma, podemos definir o erro de aproximação do par $(x_j, y_j)$ como

$$r_j(\theta) = g(x_j; \theta) - y_j$$

e minimizar a função objetivo

$$f(\theta) =  \frac{1}{2}\sum^m_j [r_j(\theta)]^2 = \frac{1}{2}\sum^m_j [  g(x_j; \theta) - y_j    ]^2.$$

Ao minimizar $f(\theta)$ encontramos o vetor de parâmetros $\theta$ que "melhor" ajusta o modelo $g(x; \theta)$ aos dados $(x_j, y_j)$. Neste caso, "melhor" é definido como o vetor $\theta$ que minimiza a soma do quadrado dos erros. Por isto o nome de método de Mínimos Quadrados.

Justificativa estatística

Suponha que os $r_j(\theta)$ são independentes e identicamente distribuídos com média 0, variância $\sigma^2$ e função densidade de probabilidade $g_\sigma$. Sob estas condições, a probabilidade de observar os $r_j(\theta)$ correspondentes ao conjunto de dados $(x_j, y_j)$, dadas a variância $\sigma$ e o vetor de parâmetros $\theta$ é:

$$p( (x, y) ; \theta ; \sigma) = \prod_j^m g_\sigma( r_j(\theta) )  = \prod_j^m g_\sigma( g(x_j; \theta) - y_i )$$

Fixadas as observações $(x_j, y_j)$ e $\sigma$, o valor "mais provável" de $\theta$ é o que maximiza $p( (x, y) ; \theta ; \sigma)$. Este estimador para $\theta$ é chamado de Estimador de Máxima Verossimilhança.

Suponha, então, que a distribuição dos $r_j$ é normal. Portanto,

$$g_\sigma(r) = \frac{1}{\sqrt{2\pi \sigma^2}} \mathrm{exp}(-\frac{r^2}{2\sigma^2}),$$

e então,

$$p( (x, y) ; \theta ; \sigma) =   \frac{1}{(2\pi\sigma^2)^{m / 2}} \mathrm{exp}(-\frac{1}{2\sigma^2}  \sum_j^m [g(x_j; \theta) - y_i )]^2  ).$$

Para maximizar $p( (x, y) ; \theta ; \sigma)$ (dado $\sigma^2$ fixo) é necessário minimizar $\sum_j^m [g(x_j; \theta) - y_i )]^2$. Portanto, a estimativa de $\theta$ obtida pelo método de Mínimos Quadrados equivale ao estimador de máxima verossimilhança de $\theta$ (supondo que os erros $r_j$ são i.i.d. com distribuição normal de média 0 e variância fixa $\sigma^2$).

Modelos Lineares

Diversos problemas podem ser modelados utilizando funções $g(x; \theta)$ lineares em $\theta$. Podemos, portanto, escrever os resíduos como $r_j(\theta) = x_j^T\theta - y_j$.

Seja $J \in \mathbb{R}^{m\times n}$ tal que

$$J = \begin{pmatrix} x_1^T \\ \dots \\ x_m^T\\ \end{pmatrix},$$

então o problema de mínimos quadrados se reduz a minimizar a função

$$f(\theta) = \frac{1}{2} || J\theta - y ||_2^2$$

Modelos lineares são vantajosos pois a função acima pode ser minimizada globalmente. Isto não é necessariamente verdade para modelos não lineares em $\theta$.

Note que a função $g(x; \theta)$ deve ser linear em $\theta$, mas não necessariamente em $x$. Logo, a função

$$g(x; \theta) = {x^{(1)}}^2  \theta_{(1)} + x^{(1)}x^{(2)}\theta_{(2)}$$

é um modelo válido para a estimação linear, pois ela equivale a construir um novo conjunto de dados $(x'_j, y_j)$, onde $x'_j = ({x^{(1)}_j}^2, x_j^{(1)}x_j^{(2)})$, e fazer a minimização linear "comum".

Exemplos usando Python

A biblioteca Scipy possui tanto a implementação do Método de Mínimos Quadrados para o caso linear quanto para o caso geral. Neste texto mostrarei como utilizar a função numpy.linalg.lstsq para estimar modelos lineares.

Iremos estimar duas funções diferentes. Ambas podem ser estimadas usando modelos lineares. A função $f_1$ é um polinômio de segundo grau e a $f_2$ é uma função de duas variáveis.

$$f_1(x) = 2x^2 + 5x + 7$$

$$f_2(x_1, x_2) = 0.3x^3 - 0.2y^2 + 3y - xy + 5$$

Estimação de polinômios 

 

Utilizaremos três modelos para estimar $f_1$:

$$g_1(x) = p_1x + p_2$$
$$g_2(x) = p_1x + p_2x^2 + p_3$$
$$g_3(x) = p_1x + p_2x^2 + p_3x^3 + p_4$$

 É esperado que o modelo $g_1$ possua erro grande, pois um polinômio de grau 1 não consegue aproximar um polinômio de grau 2. Por outro lado, tanto $g_2$ quanto $g_3$ deverão obter boas aproximações, sendo que em $g_3$ o componente $x^3$ deverá ser próximo de 0.

Primeiramente, iremos gerar nosso conjunto de observações de $f_1$. Em cada observação também adicionaremos um pequeno erro.

x_data = [1, 4, 7, 8, -2, 5, 6, 12, -3, -5, -10]
y_data = [2 * x*x + 5*x + 7 - (random.random() * 5 - 2.5) for x in x_data]

Agora, iremos definir a matrix $J$ para cada um dos modelos e aplicar a função lstsq para estimar os parâmetros de cada modelo.

m =  len(x_data)
x2 = np.array([x*x for x in x_data]).reshape((m, 1))
x3 = np.array([x*x*x for x in x_data]).reshape((m, 1))
x_data = np.reshape(x_data, (m, 1))
y_data = np.reshape(y_data, (m, 1))

ones = np.ones(m).reshape(m, 1)

J1 = np.hstack([x_data, ones])
J2 = np.hstack([x2, x_data, ones])
J3 = np.hstack([x3, x2, x_data, ones])

theta_g_1 = np.linalg.lstsq(J1, y_data )
theta_g_2 = np.linalg.lstsq(J2, y_data )
theta_g_3 = np.linalg.lstsq(J3, y_data )

Os resultados obtidos podem ser vistos nas figuras abaixo:

Como esperado, a estimação linear foi bem ruim e a quadrática deu resultados quase perfeitos. O polinômio de grau 3 também teve resultados bons e apresentou um coeficiente bem pequeno para o termo $x^3$.

Estimação (linear) de funções de várias variáveis

Utilizaremos 2 modelos para estimar a função $f_2$:

$$g_1(x, y) = \theta_1 x + \theta_2 y + \theta_3 xy + \theta_4$$
$$g_2(x, y) = \theta_1 x^2 + \theta_2 x y + \theta_3 y^2 + \theta_4 x + \theta_5 y + \theta_6$$
$$g_3(x, y) = \theta_1 x^3 + \theta_2 x^2 y + \theta_3 x y^2 + \theta_4 y^3 + \theta_5 x^2 + \theta_6y^2 + \theta_7xy + \theta_8x + \theta_9y + \theta_{10}$$  

Primeiramente, geramos as observações usando o código abaixo. Note que é necessário um número muito maior de pontos para estimar funções mais complexas.


data = np.array([(random.randint(-50, 50), random.randint(-50, 50)) for i in range(100)])
y_data = [0.3*x**3 - 0.2*y**2 + 3*y - x*y + 5 - (random.random() * 5 - 2.5) for (x, y) in data]

Após, geramos os vetores com as componentes de cada modelo e as matrizes $J_1, J2$ e $J_3$.


m = data.shape[0]
x3  = np.array([x**3 for (x, y) in data]).reshape((m, 1))
x2  = np.array([x**2 for (x, y) in data]).reshape((m, 1))
x1   = np.array([x for (x, y) in data]).reshape((m, 1))
y3  = np.array([y**3 for (x, y) in data]).reshape((m, 1))
y2  = np.array([y**2 for (x, y) in data]).reshape((m, 1))
y1   = np.array([y for (x, y) in data]).reshape((m, 1))
x2y = np.array([x**2 * y for (x, y) in data]).reshape((m, 1))
y2x = np.array([x * y**2 for (x, y) in data]).reshape((m, 1))
xy  = np.array([x*y for (x, y) in data]).reshape((m, 1))
ones= np.ones(m).reshape((m, 1))

J_1 = np.hstack([x1, y1, xy, ones])
J_2 = np.hstack([x2, xy, y2, x1, y1, ones])
J_3 = np.hstack([x3, x2y, y2x, y3, x2, y2, xy, x1, y1, ones])

theta_g_1, err1, _, _ = np.linalg.lstsq(J_1, y_data)
theta_g_2, err2, _, _ = np.linalg.lstsq(J_2, y_data)
theta_g_3, err3, _, _ = np.linalg.lstsq(J_3, y_data)

Os resultados obtidos foram os seguintes:

 Erro: 3.73790848e+09

 Erro: 3.64574253e+09

Erro:  204.59062094

Os primeiros dois modelos, mais simples, resultaram em erros grandes pois eles não possuem algumas componentes do modelo original. Já o terceiro modelo obteve resultados muito bons, porém ele contém uma quantidade muito maior de parâmetros que o modelo original. Isto ocorre pois o modelo original, em geral, não é conhecido e acabamos incluindo componentes irrelevantes para a estimação de $y_i$.

Uma desvantagem do método de Mínimos Quadrados é que a minimização feita dificilmente zera algum $\theta_i$, pois ele pode super ajustar a curva utilizando pequenos coeficientes em componentes que podem ser irrelevantes. Desta maneira, o método não é capaz de detectar quais parâmetros/características são relevantes para um dado problema.

Comentários finais

A estimação utilizando o método de Mínimos Quadrados otimiza os parâmetros de um modelo de modo a minimizar a soma dos erros ao quadrado. Este método encontra mínimos globais se o modelo utilizado for linear (nos parâmetros a serem estimados) e é vantajoso quando se conhece ou se tem alguma intuição sobre a forma analítica do modelo original dos dados.
Atualização: veja a parte 2 deste texto, onde apresento um método de solução direta do problema e uma formulação regularizada chamada Ridge Regression que evita overfitting.