Equação Normal: Resolvendo Regressão numa Conta Só

No post anterior eu deixei um gancho: dá pra resolver regressão linear numa conta só, sem ficar iterando gradiente. É exatamente isso que as aulas 2b e 2c mostram, e de quebra o professor compara o resultado contra mais cinco tipos de regressor diferentes.

Zerando a derivada em vez de andar até ela

O gradiente descendente funciona porque, a cada passo, ele anda um pouquinho na direção que reduz o erro. Mas pra regressão linear, a função de erro (a soma dos quadrados que eu já vinha usando) tem uma propriedade especial: ela é uma parábola em relação aos pesos, uma tigela lisa sem vales falsos. E uma tigela lisa tem um único ponto onde a derivada é exatamente zero, o fundo dela. Em vez de caminhar até lá passo a passo, dá pra calcular esse ponto direto.

O Bishop escreve o modelo linear numa forma mais geral que a minha, com uma matriz Φ (a matriz de design, cada linha é um paciente, cada coluna uma variável de entrada, mais uma coluna de 1's pro termo de bias). A função de erro é

$E_D(\mathbf{w}) = \frac{1}{2}\sum_{n=1}^{N}\{t_n - \mathbf{w}^{\mathsf{T}}\boldsymbol{\phi}(\mathbf{x}_n)\}^2$

Igualando o gradiente dessa conta a zero e isolando w, sobra

$\mathbf{w}_{\text{ML}} = (\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\mathbf{t}$

Essa é a equação normal. Uma matriz, duas multiplicações e uma inversão, e pronto, o ajuste ótimo sai direto, sem taxa de aprendizado, sem escolher quantas iterações, sem nenhum dos cuidados que o gradiente descendente exige.

O truque de embutir o bias

Repara que a fórmula acima não tem um + b separado, o bias tá dentro do próprio w (o Bishop chama de $w_0$). Isso só funciona porque ele define uma "função de base" falsa, $\phi_0(\mathbf{x}) = 1$, uma coluna inteira de 1's, então multiplicar por $w_0$ dá exatamente $w_0$ pra todo paciente, o mesmo efeito de somar um bias fixo.

O código faz isso na unha, colando uma coluna de 1's na frente de X:

def include_bias(X):
    return np.hstack((np.ones((X.shape[0], 1)), X))

Isso muda a arquitetura da classe: em vez de guardar coefs_ e intercept_ como duas coisas separadas (como no post anterior), agora existe um único vetor w_, onde w_[0] é o bias e o resto são os coeficientes de cada variável.

Três versões, uma comparação real

O notebook roda três variações da mesma classe LinearRegressor, todas no mesmo dataset de diabetes (agora usando o train_test_split de verdade do scikit-learn, então os números não batem exatamente com os do post anterior, que usava minha própria implementação na unha):

A equação normal ganha das duas, sem eu precisar escolher learning_rate nem max_iter. Não é coincidência: gradiente descendente é uma forma de aproximar essa mesma resposta iterando, e com iterações e taxa de aprendizado suficientes ele converge pro mesmo lugar. A equação normal só pula direto pro final.

A pseudo-inversa é o pinv

$(\Phi^{\mathsf{T}}\Phi)^{-1}\Phi^{\mathsf{T}}$ tem nome próprio: pseudo-inversa de Moore-Penrose, denotada $\Phi^\dagger$. É uma generalização de "matriz inversa" pra matrizes que não são quadradas (o que é sempre o caso aqui: $\Phi$ tem uma linha por paciente e uma coluna por variável, quase nunca são iguais). np.linalg.pinv(X) calcula exatamente essa conta, então a linha toda vira

class LinearRegressor(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        X = include_bias(X)
        self.w_ = np.linalg.pinv(X) @ y
        return self

    def predict(self, X):
        X = include_bias(X)
        return (X @ self.w_).reshape(X.shape[0],)

Uma classe inteira de regressão linear, resolvida numa linha. E o Bishop já avisa o risco: se duas colunas de entrada forem muito parecidas entre si (colinearidade), $\Phi^{\mathsf{T}}\Phi$ fica perto de singular e a conta fica numericamente instável. Isso vai voltar quando a matéria chegar em seleção de características.

Bate com o scikit-learn, na casa decimal

A prova de que a conta tá certa: rodei a mesma pseudo-inversa e o LinearRegression() pronto do scikit-learn, lado a lado.

regressor = LinearRegressor()  # a classe com pinv, acima
regressor.fit(X_train, y_train)
print(mean_squared_error(y_train, regressor.predict(X_train)))

from sklearn.linear_model import LinearRegression
sk_regressor = LinearRegression().fit(X_train, y_train)
print(mean_squared_error(y_train, sk_regressor.predict(X_train)))

Saída: minha versão: MSE treino 2868.5497028355776. LinearRegression do scikit-learn: MSE treino 2868.549702835577.

A diferença aparece só na última casa decimal, ruído de ponto flutuante, não de método. LinearRegression() faz exatamente essa conta por baixo do capô.

O zoológico de regressores

A aula 2c pega esse mesmo regressor da equação normal e mede o MSE também no conjunto de teste (os 20% que ficaram de fora do treino), e depois compara contra mais cinco tipos de modelo bem diferentes entre si, todos com os parâmetros padrão do scikit-learn:

Eu ainda não vou explicar como cada um desses modelos funciona por dentro (KNN, árvore de decisão e floresta aleatória têm aula própria mais na frente da matéria, e é lá que eu volto neles com calma). Mas dá pra tirar três lições só olhando a tabela:

1. MSE de treino baixo não significa nada sozinho. A árvore de decisão zerou o erro de treino (ela literalmente decorou cada paciente) e foi a pior de todas no teste. Isso é overfitting em estado puro, o mesmo fenômeno que o Bishop mostrou lá no capítulo 1 com o polinômio de grau 9.
2. O vencedor no teste não foi o modelo mais "exato" no treino. A floresta aleatória acerta o treino quase tão bem quanto a árvore cheia (2530.82, bem perto do zero absurdo da árvore), mas sem exagerar, e por isso generaliza melhor: 2785.98 no teste, o menor MSE de teste da tabela inteira. Uma floresta é várias árvores treinadas em pedaços diferentes dos dados, com a previsão final sendo a média de todas, e essa média cancela boa parte do exagero que cada árvore individual comete.
3. O SGDRegressor bateu a equação normal exata no teste (2863.35 contra 2900.19), mesmo com um MSE de treino um pouco pior. Não é coincidência: como eu já vi na playlist da especialização, o SGDRegressor do scikit-learn vem com regularização L2 ligada por padrão. Aqui essa regularização, sem eu pedir, acabou ajudando a generalizar melhor.

O LinearSVR ficou visivelmente pior que todo o resto, treino e teste, mas isso é mais sobre os hiperparâmetros padrão dele não servirem bem pra esse dataset do que sobre o método em si, história pra outro dia.

Fechando

Aplicação Prática

Uso o modelo da equação normal (o mesmo LinearRegressor com pinv de mais acima) e olho pro conjunto de teste, os 89 pacientes que o modelo nunca viu durante o ajuste, pra visualizar o que aquele MSE de 2900.19 realmente significa caso a caso.

regressor = LinearRegressor()
regressor.fit(X_train, y_train)
y_pred_test = regressor.predict(X_test)

Repara que a nuvem é bem parecida com a do post anterior (que era em cima do treino), sem ficar visivelmente pior no teste. Isso confirma numericamente o que a tabela já mostrou: MSE treino 2868.55 contra MSE teste 2900.19, uma diferença pequena. O modelo não decorou o treino, ele generalizou de verdade, só que sendo um modelo linear simples, ainda erra bastante caso a caso, o mesmo teto que eu já via no post anterior.