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))

Duas linhas, então vale a cena por trás delas. np.ones((X.shape[0], 1)) cria uma coluna nova, uma casinha pra cada paciente, todas preenchidas com o número 1. np.hstack (de "horizontal stack") gruda essa coluna bem na frente de X: mesma quantidade de pacientes de sempre, só que agora com mais uma "variável" artificial, sempre igual a 1 pra todo mundo, sem exceção.

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] multiplica sempre a coluna de 1's (então na prática funciona como 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],)

Vale parar em cada linha do fit, mas trocando a lupa da matemática pela lupa da intuição. Lembra do post anterior que o gradiente descendente aprendia por tentativa e erro, feito um aluno chutando resposta, vendo o quanto errou, corrigindo um pouquinho, chutando de novo, centenas de vezes até o erro parar de cair? A equação normal não chuta nada, ela calcula a resposta certa direto:

• X = include_bias(X): cola a coluna de 1, a mesma função de acima.
• np.linalg.pinv(X): aqui mora o pulo do gato. Pensa em cada um dos 353 pacientes do treino como alguém dando um palpite de como os 11 pesos deveriam ser, "com esses pesos aí, minha previsão bateria certinha comigo". São 353 palpites competindo por só 11 números, então é impossível cravar todo mundo ao mesmo tempo, sempre sobra alguém insatisfeito. pinv é o juiz que olha pra todos os palpites de uma vez só e acha o conjunto de pesos que deixa a "reclamação total" (a soma dos erros ao quadrado, a mesma conta que o gradiente descendente vinha reduzindo aos poucos, chute por chute) a menor possível. A diferença é que ele já entrega esse veredito final numa conta só, sem precisar chutar e corrigir centenas de vezes.
• np.linalg.pinv(X) @ y: essa multiplicação é o veredito do juiz virando número de verdade, o vetor w_ com os 11 pesos escolhidos.

O predict faz o caminho inverso: aplica o mesmo include_bias (pro paciente novo entrar com a mesma "forma" dos que treinaram o modelo) e multiplica pelos pesos já prontos pra chegar na previsão.

Uma classe inteira de regressão linear, resolvida numa linha de treino, sem nenhum chute. 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.