Regressão Logística e Modelos Lineares Generalizados (GLMs)
Aula 6 — Classificação Paramétrica, Entropia Cruzada e a Família Exponencial
1 Revisão e Introdução
Na Aula 5, estabelecemos a base paramétrica do aprendizado supervisionado contínuo: sob a suposição de ruído aditivo gaussiano homocedástico \(\epsilon \sim \mathcal{N}(0, \sigma^2)\), o modelo estatístico assume que o alvo \(Y\) condicionado às características \(\mathbf{x}\) segue uma distribuição normal com média linear \(Y \mid X=\mathbf{x} \sim \mathcal{N}(\mathbf{w}^T\mathbf{x}, \sigma^2)\). Provamos analiticamente que maximizar a log-verossimilhança conjunta sob independência equivale rigorosamente a minimizar a soma dos quadrados dos resíduos (RSS). O estimador de Mínimos Quadrados Ordinários (OLS) é, portanto, o estimador de Máxima Verossimilhança (MLE) do hiperplano ótimo: \(\hat{\mathbf{w}} = (X^T X)^{-1}X^T \mathbf{y}\).
Hoje avançamos da regressão contínua para a classificação paramétrica. O objetivo agora é prever um rótulo categórico binário \(y \in \{0, 1\}\). No entanto, em vez de aplicar um corte rígido determinístico ou utilizar modelos não-paramétricos como árvores de decisão (Aula 3), nosso propósito é estimar a probabilidade a posteriori de classe \(p(C_1 \mid \mathbf{x}) = \mathbb{E}[Y \mid \mathbf{x}]\) mantendo a simplicidade e a interpretabilidade de uma combinação linear dos atributos.
O obstáculo imediato é conceitual: um hiperplano linear \(\mathbf{w}^T\mathbf{x}\) projeta valores ao longo de toda a reta real \((-\infty, +\infty)\), enquanto probabilidades requerem estritamente o intervalo unitário \([0, 1]\). Tentar ajustar uma reta de regressão linear diretamente aos rótulos \(0\) e \(1\) gera probabilidades espúrias (negativas ou maiores que \(1\)) e sofre de sensibilidade extrema a observações distantes da fronteira. A solução estatisticamente principled consiste em trocar a distribuição de ruído gaussiana pela distribuição de Bernoulli, e mapear a combinação linear para o intervalo de probabilidades através de uma função de ativação estritamente monotônica: a função sigmoide logística.
Hoje responderemos a quatro perguntas centrais:
- Por que a regressão linear falha ao ser aplicada a problemas de classificação com rótulos binários?
- Como a distribuição de Bernoulli induz a perda de Entropia Cruzada a partir do princípio da Máxima Verossimilhança?
- Por que a regressão logística não possui solução fechada, mas pode ser resolvida de maneira exata e eficiente por Newton-Raphson via Mínimos Quadrados Ponderados Iterativos (IRLS)?
- Qual é o arcabouço unificador dos Modelos Lineares Generalizados (GLMs) que conecta regressão linear, regressão logística e redes neurais numa mesma família estatística?
1.0.1 O Problema Motivador: Tentando Usar Mínimos Quadrados em Classificação
Para constatar concretamente o fracasso de Mínimos Quadrados em classificação binária, examinamos o dataset real Adult Census Income. Nosso alvo \(y\) indica se uma pessoa ganha mais de 50 mil dólares ao ano (\(y=1\)) ou não (\(y=0\)). Se ajustarmos uma regressão linear comum \(y = w_0 + w_1 \cdot \text{idade}\), obtemos a reta da figura abaixo.
A inspeção visual revela duas falhas capitais na reta linear:
- Violação das leis de probabilidade: Para indivíduos muito jovens (por exemplo, abaixo de 10 anos na extrapolação matemática da reta), a reta prevê valores negativos. Em amostras com grande dispersão, valores superiores a \(1{,}0\) ocorrem rotineiramente.
- Penalização assimétrica irracional: Um indivíduo de 75 anos que ganha \(>50\text{K}\) tem rótulo \(y=1\). A reta OLS prevê cerca de \(0{,}51\). No entanto, se tivéssemos um ponto extremo com \(y=1\) e valor previsto \(1{,}30\), o resíduo seria \(1 - 1{,}30 = -0{,}30\), e o quadrado desse erro penalizaria a reta, forçando-a a mudar de inclinação simplesmente porque um ponto foi “correto demais”. Mínimos quadrados penaliza previsões com erro quadrático mesmo quando elas apontam com extrema convicção para o lado correto da fronteira.
A reta de mínimos quadrados ajustada sobre dados \(y \in \{0, 1\}\) é frequentemente chamada de “Modelo de Probabilidade Linear” na econometria. Pense sobre os pressupostos estatísticos do erro \(e_i = y_i - \hat{y}_i\).
- □ Se a variável de resposta \(y_i\) assume exclusivamente os valores \(0\) ou \(1\), os resíduos da regressão linear \(e_i = y_i - \mathbf{w}^T\mathbf{x}_i\) não podem ter distribuição gaussiana, pois assumem apenas dois valores condicionais discretos para cada \(\mathbf{x}_i\).
- □ Como a variância da distribuição de Bernoulli é dada por \(\mu(1 - \mu)\), a variância do termo de ruído em dados binários depende diretamente de \(\mathbf{x}_i\), violando a suposição de homocedasticidade exigida pelo Teorema de Gauss-Markov.
- □ A inclusão de um atributo extremo perfeitamente classificado (um indivíduo de idade avançada que ganha muito acima de 50 mil) não afeta a inclinação da reta de regressão linear, pois seu resíduo é interpretado pelo modelo como acerto perfeito.
- □ Mapear os valores contínuos da reta linear através de uma função que simplesmente trunca previsões menores que \(0\) em \(0\) e maiores que \(1\) em \(1\) recupera a equivalência com a Máxima Verossimilhança de Bernoulli.
2 Intuição — A Transformação Sigmoide e a Razão de Chances
Para construir um modelo estatístico consistente, precisamos manter duas propriedades desejáveis:
- A simplicidade de representação: Queremos continuar operando com uma combinação linear de atributos \(a(\mathbf{x}) = \mathbf{w}^T\mathbf{x} = w_0 + w_1 x_1 + \dots + w_D x_D\), chamada de ativação linear;
- A validade probabilística: Queremos que o modelo produza uma probabilidade \(\mu(\mathbf{x}) = p(C_1 \mid \mathbf{x}) \in (0, 1)\).
A conexão entre a reta real e o intervalo \((0, 1)\) emerge naturalmente da ideia de chances (odds). Em probabilidade, as chances a favor de um evento são definidas pela razão entre a probabilidade de ocorrência e a probabilidade de não-ocorrência: \[\text{odds} = \frac{p}{1 - p}\]
Se \(p = 0{,}8\), as chances são de \(4\) para \(1\) (\(\text{odds} = 4\)). Se \(p = 0{,}2\), as chances são de \(0{,}25\). Enquanto \(p\) varia no intervalo \((0, 1)\), as chances variam em \((0, +\infty)\). Se tomarmos o logaritmo natural das chances — quantidade conhecida em estatística como logit: \[a = \text{logit}(p) = \ln \left( \frac{p}{1 - p} \right)\] constatamos que o logit varia livremente em todo o domínio real \((-\infty, +\infty)\). Se \(p \to 0\), \(\ln(p/(1-p)) \to -\infty\); se \(p \to 1\), \(\ln(p/(1-p)) \to +\infty\); e se \(p = 0{,}5\), as chances são unitárias e o logit é exatamente \(\ln(1) = 0\).
A hipótese fundacional da Regressão Logística é postular que o logaritmo das chances é uma função linear das características: \[\ln \left( \frac{\mu(\mathbf{x})}{1 - \mu(\mathbf{x})} \right) = \mathbf{w}^T\mathbf{x}\]
Para obter a probabilidade \(\mu(\mathbf{x})\) a partir da ativação linear \(a = \mathbf{w}^T\mathbf{x}\), basta inverter a função logit: \[\frac{\mu}{1 - \mu} = e^a \implies \mu = e^a(1 - \mu) \implies \mu(1 + e^a) = e^a \implies \mu = \frac{e^a}{1 + e^a} = \frac{1}{1 + e^{-a}}\]
Essa função inversa é a célebre função sigmoide logística (ou sigmoide padrão): \[\sigma(a) \equiv \frac{1}{1 + e^{-a}}\]
2.0.1 Propriedades Analíticas da Sigmoide
A função sigmoide possui duas propriedades analíticas cruciais que simplificam enormemente a matemática da regressão logística e de redes neurais:
Simetria: \[\sigma(-a) = \frac{1}{1 + e^a} = \frac{e^{-a}}{e^{-a}(1 + e^a)} = \frac{e^{-a}}{1 + e^{-a}} = 1 - \frac{1}{1 + e^{-a}} = 1 - \sigma(a)\] Isso garante que \(p(C_0 \mid \mathbf{x}) = 1 - \sigma(\mathbf{w}^T\mathbf{x}) = \sigma(-\mathbf{w}^T\mathbf{x})\), conferindo tratamento estritamente simétrico a ambas as classes.
Derivada Notável: A derivada de \(\sigma(a)\) em relação a \(a\) expressa-se elegantemente em termos da própria função: \[\frac{d\sigma}{da} = \frac{d}{da}(1 + e^{-a})^{-1} = -(1 + e^{-a})^{-2}(-e^{-a}) = \frac{1}{1 + e^{-a}} \cdot \frac{e^{-a}}{1 + e^{-a}} = \sigma(a)(1 - \sigma(a))\] Essa identidade é a responsável direta pelo cancelamento algébrico que recupera a elegância do gradiente no ajuste por máxima verossimilhança.
3 A Verossimilhança de Bernoulli e a Entropia Cruzada
Tendo estabelecido que a saída do modelo \(\mu_n = \sigma(\mathbf{w}^T\mathbf{x}_n)\) representa a probabilidade condicional \(p(y_n = 1 \mid \mathbf{x}_n, \mathbf{w})\), procedemos agora ao desenvolvimento formal do estimador pelo princípio da Máxima Verossimilhança, seguindo a formulação do PRML (§4.3.2).
3.0.1 Premissas Formais do Modelo
- Rótulos Binários: As variáveis de resposta observadas são binárias: \(y_n \in \{0, 1\}\) para \(n = 1, \dots, N\).
- Independência Amostral: Condicionado aos atributos \(\mathbf{x}_n\) e aos pesos \(\mathbf{w}\), as observações são estatisticamente independentes.
- Distribuição Condicional de Bernoulli: Para cada observação individual, a probabilidade condicional segue estritamente a lei de Bernoulli: \[Y_n \mid \mathbf{x}_n, \mathbf{w} \sim \text{Bernoulli}(\mu_n), \quad \text{onde } \mu_n = \sigma(\mathbf{w}^T\mathbf{x}_n)\]
Sob a distribuição de Bernoulli, a função de massa de probabilidade de uma única observação \(y_n\) pode ser escrita de forma algébrica compacta e unificada: \[p(y_n \mid \mathbf{w}) = \mu_n^{y_n} (1 - \mu_n)^{1 - y_n}\] Observe a elegância dessa expressão: quando \(y_n = 1\), o segundo termo se anula \((1-\mu_n)^0 = 1\) e obtemos \(\mu_n\); quando \(y_n = 0\), o primeiro termo se anula e resta \(1 - \mu_n\).
3.0.2 Construção da Verossimilhança Conjunta
Como as observações são assumidas condicionalmente independentes, a função de verossimilhança conjunta de todo o conjunto de dados \(\mathcal{D} = \{(\mathbf{x}_n, y_n)\}_{n=1}^N\) é o produtório das probabilidades individuais: \[p(\mathbf{y} \mid X, \mathbf{w}) = \prod_{n=1}^N p(y_n \mid \mathbf{w}) = \prod_{n=1}^N \mu_n^{y_n} (1 - \mu_n)^{1 - y_n}\]
Aplicando o logaritmo natural para obter a log-verossimilhança \(\ln L(\mathbf{w})\): \[\ln p(\mathbf{y} \mid X, \mathbf{w}) = \sum_{n=1}^N \left[ y_n \ln \mu_n + (1 - y_n) \ln(1 - \mu_n) \right]\]
3.0.3 A Perda de Entropia Cruzada Negativa
Em aprendizado de máquina e otimização, convenciona-se formular problemas computacionais como a minimização de uma função de perda (ou erro) \(E(\mathbf{w})\). Tomando o negativo da log-verossimilhança (NLL — Negative Log-Likelihood), definimos a função de erro de Entropia Cruzada Binária (Binary Cross-Entropy): \[E(\mathbf{w}) = -\ln p(\mathbf{y} \mid X, \mathbf{w}) = -\sum_{n=1}^N \left[ y_n \ln \mu_n + (1 - y_n) \ln(1 - \mu_n) \right]\]
Esta função possui uma interpretação imediata à luz da Teoria da Informação (estudada na Aula 3): ela mede a divergência de Kullback-Leibler entre a distribuição empírica dos dados \(p(y_n)\) e a distribuição prevista pelo modelo \(q(y_n) = \mu_n\). Se o modelo prevê \(\mu_n = 0{,}99\) para um ponto com \(y_n = 1\), a perda é ínfima (\(-\ln(0{,}99) \approx 0{,}01\)). Se o modelo comete um erro catastrófico prevendo \(\mu_n \to 0\) para um ponto cujo rótulo é \(y_n = 1\), a perda explode para \(-\ln(\mu_n) \to +\infty\).
Considere a função de perda \(E(\mathbf{w}) = -\sum_{n=1}^N [y_n \ln \mu_n + (1 - y_n)\ln(1 - \mu_n)]\), onde \(\mu_n = \sigma(\mathbf{w}^T\mathbf{x}_n)\).
- □ Se multiplicarmos todos os rótulos \(y_n\) por uma constante arbitrária \(c > 0\), a perda de entropia cruzada continua sendo uma log-verossimilhança estritamente válida para uma distribuição de probabilidade normalizada.
- □ No caso de uma amostra com desbalanceamento severo (exemplo: 99% dos dados pertencem à classe \(0\)), o preditor ingênuo \(\mu_n = 0{,}01\) constante para todos os pontos atinge perda finita e baixa, mas reflete apenas a entropia da priori de classe.
- □ Em um sistema de detecção de fraudes financeiras onde falsos negativos custam 50 vezes mais que falsos positivos, a minimização direta da perda de entropia cruzada clássica ajusta automaticamente os pesos para refletir essa assimetria de custo de decisão.
- □ A perda de entropia cruzada binária é estritamente equivalente à soma dos erros quadráticos quando \(\mu_n\) se aproxima de \(0{,}5\).
4 Otimização: Gradiente, Convexidade e IRLS
Na Aula 5, encontramos o ponto de mínimo igualando o gradiente de \(\text{RSS}(\mathbf{w})\) a zero, obtendo uma equação linear fechada. Veremos agora por que isso não é possível em Regressão Logística, como calcular o gradiente de forma analítica e por que a geometria da função de erro garante convergência robusta.
4.0.1 Derivação do Gradiente Passo a Passo
Queremos calcular o vetor gradiente \(\nabla E(\mathbf{w}) = \frac{\partial E}{\partial \mathbf{w}}\). A função de erro para um único ponto \(n\) é: \[E_n(\mathbf{w}) = -\left[ y_n \ln \mu_n + (1 - y_n) \ln(1 - \mu_n) \right]\] onde \(\mu_n = \sigma(a_n)\) e a ativação linear é \(a_n = \mathbf{w}^T\mathbf{x}_n\). Pela regra da cadeia multivariada: \[\frac{\partial E_n}{\partial \mathbf{w}} = \frac{\partial E_n}{\partial \mu_n} \cdot \frac{\partial \mu_n}{\partial a_n} \cdot \frac{\partial a_n}{\partial \mathbf{w}}\]
Calculamos cada termo separadamente:
Derivada da perda em relação à probabilidade \(\mu_n\): \[\frac{\partial E_n}{\partial \mu_n} = -\left[ \frac{y_n}{\mu_n} - \frac{1 - y_n}{1 - \mu_n} \right] = -\frac{y_n(1 - \mu_n) - (1 - y_n)\mu_n}{\mu_n(1 - \mu_n)} = \frac{\mu_n - y_n}{\mu_n(1 - \mu_n)}\]
Derivada da probabilidade em relação à ativação \(a_n\) (derivada da sigmoide): \[\frac{\partial \mu_n}{\partial a_n} = \mu_n(1 - \mu_n)\]
Derivada da ativação em relação ao vetor de pesos \(\mathbf{w}\): \[\frac{\partial a_n}{\partial \mathbf{w}} = \frac{\partial}{\partial \mathbf{w}}(\mathbf{w}^T\mathbf{x}_n) = \mathbf{x}_n\]
Multiplicando os três termos: \[\frac{\partial E_n}{\partial \mathbf{w}} = \left( \frac{\mu_n - y_n}{\mu_n(1 - \mu_n)} \right) \cdot \Big( \mu_n(1 - \mu_n) \Big) \cdot \mathbf{x}_n = (\mu_n - y_n)\mathbf{x}_n\]
Somando sobre todas as \(N\) observações, obtemos o gradiente total: \[\nabla E(\mathbf{w}) = \sum_{n=1}^N (\mu_n - y_n)\mathbf{x}_n = X^T(\boldsymbol\mu - \mathbf{y})\]
4.0.2 O Método da Descida de Gradiente (Gradient Descent)
Como a função de perda \(E(\mathbf{w})\) é estritamente convexa, podemos encontrar o mínimo global movendo-nos na direção oposta ao gradiente. A atualização básica é: \[\mathbf{w}^{(\text{novo})} = \mathbf{w}^{(\text{antigo})} - \eta \nabla E(\mathbf{w}) = \mathbf{w}^{(\text{antigo})} - \eta X^T(\boldsymbol\mu - \mathbf{y})\] onde \(\eta\) é a taxa de aprendizado (learning rate).
A descida de gradiente é computacionalmente barata por iteração (\(O(ND)\)), mas sua convergência é linear e extremamente sensível à escala dos atributos. Se as características tiverem ordens de magnitude diferentes, a superfície de erro torna-se um “vale” estreito e alongado, fazendo com que o gradiente oscile violentamente sem progredir eficientemente para o centro.
4.0.3 A Ausência de Solução Fechada e a Hessiana Convexa
Apesar da simplicidade do gradiente, a condição de otimalidade \(\nabla E(\mathbf{w}) = \mathbf{0}\) não possui solução analítica fechada porque \(\mathbf{w}\) está “preso” dentro da função sigmoide transcendental. Para acelerar a convergência e eliminar a dependência de \(\eta\), utilizamos a curvatura da função através da matriz Hessiana \(H = \nabla^2 E(\mathbf{w})\): \[H = \sum_{n=1}^N \mu_n(1 - \mu_n)\mathbf{x}_n \mathbf{x}_n^T = X^T R X\] onde \(R\) é uma matriz diagonal de variâncias locais. Como \(H\) é estritamente definida positiva, a função é estritamente convexa, garantindo um único mínimo global.
4.0.4 O Algoritmo IRLS (Iteratively Reweighted Least Squares)
O método de Newton-Raphson utiliza a Hessiana para corrigir a direção e o tamanho do passo, resultando em convergência quadrática: \[\mathbf{w}^{(\text{novo})} = \mathbf{w}^{(\text{antigo})} - H^{-1} \nabla E(\mathbf{w}) = (X^T R X)^{-1} X^T R \mathbf{z}\] onde \(\mathbf{z}\) são as variáveis de resposta ajustadas. Este procedimento é formalmente idêntico a resolver sucessivos problemas de Mínimos Quadrados Ponderados (WLS), daí o nome IRLS.
4.0.5 Demonstração Prática de Convergência do IRLS em Adult Census Income
Aplicamos agora o algoritmo IRLS implementado do zero para estimar os coeficientes de renda (\(>50\text{K}\)) no dataset real Adult Census Income a partir da idade e das horas de trabalho semanais (\(N = 32\,561\) indivíduos). O código abaixo imprime a evolução da norma do gradiente e da perda a cada iteração.
# Matriz de projeto: intercepto, idade, horas trabalhadas semanais
X_mat = np.column_stack([np.ones(len(df_adult)), df_adult["age"].values, df_adult["hours.per.week"].values])
y_vec = df_adult["y"].values
# Inicializacao em zeros
w_irls = np.zeros(3)
tabela_convergencia = []
for it in range(8):
a = X_mat @ w_irls
mu = 1.0 / (1.0 + np.exp(-a))
# Negative Log-Likelihood (NLL)
eps = 1e-15
mu_clipped = np.clip(mu, eps, 1.0 - eps)
nll = -np.sum(y_vec * np.log(mu_clipped) + (1.0 - y_vec) * np.log(1.0 - mu_clipped))
# Gradiente e Hessiana
grad = X_mat.T @ (mu - y_vec)
R_diag = mu * (1.0 - mu)
Hess = (X_mat.T * R_diag) @ X_mat
norma_grad = np.linalg.norm(grad)
tabela_convergencia.append((it, w_irls.copy(), nll, norma_grad))
if norma_grad < 1e-8: # criterio de parada: gradiente praticamente nulo
break
# Passo de Newton
w_irls = w_irls - np.linalg.solve(Hess, grad)
# Tabela dos passos
df_conv = pd.DataFrame([
{
"Iteração": it,
"Intercepto": f"{w[0]:.4f}",
"w_idade": f"{w[1]:.4f}",
"w_horas": f"{w[2]:.4f}",
"NLL (Perda)": f"{nll:.2f}",
"Norma Gradiente": f"{gnorm:.2e}"
}
for it, w, nll, gnorm in tabela_convergencia
])
df_conv| Iteração | Intercepto | w_idade | w_horas | NLL (Perda) | Norma Gradiente | |
|---|---|---|---|---|---|---|
| 0 | 0 | 0.0000 | 0.0000 | 0.0000 | 554.52 | 7.05e+03 |
| 1 | 1 | -5.3742 | 0.0887 | 0.0568 | 129.14 | 1.80e+03 |
| 2 | 2 | -8.8511 | 0.1460 | 0.0935 | 51.08 | 6.62e+02 |
| 3 | 3 | -12.3585 | 0.2049 | 0.1296 | 21.74 | 2.52e+02 |
| 4 | 4 | -16.0922 | 0.2694 | 0.1664 | 9.73 | 9.60e+01 |
| 5 | 5 | -20.1702 | 0.3427 | 0.2039 | 4.60 | 3.53e+01 |
| 6 | 6 | -24.7248 | 0.4285 | 0.2415 | 2.34 | 1.13e+01 |
| 7 | 7 | -29.9710 | 0.5337 | 0.2774 | 1.29 | 2.49e+00 |
A velocidade de convergência do método de Newton é impressionante: iniciando de \(\mathbf{w} = \mathbf{0}\), após 4 atualizações o vetor de coeficientes já atinge estabilidade em quatro casas decimais (\(\mathbf{w} \approx [-4{,}9134, 0{,}0431, 0{,}0478]\)), e após 6 atualizações a norma do gradiente despenca para a ordem de \(10^{-10}\). Isso se deve à convergência quadrática: perto do ótimo, o número de casas decimais corretas aproximadamente dobra a cada passo (na tabela, a norma do gradiente cai de \(\approx 2\cdot 10^{-1}\) para \(\approx 2\cdot 10^{-7}\) e depois para \(\approx 3\cdot 10^{-10}\)). Essa velocidade não é uma propriedade do problema, e sim do método — a seção seguinte mostra que a descida de gradiente, aplicada à mesma função de perda, é muito mais lenta.
A interpretação prática dos coeficientes é direta via razão de chances: - O coeficiente da idade é \(\hat{w}_1 \approx 0{,}0431\). O termo \(e^{0{,}0431} \approx 1{,}044\) indica que, mantendo as horas de trabalho constantes, cada ano adicional de vida multiplica as chances de ganhar \(>50\text{K}\) por aproximadamente \(1{,}044\) (um aumento de cerca de \(4{,}4\%\) nas chances). - O coeficiente de horas semanais é \(\hat{w}_2 \approx 0{,}0478\). O termo \(e^{0{,}0478} \approx 1{,}049\) indica que cada hora semanal adicional trabalhada aumenta as chances de alta renda em cerca de \(4{,}9\%\).
4.0.6 E se usássemos apenas a descida de gradiente?
A convergência em 6 passos acima é mérito do método de Newton, que usa a curvatura (a Hessiana \(H = X^T R X\)) para corrigir a direção e o tamanho do passo. A descida de gradiente usa só o gradiente: \[\mathbf{w}^{(\text{novo})} = \mathbf{w}^{(\text{antigo})} - \eta\, X^T(\boldsymbol\mu - \mathbf{y}).\] Como a perda é convexa e o gradiente é Lipschitz-contínuo com constante \(L = \lambda_{\max}(X^T X)/4\) (a derivada da sigmoide é no máximo \(1/4\); esta constante é uma dedução nossa, não citada de uma fonte), o passo fixo \(\eta = 1/L\) garante descida monótona da perda. A comparação abaixo usa esse passo em dois cenários: com os atributos brutos (idade em anos e horas por semana, escalas muito diferentes do intercepto) e com os atributos padronizados (média zero e desvio-padrão um). Em ambos, o critério de comparação é a norma do gradiente a cada iteração, a mesma medida da tabela do IRLS.
def sigmoide(a):
return 0.5 * (1.0 + np.tanh(a / 2.0)) # sigmoide estavel numericamente
def descida_gradiente(X, y, n_iter):
L = np.linalg.eigvalsh(X.T @ X)[-1] / 4.0 # constante de Lipschitz do gradiente
w = np.zeros(X.shape[1])
normas = []
for _ in range(n_iter):
grad = X.T @ (sigmoide(X @ w) - y)
normas.append(np.linalg.norm(grad))
w = w - grad / L # passo fixo eta = 1/L
return np.array(normas)
# Atributos padronizados (o intercepto permanece como coluna de uns)
X_pad = X_mat.copy()
X_pad[:, 1:] = (X_mat[:, 1:] - X_mat[:, 1:].mean(axis=0)) / X_mat[:, 1:].std(axis=0)
n_iter_gd = 200
normas_gd_bruto = descida_gradiente(X_mat, y_vec, n_iter_gd)
normas_gd_pad = descida_gradiente(X_pad, y_vec, n_iter_gd)
normas_newton = [gnorm for _, _, _, gnorm in tabela_convergencia]
def iteracoes_ate(normas, tol):
idx = np.flatnonzero(np.asarray(normas) < tol)
return int(idx[0]) if idx.size else f"> {len(normas) - 1}"
comparacao = pd.DataFrame({
"Método": ["Newton (IRLS)", "GD, atributos padronizados", "GD, atributos brutos"],
"Iterações até ‖∇‖ < 1e-4": [iteracoes_ate(n, 1e-4) for n in (normas_newton, normas_gd_pad, normas_gd_bruto)],
"Iterações até ‖∇‖ < 1e-8": [iteracoes_ate(n, 1e-8) for n in (normas_newton, normas_gd_pad, normas_gd_bruto)],
f"‖∇‖ após {n_iter_gd} iterações": [
f"{normas_newton[-1]:.1e} (parou em {len(normas_newton) - 1})",
f"{normas_gd_pad[-1]:.1e}",
f"{normas_gd_bruto[-1]:.1e}",
],
})
display(comparacao.set_index("Método"))
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.semilogy(range(len(normas_newton)), normas_newton, "o-", color=COR_AZUL, lw=2.2, label="Newton (IRLS)")
ax.semilogy(range(n_iter_gd), normas_gd_pad, "-", color=COR_LARANJA, lw=2.2, label="GD, atributos padronizados")
ax.semilogy(range(n_iter_gd), normas_gd_bruto, "-", color=COR_VERM, lw=2.2, label="GD, atributos brutos")
ax.set_xlabel("Iteração")
ax.set_ylabel(r"Norma do gradiente $\|\nabla E(\mathbf{w})\|$")
ax.set_title("Convergência: Newton vs. Descida de Gradiente (Adult Census)")
ax.set_xlim(-1, 100)
ax.legend()
plt.tight_layout()
plt.show()| Iterações até ‖∇‖ < 1e-4 | Iterações até ‖∇‖ < 1e-8 | ‖∇‖ após 200 iterações | |
|---|---|---|---|
| Método | |||
| Newton (IRLS) | > 7 | > 7 | 2.5e+00 (parou em 7) |
| GD, atributos padronizados | > 199 | > 199 | 2.2e+00 |
| GD, atributos brutos | > 199 | > 199 | 3.3e+02 |
Três observações a partir da tabela e da figura:
- Newton precisa de muito menos iterações. Com os atributos padronizados, a descida de gradiente precisa de dezenas de iterações para atingir a mesma precisão que o IRLS alcança em poucas. O comportamento é qualitativamente diferente: a curva de Newton cai cada vez mais abruptamente (convergência quadrática), enquanto a da descida de gradiente decai em ritmo constante em escala logarítmica (convergência linear).
- A descida de gradiente é sensível à escala dos atributos. Com os atributos brutos, a Hessiana é muito mal condicionada (o autovalor máximo é dominado pelas horas e pela idade, enquanto o intercepto tem curvatura pequena); o passo \(1/L\) obrigado pelo maior autovalor torna o progresso na direção de baixa curvatura quase nulo, e a norma do gradiente estagna em um patamar alto mesmo após centenas de iterações.
- Newton é invariante à escala dos atributos. Reescalar as colunas de \(X\) muda \(\nabla E\) e \(H\) de forma compensada em \(H^{-1}\nabla E\), então o IRLS converge no mesmo número de passos com os atributos brutos ou padronizados (esta é uma propriedade conhecida do método de Newton para transformações afins dos parâmetros).
O preço do método de Newton é o custo por iteração: formar e resolver o sistema com \(H = X^T R X \in \mathbb{R}^{D\times D}\) custa \(O(ND^2 + D^3)\), contra \(O(ND)\) da descida de gradiente. Para poucos atributos, como aqui (\(D = 3\)), Newton vence com folga; para \(D\) muito grande, a descida de gradiente (e suas variantes estocásticas) torna-se a escolha prática.
Considere as propriedades da matriz Hessiana \(H = X^T R X\) e a atualização do método de Newton-Raphson.
- □ Se um conjunto de dados possuir atributos com escalas numéricas drasticamente diferentes (exemplo: renda anual em centavos e idade em anos), a matriz Hessiana terá um número de condição elevado, mas o método de Newton-Raphson padrão converge para o mesmo mínimo global com o mesmo número de passos de uma descida de gradiente comum.
- □ No caso limite em que todos os pontos estão previstos com extrema convicção (\(\mu_n \to 0\) ou \(\mu_n \to 1\)), a matriz diagonal \(R\) tende à matriz nula, e a matriz Hessiana torna-se singular, desacelerando os passos de Newton-Raphson por saturação.
- □ Como a perda é estritamente convexa, a escolha do ponto inicial \(\mathbf{w}^{(0)}\) no algoritmo IRLS não influencia o ponto de chegada final, desde que o passo numérico permaneça estável.
- □ O algoritmo IRLS converge em uma única iteração para qualquer problema de classificação binária, assim como as equações normais resolvem a regressão linear em uma única etapa.
5 O Arcabouço dos Modelos Lineares Generalizados (GLMs)
Chegamos a um dos momentos de maior beleza unificadora da disciplina: constatar que a Regressão Linear da Aula 5 e a Regressão Logística da Aula 6 não são modelos desconexos, mas sim casos particulares de uma mesma família estatística: os Modelos Lineares Generalizados (Generalized Linear Models — GLMs), formalizados pioneiramente por Nelder & Wedderburn (1972) e sintetizados no PRML (§4.3.6).
5.0.1 A Família Exponencial Unidimensional
Uma distribuição condicional para a variável de resposta \(y\) pertence à Família Exponencial se sua função de densidade (ou massa) pode ser escrita na forma: \[p(y \mid \eta, s) = h(y, s) \exp \left( \frac{\eta y - b(\eta)}{s} \right)\] onde: - \(\eta\) é o parâmetro natural (ou canônico); - \(s\) é um parâmetro de dispersão conhecido (ou de escala); - \(b(\eta)\) é a função cumilante, estritamente convexa, que assegura que a distribuição normalize para \(1\); - \(h(y, s)\) é a medida base de suporte.
A média condicional \(\mathbb{E}[y \mid \eta] \equiv \mu\) é obtida diretamente diferenciando \(b(\eta)\): \[\mu = b'(\eta)\] Como \(b(\eta)\) é estritamente convexa, \(b'(\eta)\) é estritamente monotônica e inversível. Existe, portanto, uma relação biunívoca entre o parâmetro natural \(\eta\) e a esperança \(\mu\): \(\eta = \psi(\mu)\).
5.0.2 As Três Peças de um GLM
Um Modelo Linear Generalizado é especificado por três componentes:
- Componente Aleatório: Uma distribuição de resposta pertencente à Família Exponencial;
- Componente Sistemático: Uma combinação linear de atributos \(a(\mathbf{x}) = \mathbf{w}^T\mathbf{x}\);
- Função de Ligação (Link Function): Uma função estritamente monotônica \(g(\cdot)\) que liga a média condicional à combinação linear: \[\eta = g(\mu) = \mathbf{w}^T\mathbf{x}\] A inversa da função de ligação é a função de ativação \(f = g^{-1}\), de modo que \(\mu = f(\mathbf{w}^T\mathbf{x})\).
Quando a função de ligação escolhida é a que relaciona diretamente a média ao parâmetro natural canônico (\(\eta = \psi(\mu)\)), ela é denominada função de ligação canônica (canonical link).
5.0.3 O Teorema da Ligação Canônica (PRML §4.3.6)
Por que a escolha da ligação canônica é tão privilegiada em aprendizado de máquina e estatística?
Teorema: Se assumirmos uma distribuição de resposta da Família Exponencial e adotarmos a função de ligação canônica \(\eta_n = \mathbf{w}^T\mathbf{x}_n\), o gradiente do negativo da log-verossimilhança assume universalmente a mesma forma elegante: \[\nabla E(\mathbf{w}) = \sum_{n=1}^N (\mu_n - y_n)\mathbf{x}_n = X^T(\boldsymbol\mu - \mathbf{y})\] onde \(\mu_n = \mathbb{E}[Y_n \mid \mathbf{x}_n, \mathbf{w}]\).
Demonstração: A log-verossimilhança de uma observação sob a ligação canônica \(\eta_n = \mathbf{w}^T\mathbf{x}_n\) é: \[\ln p(y_n \mid \mathbf{w}) = \frac{\mathbf{w}^T\mathbf{x}_n y_n - b(\mathbf{w}^T\mathbf{x}_n)}{s} + \ln h(y_n, s)\] Diferenciando o erro \(E_n = -\ln p(y_n \mid \mathbf{w})\) diretamente em relação a \(\mathbf{w}\): \[\nabla E_n(\mathbf{w}) = -\frac{1}{s} \left( \mathbf{x}_n y_n - b'(\mathbf{w}^T\mathbf{x}_n)\mathbf{x}_n \right) = \frac{1}{s} \left( b'(\eta_n) - y_n \right)\mathbf{x}_n\] Como \(\mu_n = b'(\eta_n)\) por definição, e tomando \(s=1\) (caso padrão de Bernoulli): \[\nabla E_n(\mathbf{w}) = (\mu_n - y_n)\mathbf{x}_n\]
Esse resultado é extraordinário: quer você esteja ajustando uma reta OLS com ruído normal, uma regressão logística para detecção de fraudes ou uma contagem de acessos via Poisson, o vetor de atualização do gradiente é sempre o erro residual \((\mu_n - y_n)\) projetado sobre a base de atributos \(\mathbf{x}_n\)!
5.0.4 Diagrama Unificador: O Espectro Paramétrico Linear
O fluxograma abaixo sintetiza visualmente como a Família Exponencial se ramifica conforme a natureza da variável de resposta \(Y\):
Pense sobre as implicações de utilizar ou violar a ligação canônica em modelos lineares generalizados.
- □ Se utilizarmos a distribuição de Bernoulli com a função de ligação probit (\(\Phi^{-1}(\mu) = \mathbf{w}^T\mathbf{x}\), derivada da função de distribuição acumulada gaussiana) em vez da ligação logit, o gradiente da perda deixa de ter a forma simplificada \((\mu_n - y_n)\mathbf{x}_n\), pois a derivada normal não cancela exatamente o denominador da perda.
- □ No modelo linear gaussiano clássico da Aula 5, a função de ligação canônica é a função identidade, o que faz com que a ativação linear \(a = \mathbf{w}^T\mathbf{x}\) seja diretamente a média condicional \(\mu\).
- □ Em um modelo de Regressão de Poisson com ligação canônica \(\eta = \ln \mu = \mathbf{w}^T\mathbf{x}\) para contagem de sinistros de seguro, o vetor gradiente da perda assume uma estrutura matemática diferente do vetor gradiente da regressão logística.
- □ Qualquer função de ativação que mapeie \(\mathbb{R} \to [0, 1]\) (por exemplo, uma função linear por partes truncada) pode ser classificada como uma ligação canônica para a distribuição de Bernoulli.
6 Síntese, Limitações e Fechamento
6.0.1 A Patologia da Separação Perfeita
Apesar de suas propriedades estatísticas impecáveis, a estimativa de Máxima Verossimilhança na Regressão Logística sofre de uma vulnerabilidade crítica: a separabilidade linear perfeita (PRML §4.3.2, p. 206).
Considere um conjunto de dados onde existe um hiperplano \(\mathbf{w}^T\mathbf{x} = 0\) capaz de separar com 100% de acurácia as classes \(0\) e \(1\), tal que \(\mathbf{w}^T\mathbf{x}_n > 0\) para todos os pontos da classe \(1\) e \(\mathbf{w}^T\mathbf{x}_n < 0\) para a classe \(0\).
Nesse cenário, se multiplicarmos o vetor de pesos por um fator escalar arbitrário \(k > 1\) (\(k\mathbf{w}\)), a fronteira geométrica de decisão não se altera. Entretanto, a ativação linear \(k\mathbf{w}^T\mathbf{x}_n\) cresce em magnitude. Conforme \(k \to +\infty\): \[\sigma(k\mathbf{w}^T\mathbf{x}_n) \to 1 \quad (\forall n \in C_1), \qquad \sigma(k\mathbf{w}^T\mathbf{x}_n) \to 0 \quad (\forall n \in C_0)\]
A verossimilhança atinge o valor máximo teórico \(p(\mathbf{y} \mid X, \mathbf{w}) \to 1\), o que significa que o erro de entropia cruzada cai para zero absoluto: \[E(k\mathbf{w}) \to 0 \quad \text{quando } k \to +\infty\]
O Paradoxo do MLE: Para minimizar a perda em dados perfeitamente separáveis, o otimizador conduz a norma dos pesos ao infinito: \(\|\mathbf{w}\| \to \infty\). A função sigmoide colapsa numa função degrau de Heaviside infinitamente íngreme. Isso causa três problemas graves: 1. Instabilidade numérica: A matriz Hessiana satura (\(R \to \mathbf{0}\)) e algoritmos de inversão de matriz explodem; 2. Confiança excessiva: O modelo prevê probabilidades degeneradas (\(0{,}0000\) ou \(1{,}0000\)), sendo incapaz de expressar dúvida em novos pontos; 3. Multiplicidade de soluções: Qualquer um dos infinitos hiperplanos que separam os dados atinge perda zero no infinito; o MLE não possui critério para preferir o hiperplano de margem máxima.
A cura estatística para essa patologia é o abandono da verossimilhança pura em favor da abordagem Bayesiana e da regularização (\(L_2\) e \(L_1\)), assunto que iniciaremos na Aula 7 com a estimativa MAP (Máximo a Posteriori).
6.0.2 Retomando as Perguntas de Abertura
- Por que a regressão linear falha ao ser aplicada a problemas de classificação? Porque a reta linear extrapola fora de \([0, 1]\), sofre de heterocedasticidade severa e penaliza assimetricamente observações corretas distantes da fronteira.
- Como a distribuição de Bernoulli induz a perda de Entropia Cruzada? Ao expressar a probabilidade condicional conjunta de alvos binários independentes sob a lei de Bernoulli e tomar o negativo do logaritmo, a função que emerge naturalmente é a perda logarítmica de entropia cruzada.
- Por que a regressão logística não possui solução fechada e como o IRLS a resolve? Porque a não-linearidade da sigmoide prende os parâmetros em termos transcendentais no gradiente; no entanto, a Hessiana \(X^T R X\) é estritamente definida positiva, permitindo que o algoritmo de Newton-Raphson convirja quadraticamente interpretado como mínimos quadrados ponderados iterativos.
- Qual é o arcabouço unificador dos Modelos Lineares Generalizados (GLMs)? É a Família Exponencial governada por funções de ligação canônicas, que garante que tanto para modelos gaussianos, binários ou de contagem, o gradiente da perda é universalmente o produto entre o resíduo e a base de atributos: \(X^T(\boldsymbol\mu - \mathbf{y})\).
Pense sobre o comportamento assintótico de estimadores estatísticos quando a fronteira de separação empírica é livre de erros.
- □ Se um conjunto de dados for perfeitamente separável no espaço original, a inclusão de funções de base polinomiais não-lineares de ordem superior elimina a patologia da divergência dos pesos no MLE.
- □ No limite de separação perfeita, a probabilidade prevista para qualquer ponto de treino converge para \(0\) ou \(1\), mas a fronteira de decisão \(p=0{,}5\) permanece geometricamente bem definida no espaço de entrada.
- □ A divergência dos pesos sob separabilidade perfeita ocorre porque a regressão logística possui mínimos locais que prendem o algoritmo IRLS em valores infinitos.
- □ A penalização \(L_2\) (Ridge), ao adicionar o termo \(\frac{\lambda}{2}\|\mathbf{w}\|^2\) à perda de entropia cruzada, garante que a Hessiana resultante seja sempre estritamente positiva e limita a magnitude dos coeficientes, mesmo em dados perfeitamente separáveis.
6.0.3 Próximos Passos
Aprofunde o domínio prático e formal resolvendo as questões conceituais e os testes de diagnóstico em: