Resumo de mínimos quadrados - parte 2

Em diversos textos anteriores abordei alguns temas sobre otimização e um texto que recebeu muita atenção foi o de Regressão Linear usando Mínimos Quadrados. No primeiro texto expliquei a motivação inicial da minimização de normas $l_2$ e apresentei alguns exemplos de regressões lineares usando o método de Mínimos quadrados. Neste texto irei um pouco mais a fundo no assunto e apresentarei a formulação básica do problema de otimização, sua solução direta e uma extensão deste método que visa controlar o erro de generalização.

Relembrando a Regressão Linear¶

O objetivo de uma regressão linear é estimar um vetor $y \in \mathbb{R}^n$ a partir das entradas de uma matrix $X \in \mathbb{R}^{n \times p}, p < n$ utilizando uma função linear $g(x) = x^t w + b$. Tipicamente, $X$ contém em suas linhas observações de um vetor aleatório $\mathbf{X}$ e $y$ contém observações de uma variável aleatória $\mathbf{y}$ que queremos estimar a partir de $\mathbf{X}$. Neste cenário podemos obter facilmente novas observações de $\mathbf{X}$, mas não de $\mathbf{y}$.

Como vimos anteriormente, estimamos o vetor de pesos $w \in \mathbb{R}^p$ minimizando a soma dos erros ao quadrado (Eq. \eqref{custo}).

\begin{equation} \min_w ||Xw - y||_2^2 \label{custo} \end{equation}

Primeiramente,

$$ ||Xw - y||_2^2 = (Xw - y)^T (Xw - y) = \\ (Xw)^TXw - (Xw)^T y - y^T Xw + y^T y = \\ w^T X^T X w - 2 w^T X^T y + ||y||_2^2. $$

Como $||y||_2^2$ é um termo constante (não depende de $w$), podemos ignorá-lo durante a minimização. Sabemos que, no ponto ótimo $w^*$, o gradiente do custo é 0. Logo,

$$ 2X^T X w^* - 2X^T y = 0 \\ 2X^T X w^* = 2X^T y \\ w^* = (X^T X)^{-1} X^T y $$

Portanto, a solução ótima pode ser determinada resolvendo um sistema de equações lineares se $(X^TX)$ for invertível, ou seja, se todas as colunas (ou linhas) desta matrix forem linearmente independentes (LI). Isto ocorre se as colunas/linhas de $X$ forem LI.

---

Prova:

Suponha, primeiramente, que as colunas de $X$ são LI mas $G = (X^TX)$ possui uma (ou mais) colunas LD e seja $j$ o índice desta coluna. Logo, o elemento $i, j$ de $G$ pode ser escrito como

$$ G_{i,j} = \sum_{k \neq j} \alpha_k G_{i, k} = \sum_{k \neq j} \alpha_k X_i^T X_k = X_i^T (\sum_{k \neq j} \alpha_k X_k) $$

Como $G_{i, j} = X_i^T X_j$, isto implica que $X_j = \sum_{k \neq j} \alpha_k X_k$, uma contradição pois as colunas de $X$ são LI. Logo, se as colunas de $X$ são LI entõ as colunas de $(X^T X)$ também devem ser LI.

Por outro lado, suponha que as colunas de $X$ são LD. Então existe um $j$ tal que $X_j = \sum_{k \neq j} \alpha_k X_k$. Vamos então analisar a coluna $G_j$

$$ G_j = \left( \begin{array}{c} X_1^T X_j \\ \dots \\ X_i^T X_j \\ \dots \end{array} \right) = \left( \begin{array}{c} X_1^T (\sum_{k \neq j} \alpha_k X_k) \\ \dots \\ X_i^T (\sum_{k \neq j} \alpha_k X_k) \\ \dots \end{array} \right) = \left( \begin{array}{c} \sum_{k \neq j} \alpha_k X_1^T X_k \\ \dots \\ \sum_{k \neq j} \alpha_k X_i^T X_k \\ \dots \end{array} \right) = \sum_{k \neq j} \alpha_k G_k $$

Logo, se a $X$ não possui rank completo, $X^T X$ também não. Portanto, uma regressão linear usando mínimos quadrados só possui solução única se a matrix $X$ possuir rank completo.

---

Apesar de este método poder ser usado para calcular a solução, outros métodos (como a decomposição QR) podem ser mais rápidos na prática e também numericamente mais estáveis, o que significa que eles cometem erros de aproximação menores.

Controlando overfitting¶

A solução tradicional usando OLS possui duas desvantagens que são especialmente custosas quando consideramos problemas com muitas variáveis.

O primeiro problema é que soluções OLS tendem a privilegiar o encaixe em pontos com erro grande. Dado que um aumento no número de variáveis implica em um aumento nas distâncias entre os pontos, um único outlier pode literalmente estragar o modelo predito. A segunda desvantagem é que a solução OLS nem sempre existe ou é única. Isto ocorre principalmente quando $p > n$, ou seja, quando temos mais variáveis que exemplos.

Estes dois pontos negativos podem ser mitigados adicionando um termo de regularização da solução. O termo mais comum (e matematicamente conveniente) de ser adicionado é a norma $\mathcal{l}_2$ ao quadrado (ou seja, $w^T w$). Este tipo de regularização é conhecido como Ridge Regression e é representado pela equação abaixo (Eq. \eqref{rigde}). Iremos mostrar que este problema sempre tem solução e argumentar que esta formulação privilegia soluções com coeficientes menores, evitando situações como a mostrada anteriormente.

\begin{equation} \min ||Xw - y||_2^2 + \lambda||w||_2^2 \label{ridge} \end{equation}

Como anteriormente, calculamos o gradiente da função de custo em um mínimo $w^*$:

$$ \nabla(w^*) ||Xw^* - y||_2^2 + \lambda||w^*||_2^2 = 0\\ 2X^T X w^* - 2X^T y + 2\lambda w^* = 0 (X^T X + \lambda I)w^* = X^T y. $$

Logo, \begin{equation} w^* = (X^T X + \lambda I)^{-1} X^T y \label{ridge-sol} \end{equation}

A matriz $(X^T X + \lambda I)$ é invertível se ela não possui nenhum autovalor $0$. Primeiramente, sabemos que $X^T X$ é positiva semi-definida, o que implica que todos seus autovalores são não negativos. Logo, tomemos um par de autovalor/vetor $(\alpha, w)$ de $X^T X$. Pela definição,

$$ (X^T X)w = \alpha w \\ (X^T X)w + \lambda w = \alpha w + \lambda w \\ (X^T X + \lambda I)w = (\alpha + \lambda) w $$

Logo, o par $(\alpha + \lambda, w)$ é autovalor/vetor de $X^T X + \lambda I$. Como $X^T X$ é PSD, $\alpha \geq 0$, logo $\alpha + \lambda > 0$, pois $\lambda >0$. Portanto, $X^T X + \lambda I$ é invertível.

---

O termo $||w||_2^2$ penaliza soluções que possuem componentes com grande magnitude. Isto significa que soluções com custo muito baixo, mas com um ou mais componentes muito maiores que os outros irão ter seu custo aumentado significativamente, enquanto soluções ligeiramente piores mas com magnitudes baixas não irão ter grande aumento no custo. O parâmetro $\lambda$ controla justamente quais tipos de soluções queremos privilegiar. Uma $\lambda$ grande privilegia soluções com baixa magnitude, enquanto um $\lambda$ pequeno interfere menos nos pesos estimados.

Vemos abaixo um exemplo

In [14]:

%matplotlib inline
import numpy as np
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import numpy.linalg

def poly2(inter):
    return np.array([np.ones(len(inter)), [x**2 for x in inter]]).T

def lstsq(X, y):
    inv1 = np.linalg.inv(np.dot(X.T, X))
    return np.dot(inv1, np.dot(X.T, y))

def ridge(X, y, lamb):
    inv1 = np.linalg.inv(np.dot(X.T, X) + lamb)
    return np.dot(inv1, np.dot(X.T, y))

def gen_data():
    inter = np.random.uniform(-50, 50, 21)
    X = poly2(inter)
    y = np.array([x**2 + 56 for x in inter])
    
    inter2 = inter.copy()
    y2 = y.copy()
    inter2[-1] = -15
    y2[-1] = -2000
    X2 = poly2(inter2)
    return inter, X, y, inter2, X2, y2

def plot_curve(inter, alpha, y):
    plt.plot(inter, y, 'o')
    ic = range(-50, 50)
    y_curve = np.dot(poly2(ic), alpha)
    plt.plot(ic, y_curve, 'r-')

inter, X, y, inter2, X2, y2 = gen_data()

plt.figure(figsize=(15,6))
plt.suptitle('Regressão linear')
plt.subplot(121)
alpha = lstsq(X, y)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = lstsq(X2, y2)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)

plt.figure(figsize=(15, 6))
plt.suptitle('Ridge com lambda=300')
plt.subplot(121)
alpha = ridge(X, y, 300)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = ridge(X2, y2, 300)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)


plt.figure(figsize=(15, 6))
plt.suptitle('Ridge com lambda=10')
plt.subplot(121)
alpha = ridge(X, y, 10)
error = mean_squared_error(y, np.dot(X, alpha))
plt.title('0 outliers: %s, erro %f'%(str(alpha), error))
plot_curve(inter, alpha, y)

plt.subplot(122)
alpha2 = ridge(X2, y2, 10)
error = mean_squared_error(y2, np.dot(X2, alpha2))
plt.title('1 outlier: %s, erro %f'%(str(alpha2), error))
plot_curve(inter2, alpha2, y2)

No caso sem outliers, a regressão linear "direta" obtém os melhores resultados, tanto visualmente quanto em termos do erro. No caso com outliers, o erro da regressão direta é menor, mas a curva está, visualmente, menos ajustada aos pontos do que as duas versões regularizadas. A versão com $\lambda$ grande "toca" todos os pontos da curva (menos o outlier) apesar de ter um erro quadrático médio pior e mesmo a versão com $\lambda$ pequeno é claramente melhor que a sem regularização. Este exemplo ilustra ambas as propriedades que comentamos antes e motiva a utilização de alguma regularização ao estimar os coeficientes de uma regressão linear. É recomendado, ao usar Ridge Regression, manter todas as variáveis em uma mesma escala (ou ao menos em uma escala parecida), já que a otimização dará preferência a soluções com norma menor. Um procedimento comum é, para cada variável $x_i$, aplicar a transformação $x_i = \frac{x_i - \overline{x_i}}{\sigma_{x_i}}$ para centralizar as variáveis e deixá-las com desvio padrão 1.

Isto conclui este novo post sobre regressão linear. Se houverem detalhes que não ficaram claros ou sugestões de outros assuntos (ou erros, por que não?), não hesite em deixar um comentário ;)

Até a próxima!