Curso de ML EN

Capítulo 16 de 37 · intermedio

Elastic net

Qué cubre este capítulo

Elastic net es regresión lineal con dos penalizaciones montadas a la vez: la penalización L1 de lasso, que empuja coeficientes a exactamente cero y te da un modelo disperso, y la penalización L2 de ridge, que mantiene los coeficientes pequeños y bien portados cuando tus features se mueven juntas. Una sola perilla, la proporción de mezcla, se desliza entre ambas. Gírala por completo hacia un lado y tienes lasso; por completo hacia el otro y tienes ridge; en cualquier punto intermedio obtienes los dos efectos al mismo tiempo.

Ese "los dos a la vez" es toda la razón de su existencia. La dispersión de lasso es maravillosa justo hasta que tus features están correlacionadas, momento en el que lasso se pone nervioso: agarra una feature de un grupo correlacionado, pone en cero las demás, y cuál agarró puede cambiar si remuestreas los datos. Ridge es estable exactamente en esa situación pero nunca te da un cero, así que te quedas con todas tus features. Elastic net es el default al que recurro cuando tengo un montón de features correlacionadas y quiero dispersión sin que lasso elija favoritas al azar. Se construye directamente sobre el capítulo de regresión lineal y se ubica entre los capítulos de ridge y lasso como su síntesis; si ya leíste esos, esta es la perilla que los conecta.

Lo construimos desde cero en NumPy — el objetivo de elastic net, descenso por coordenadas con un soft-threshold para la parte L1 y un factor de encogimiento para la parte L2, predict, R² — y lo comparamos contra sklearn.linear_model.ElasticNet. Los datos reales son el dataset de diabetes: 442 pacientes, diez features estandarizadas, y un grupo de colesterol sérico lo bastante correlacionado para que el efecto de agrupamiento se vea a simple vista.

Un poco de historia

Las dos penalizaciones llegaron primero, con décadas de diferencia. La regresión ridge es la más vieja: Arthur Hoerl y Robert Kennard la publicaron en 1970 como un arreglo para la multicolinealidad, agregando una penalización sobre los coeficientes al cuadrado para que XXX^{\top}X volviera a ser invertible y se amortiguaran los bandazos salvajes de coeficientes que producen las features correlacionadas. Encoge, pero nunca selecciona: cada coeficiente se queda en el modelo, solo que más pequeño. Luego, en 1996, Robert Tibshirani publicó el lasso, cambiando la penalización cuadrática por una de valor absoluto, y ese único cambio compró la dispersión: la penalización de valor absoluto tiene una esquina en cero, y la optimización desliza los coeficientes justo hacia esa esquina, así que el modelo hace selección de features gratis.

Lasso era el emocionante, y la gente encontró sus límites rápido. Puede seleccionar a lo más nn features cuando tienes más features que muestras, lo cual es un problema exactamente en los datasets anchos donde más quieres selección. Y con un grupo de features correlacionadas tiende a quedarse con una y descartar el resto de forma más o menos arbitraria, lo que vuelve inestable el conjunto seleccionado. Hui Zou y Trevor Hastie le pusieron nombre y solución en 2005, en "Regularization and variable selection via the elastic net": suma ambas penalizaciones, obtén dispersión de la parte L1 y selección agrupada y estable de la parte L2. Acuñaron el efecto de agrupamiento como la propiedad de que las features correlacionadas reciben coeficientes similares en lugar de que una sobreviva y las otras mueran. Cinco años después Friedman, Hastie y Tibshirani le dieron el solver rápido de descenso por coordenadas (el algoritmo glmnet) que scikit-learn sigue usando, y del que abajo construimos una versión pequeña.

La intuición

Empieza con el problema para el que elastic net fue construido. Aquí están dos de las features de diabetes, s1 y s2 — colesterol sérico total y LDL, dos mediciones de sangre que suben y bajan juntas. Graficadas una contra la otra son casi una línea: su correlación es 0.897.

Ahora piensa en qué hace un ajuste con dos features así de cercanas. Cargan casi la misma información, así que el modelo puede darle el crédito a s1, o a s2, o repartirlo de cualquier manera intermedia, y el error de entrenamiento apenas lo nota. Mínimos cuadrados ordinarios elige una de esas infinitas reparticiones y te da coeficientes enormes de signos opuestos que se cancelan: inestable e ilegible. Lasso elige una esquina: se queda con una de ellas y pone en cero la otra, lo cual al menos es disperso pero tira una feature real y depende del ruido para decidir cuál muere. Ridge se queda con ambas y las hace pequeñas y similares, lo cual es estable pero nunca disperso.

La jugada de elastic net es hacer ambas cosas al mismo tiempo. La parte L2 jala las features correlacionadas una hacia la otra para que compartan el coeficiente como lo hace ridge — ese es el efecto de agrupamiento — mientras la parte L1 sigue libre de poner en cero las features que no aportan nada. Obtienes un modelo disperso donde las features correlacionadas entran y salen como grupo, no de una en una por volado. Esa es la imagen que hay que retener: la proporción de mezcla decide cuánto "mantén juntas las features correlacionadas" compras al costo de cuánto "haz el modelo disperso".

La matemática

Todo es el modelo lineal del capítulo de regresión — las predicciones son una suma ponderada de features más un intercepto:

y^i=wxi+b,y^=Xw+b\hat{y}_i = w^{\top} x_i + b, \qquad \hat{y} = Xw + b

con wRdw \in \mathbb{R}^{d} el vector de coeficientes y bb el intercepto. Lo que cambia es la cosa que minimizamos. Mínimos cuadrados ordinarios minimiza solo el error; elastic net agrega una penalización sobre el tamaño de los coeficientes:

L(w,b)=12ni=1n(y^iyi)2+λ[αw1+1α2w22]L(w, b) = \frac{1}{2n}\sum_{i=1}^{n}\left(\hat{y}_i - y_i\right)^2 + \lambda\left[\alpha\lVert w\rVert_1 + \frac{1-\alpha}{2}\lVert w\rVert_2^2\right]

Dos perillas. La fuerza de penalización λ0\lambda \ge 0 fija qué tan duro empuja la penalización completa; en λ=0\lambda = 0 se desvanece y estás de vuelta en mínimos cuadrados. La proporción de mezcla α[0,1]\alpha \in [0, 1] reparte ese empuje entre las dos normas — la norma L1 w1=jwj\lVert w\rVert_1 = \sum_j |w_j| y la norma L2 al cuadrado w22=jwj2\lVert w\rVert_2^2 = \sum_j w_j^2. El 12n\frac{1}{2n} sobre el error y el 12\frac{1}{2} sobre el término L2 son conveniencias para que la derivada de abajo salga limpia; también empatan exactamente con el objetivo de scikit-learn, así que los dos ajustes son comparables coeficiente por coeficiente.

La proporción de mezcla es toda la historia. Pon α=1\alpha = 1 y el término L2 desaparece por completo:

L(w,b)=12ni(y^iyi)2+λw1L(w, b) = \frac{1}{2n}\sum_i\left(\hat{y}_i - y_i\right)^2 + \lambda\lVert w\rVert_1

que es el lasso. Pon α=0\alpha = 0 y el término L1 desaparece:

L(w,b)=12ni(y^iyi)2+λ2w22L(w, b) = \frac{1}{2n}\sum_i\left(\hat{y}_i - y_i\right)^2 + \frac{\lambda}{2}\lVert w\rVert_2^2

que es la regresión ridge. Elastic net es la combinación convexa de las dos, y esos son sus extremos.

Para ajustarlo no podemos simplemente dar un paso de gradiente, porque la norma L1 no tiene derivada en cero — esa esquina es justamente lo que produce la dispersión, y también lo que rompe el descenso por gradiente de a pie. El truco es el descenso por coordenadas: congela todos los coeficientes menos uno y minimiza la pérdida a lo largo de esa única coordenada, que es un problema unidimensional con solución de forma cerrada. Trabajar la derivada para el coeficiente jj, con features estandarizadas para que la curvatura 1nixij2\frac{1}{n}\sum_i x_{ij}^2 valga 1, da la actualización

wjS ⁣(ρj,  λα)1+λ(1α),ρj=1nxj ⁣(yy^+xjwj)w_j \leftarrow \frac{S\!\left(\rho_j,\; \lambda\alpha\right)}{1 + \lambda(1-\alpha)}, \qquad \rho_j = \frac{1}{n}\, x_j^{\top}\!\left(y - \hat{y} + x_j w_j\right)

donde ρj\rho_j es la correlación entre la feature jj y el residual con la propia contribución de jj sumada de vuelta, y SS es el operador soft-threshold:

S(z,γ)=sign(z)max ⁣(zγ,  0)S(z, \gamma) = \operatorname{sign}(z)\,\max\!\left(|z| - \gamma,\; 0\right)

Lee esa actualización como las dos penalizaciones actuando en secuencia. El soft-threshold en el numerador es la parte L1: encoge ρj\rho_j hacia cero por λα\lambda\alpha y fija en exactamente cero cualquier cosa más pequeña que eso — esto es lo que crea la dispersión. El denominador 1+λ(1α)1 + \lambda(1-\alpha) es la parte L2: un factor constante que encoge cada coeficiente hacia cero sin llegar nunca a él — esto es lo que estabiliza los grupos correlacionados. Cicla esa actualización sobre todos los coeficientes hasta que dejen de moverse, y ya ajustaste elastic net.

En qué es bueno, en qué no

El argumento a favor de elastic net es que te da la dispersión de lasso sin la fragilidad de lasso. Cuando las features están correlacionadas — y en datos reales normalmente lo están — la selección de lasso es inestable: el grupo de features que conserva puede cambiar cuando remuestreas, porque la penalización L1 no tiene razón alguna para preferir a un miembro de un grupo correlacionado sobre otro. La parte L2 arregla eso amarrando los coeficientes correlacionados entre sí, así que obtienes un modelo disperso cuyo conjunto seleccionado de verdad se queda quieto. También maneja el caso ancho, más features que muestras, donde el lasso puro se topa con un tope de nn features seleccionadas y elastic net no. Dos perillas en lugar de una es el precio, y validación cruzada sobre una rejilla de λ\lambda y α\alpha es la forma estándar de pagarlo.

El argumento en contra es que la regularización no es comida gratis y elastic net no siempre se gana su lugar. En estos datos de diabetes, con solo diez features y muestras de sobra, la penalización apenas mueve la exactitud predictiva — verás el R² de prueba aterrizar a mediados de los treintas, regularices o no, porque aquí no hay overfitting que arreglar. La recompensa no es un mejor score; es un modelo cuyos coeficientes puedes confiar y leer. Y las dos perillas interactúan de maneras que no siempre son intuitivas: el mismo λ\lambda significa cantidades muy distintas de dispersión a distintos α\alpha, así que no puedes afinar una y olvidarte de la otra. Recurre a él cuando tengas muchas features correlacionadas y quieras un modelo disperso y estable. No recurras a él esperando una exactitud que de otro modo no tenías.

Los datos

El dataset de diabetes viene incluido con scikit-learn, así que no hay nada que descargar. Tiene 442 pacientes, diez features de línea base — edad, sexo, índice de masa corporal, presión arterial y seis mediciones de suero sanguíneo etiquetadas s1 a s6 — y un target que mide la progresión de la enfermedad un año después de la línea base. Las features ya vienen estandarizadas: cada columna tiene media cero y está escalada, que es exactamente lo que un modelo penalizado necesita, ya que una penalización sobre el tamaño de los coeficientes solo tiene sentido cuando las features están en la misma escala.

La feature más fuerte por sí sola es el índice de masa corporal, y la relación es la que adivinarías — a mayor BMI, progresión más rápida, con mucha dispersión alrededor. Aquí está contra el target:

Lo interesante no es el BMI, sin embargo — es el grupo de suero. Las seis features s son distintas lecturas de química sanguínea, varias de ellas fuertemente correlacionadas, y s1 y s2 en particular se sientan en ese 0.897 que viste arriba. Ese cúmulo correlacionado es lo que hace de este un buen dataset para elastic net: es precisamente donde lasso y elastic net toman caminos separados.

Constrúyelo, una función a la vez

NumPy puro, en el orden en que lo escribirías: el modelo y su score primero, luego las dos piezas de la actualización, luego el loop que las ejecuta. Se asume que las features entran estandarizadas, que es lo que permite que la actualización por coordenada se mantenga así de simple.

El modelo y el score no cambian respecto a la regresión lineal ordinaria — la regularización cambia cómo elegimos los coeficientes, nunca cómo los usamos. La predicción es las features por los pesos más el intercepto:

def predict(X, w, b):
    """The linear model: yhat = Xw + b.

    Same one-liner as ordinary least squares. Elastic net changes how we CHOOSE
    w and b, never how we use them. X is (n, d), w is (d,), b is a scalar.
    """
    return X @ w + b

Y R² es la fracción de varianza explicada, el número que reportaremos — uno menos nuestro error cuadrático sobre el error de siempre adivinar la media:

def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of variance explained.

    1 minus (our squared error / the error of always guessing the mean). 1.0 is
    perfect, 0.0 is no better than the mean, negative is worse than the mean.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot

Ahora la pieza que hace a elastic net diferente de ridge: el operador soft-threshold. Esta es la penalización L1 en código. Encoge su entrada hacia cero por un margen y fija en exactamente cero cualquier cosa dentro de ese margen. La penalización suave de ridge solo puede escalar un coeficiente hacia abajo; esta puede apagarlo:

def soft_threshold(z, gamma):
    """The L1 operator: shrink z toward zero by gamma, and clamp at zero.

    S(z, gamma) = sign(z) * max(|z| - gamma, 0). This is the whole reason
    elastic net (and lasso) produce exact zeros: if the unpenalized pull on a
    coefficient is weaker than gamma, the coefficient is set to precisely 0, not
    to something small. Ridge's smooth L2 penalty can't do this — it scales
    coefficients down but never crosses zero.
    """
    return np.sign(z) * max(abs(z) - gamma, 0.0)

Vale la pena escribir el objetivo completo, tanto porque es lo que estamos minimizando como porque es la forma honesta de verificar un ajuste — la pérdida tiene que bajar en cada barrido. Es el error cuadrático medio (a la mitad) más la penalización mezclada, con α\alpha repartiendo la penalización entre los términos L1 y L2:

def elastic_net_objective(X, y, w, b, lam, alpha):
    """The full elastic-net loss: MSE half plus the blended penalty.

    (1 / 2n) * ||y - Xw - b||^2  +  lambda * [ alpha * ||w||_1
                                             + (1 - alpha) / 2 * ||w||^2 ]

    The 1/2n scaling on the error and the 1/2 on the L2 term are conveniences
    that make the coordinate-descent update fall out clean; they also match
    scikit-learn's objective so the two fits can be compared coefficient for
    coefficient. alpha = 1 zeroes the L2 term (pure lasso); alpha = 0 zeroes the
    L1 term (pure ridge).
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y
    mse_half = float(np.sum(resid ** 2)) / (2.0 * n)
    l1 = float(np.sum(np.abs(w)))
    l2 = float(np.sum(w ** 2))
    penalty = lam * (alpha * l1 + (1.0 - alpha) / 2.0 * l2)
    return mse_half + penalty

Y aquí está el ajuste. El descenso por coordenadas cicla sobre los coeficientes uno a la vez; para cada uno calcula ρj\rho_j, la correlación entre esa feature y el residual con la propia contribución de la feature sumada de vuelta, luego aplica la actualización de la sección de matemática — soft-threshold para la parte L1, dividir entre el factor de encogimiento para la parte L2. El intercepto se queda fuera de la penalización, así que centramos el target y leemos el intercepto como su media. Se detiene cuando un barrido completo apenas mueve algo:

def coordinate_descent(X, y, lam, alpha, n_iters=1000, tol=1e-9):
    """Fit elastic net by cycling over coefficients, one closed-form step each.

    Standard trick for a non-differentiable (L1) penalty: you can't take a
    single gradient step, but you CAN minimize the loss exactly along one
    coordinate at a time while the others are held fixed. For coefficient j that
    1-D problem has a closed form,

        rho_j = (1/n) * x_j . (y - yhat + x_j * w_j)     # correlation with the
                                                         # partial residual
        w_j   = soft_threshold(rho_j, lambda * alpha)    # L1: shrink and clamp
                / (1 + lambda * (1 - alpha))             # L2: shrink factor

    Read the update as two effects in sequence. The soft-threshold is the L1
    part: it can push w_j to exactly zero. The denominator is the L2 part: a
    constant shrink that pulls every coefficient toward zero without ever
    reaching it — this is what stabilizes the fit when features are correlated,
    so a group of correlated features gets kept together rather than one being
    picked arbitrarily. Because features are standardized, the curvature term
    (1/n) * sum(x_j^2) is 1, which is why it doesn't appear.

    y is centered here, so the intercept is just the target mean; we return it
    alongside the weights. We cycle until the largest coefficient change in a
    full sweep is a negligible fraction of the largest coefficient.
    """
    n, d = X.shape
    b = float(y.mean())
    yc = y - b                                 # center: intercept is out of the penalty
    w = np.zeros(d)
    denom = 1.0 + lam * (1.0 - alpha)          # the L2 shrink factor
    for _ in range(n_iters):
        max_change = 0.0
        max_coef = 0.0
        for j in range(d):
            resid = yc - X @ w                 # residual with the current w
            rho = (1.0 / n) * X[:, j] @ (resid + X[:, j] * w[j])
            w_j = soft_threshold(rho, lam * alpha) / denom
            max_change = max(max_change, abs(w_j - w[j]))
            max_coef = max(max_coef, abs(w_j))
            w[j] = w_j
        if max_coef > 0.0 and max_change / max_coef < tol:
            break
    return w, b

Ese es todo el algoritmo. Dos efectos, una línea de actualización, envueltos en un loop.

Míralo trabajar

Este es el punto del capítulo. Fijamos la fuerza de penalización en λ=1\lambda = 1 y barremos la proporción de mezcla α\alpha de 0 a 1 — de ridge puro en el extremo izquierdo a lasso puro en el derecho — reajustando el modelo completo desde cero en cada paso. Cada frame es un ajuste real; las barras son los diez coeficientes de ese ajuste, azul para positivos y naranja para negativos, con la línea punteada en cero. El caption nombra α\alpha, dice en qué régimen estás y cuenta cuántos coeficientes siguen siendo distintos de cero.

Presiona play y mira cambiar la forma. En α=0\alpha = 0, ridge, las diez barras están arriba — nada es cero, todo es pequeño y estable. Conforme α\alpha sube, la parte L1 empieza a morder: las barras se encogen, y una por una las débiles se van de golpe a exactamente cero y desaparecen. Para α=1\alpha = 1, lasso, tres de ellas ya no están y te quedan siete. Reinicia y córrelo las veces que quieras; es determinista, el mismo barrido cada vez.

Lo que hay que vigilar es s2, la feature de LDL — la que está correlacionada en 0.897 con s1. Sigue su barra conforme α\alpha sube. Bajo ridge y a lo largo de la mayor parte del rango de elastic net se mantiene arriba, compartiendo la carga con s1. En el extremo derecho, bajo lasso puro, cae a cero: lasso decidió que el par correlacionado solo necesita un miembro y soltó a s2. Esa única barra desapareciendo es el efecto de agrupamiento, o más bien su ausencia — es exactamente lo que elastic net está construido para prevenir, y es la razón por la que el medio de este barrido es donde quieres vivir cuando tus features están correlacionadas.

La implementación completa

El archivo entero, sin librería — el modelo, R², el soft-threshold, el objetivo y el descenso por coordenadas. Este es el código que de verdad corrió la animación de arriba:

"""Elastic net regression, built from scratch.

Elastic net is linear regression with two penalties stapled on: an L1 term that
drives coefficients to exactly zero (sparsity, lasso's trick) and an L2 term
that keeps them small and stable when features move together (ridge's trick). A
mix ratio alpha slides between the two — alpha = 1 is pure lasso, alpha = 0 is
pure ridge, and everything in between is elastic net.

We fit it by coordinate descent: cycle over one coefficient at a time and solve
its 1-D subproblem in closed form. That subproblem is a soft-threshold (the L1
part, which can snap a coefficient to zero) divided by a shrink factor (the L2
part, which never quite does). Pure NumPy — no ML library in this file. The
`# region:` markers are what the chapter's include directives pull in.

Inputs are assumed standardized: each feature column mean-zero and unit
variance, so the per-coordinate curvature (1/n) * sum(x_j^2) is 1 and every
coefficient is penalized on the same footing. Standardizing is not optional for
a penalized model — without it lambda punishes big-scale features less than
small-scale ones purely by accident.
"""

import numpy as np
import pandas as pd


# region: predict
def predict(X, w, b):
    """The linear model: yhat = Xw + b.

    Same one-liner as ordinary least squares. Elastic net changes how we CHOOSE
    w and b, never how we use them. X is (n, d), w is (d,), b is a scalar.
    """
    return X @ w + b
# endregion


# region: r2
def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of variance explained.

    1 minus (our squared error / the error of always guessing the mean). 1.0 is
    perfect, 0.0 is no better than the mean, negative is worse than the mean.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot
# endregion


# region: soft_threshold
def soft_threshold(z, gamma):
    """The L1 operator: shrink z toward zero by gamma, and clamp at zero.

    S(z, gamma) = sign(z) * max(|z| - gamma, 0). This is the whole reason
    elastic net (and lasso) produce exact zeros: if the unpenalized pull on a
    coefficient is weaker than gamma, the coefficient is set to precisely 0, not
    to something small. Ridge's smooth L2 penalty can't do this — it scales
    coefficients down but never crosses zero.
    """
    return np.sign(z) * max(abs(z) - gamma, 0.0)
# endregion


# region: objective
def elastic_net_objective(X, y, w, b, lam, alpha):
    """The full elastic-net loss: MSE half plus the blended penalty.

    (1 / 2n) * ||y - Xw - b||^2  +  lambda * [ alpha * ||w||_1
                                             + (1 - alpha) / 2 * ||w||^2 ]

    The 1/2n scaling on the error and the 1/2 on the L2 term are conveniences
    that make the coordinate-descent update fall out clean; they also match
    scikit-learn's objective so the two fits can be compared coefficient for
    coefficient. alpha = 1 zeroes the L2 term (pure lasso); alpha = 0 zeroes the
    L1 term (pure ridge).
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y
    mse_half = float(np.sum(resid ** 2)) / (2.0 * n)
    l1 = float(np.sum(np.abs(w)))
    l2 = float(np.sum(w ** 2))
    penalty = lam * (alpha * l1 + (1.0 - alpha) / 2.0 * l2)
    return mse_half + penalty
# endregion


# region: coordinate_descent
def coordinate_descent(X, y, lam, alpha, n_iters=1000, tol=1e-9):
    """Fit elastic net by cycling over coefficients, one closed-form step each.

    Standard trick for a non-differentiable (L1) penalty: you can't take a
    single gradient step, but you CAN minimize the loss exactly along one
    coordinate at a time while the others are held fixed. For coefficient j that
    1-D problem has a closed form,

        rho_j = (1/n) * x_j . (y - yhat + x_j * w_j)     # correlation with the
                                                         # partial residual
        w_j   = soft_threshold(rho_j, lambda * alpha)    # L1: shrink and clamp
                / (1 + lambda * (1 - alpha))             # L2: shrink factor

    Read the update as two effects in sequence. The soft-threshold is the L1
    part: it can push w_j to exactly zero. The denominator is the L2 part: a
    constant shrink that pulls every coefficient toward zero without ever
    reaching it — this is what stabilizes the fit when features are correlated,
    so a group of correlated features gets kept together rather than one being
    picked arbitrarily. Because features are standardized, the curvature term
    (1/n) * sum(x_j^2) is 1, which is why it doesn't appear.

    y is centered here, so the intercept is just the target mean; we return it
    alongside the weights. We cycle until the largest coefficient change in a
    full sweep is a negligible fraction of the largest coefficient.
    """
    n, d = X.shape
    b = float(y.mean())
    yc = y - b                                 # center: intercept is out of the penalty
    w = np.zeros(d)
    denom = 1.0 + lam * (1.0 - alpha)          # the L2 shrink factor
    for _ in range(n_iters):
        max_change = 0.0
        max_coef = 0.0
        for j in range(d):
            resid = yc - X @ w                 # residual with the current w
            rho = (1.0 / n) * X[:, j] @ (resid + X[:, j] * w[j])
            w_j = soft_threshold(rho, lam * alpha) / denom
            max_change = max(max_change, abs(w_j - w[j]))
            max_coef = max(max_coef, abs(w_j))
            w[j] = w_j
        if max_coef > 0.0 and max_change / max_coef < tol:
            break
    return w, b
# endregion


def standardize(X_train, X_other=None):
    """Z-score using the training mean/std so test data can't peek.

    Returns the standardized train matrix (and, if given, the other matrix
    scaled by the SAME statistics). Population std matches the (1/n) curvature
    the coordinate-descent update assumes.
    """
    mu = X_train.mean(axis=0)
    sd = X_train.std(axis=0)
    Xtr = (X_train - mu) / sd
    if X_other is None:
        return Xtr, mu, sd
    return Xtr, (X_other - mu) / sd, mu, sd


def load_data(path="../data/diabetes.csv"):
    """Diabetes: 442 patients, 10 standardized features, one-year progression.

    Returns (X, y, feature_names). y is a quantitative measure of disease
    progression a year after baseline. The serum-cholesterol features (s1..s6)
    are strongly correlated — s1 and s2 sit at 0.90 — which is exactly the mess
    elastic net is built for.
    """
    df = pd.read_csv(path)
    target = "target"
    features = [c for c in df.columns if c != target]
    X = df[features].to_numpy(float)
    y = df[target].to_numpy(float)
    return X, y, features

La versión de librería

Nadie escribe descenso por coordenadas a mano en producción. sklearn.linear_model.ElasticNet es el mismo modelo con el mismo objetivo, ajustado por el mismo algoritmo — el descenso por coordenadas de glmnet — solo que escrito en Cython y endurecido. Lo único que hay que hacer bien es la nomenclatura, porque los parámetros de scikit-learn chocan con la matemática. Lo que nosotros llamamos λ\lambda, la fuerza de penalización, sklearn lo llama alpha. Lo que llamamos α\alpha, la mezcla L1/L2, sklearn lo llama l1_ratio. Así que nuestro λ=1\lambda = 1, α=0.5\alpha = 0.5 es el alpha=1.0, l1_ratio=0.5 de sklearn:

def sklearn_fit(X_train, y_train, lam, alpha):
    """Fit elastic net with sklearn. `lam` is our penalty strength (sklearn's
    `alpha`); `alpha` is our L1/L2 mix (sklearn's `l1_ratio`). Returns
    (weights, intercept) so the coefficients line up one for one with ours."""
    model = ElasticNet(
        alpha=lam, l1_ratio=alpha,
        fit_intercept=True, max_iter=200000, tol=1e-12,
    )
    model.fit(X_train, y_train)
    return model.coef_, float(model.intercept_)

Evaluar es ajustar en train y reportar R² en el conjunto de prueba apartado — el número del cara a cara:

def sklearn_score(X_train, y_train, X_test, y_test, lam, alpha):
    """Fit on train, report test R^2 — the face-off number."""
    model = ElasticNet(
        alpha=lam, l1_ratio=alpha,
        fit_intercept=True, max_iter=200000, tol=1e-12,
    )
    model.fit(X_train, y_train)
    return float(model.score(X_test, y_test))

Como es el mismo objetivo resuelto por el mismo algoritmo sobre los mismos datos estandarizados, aterriza en los mismos coeficientes que nosotros. Las únicas diferencias que importan en la práctica son las que aquí no estamos ejercitando: el l1_ratio de sklearn tiene que ser estrictamente positivo (el extremo de ridge puro usa Ridge en su lugar), y por default no estandariza nada, así que tú eres responsable de escalar tus features antes de entregárselas — cosa que ya hicimos.

Desde cero contra librería

Divide los datos de diabetes 80/20 — 353 pacientes para entrenar, 89 apartados —, ajusta elastic net en λ=1\lambda = 1, α=0.5\alpha = 0.5 sobre la mitad de entrenamiento con ambas implementaciones, y evalúa sobre la mitad apartada:

Ambas barras miden lo mismo: R² de prueba de 0.3653. No parecido — idéntico. Nuestro descenso por coordenadas empata los coeficientes de sklearn hasta catorce decimales, porque hay un solo minimizador de un objetivo convexo y todo solver correcto lo encuentra. Las cifras exactas viven en results.json, regenerado cada vez que el código cambia, para que la prosa y la gráfica no puedan desviarse de lo que el código produjo. En esta división el ajuste conserva nueve de las diez features — puso una en cero — y ambas implementaciones ponen en cero la misma.

Para ver cómo se ve ese R² como predicciones, grafica el ajuste contra la verdad sobre los pacientes apartados — progresión predicha en la vertical, real en la horizontal, con la diagonal punteada marcando un ajuste perfecto. Los puntos sobre la línea están exactamente bien; la dispersión vertical fuera de ella es el error que el R² está midiendo:

La nube sigue la diagonal — el ajuste tiene la tendencia correcta — pero es una nube ancha, y está más aplanada que la línea: los valores reales bajos se sobre-predicen y los altos se sub-predicen, la regresión hacia la media que esperarías de un modelo que explica un tercio de la varianza. Un R² de 0.37 es modesto, y ese es el titular honesto: en datos así de pequeños y así de poco sobreajustados, la penalización no te compra exactitud. Lo que compra es la estructura de los coeficientes, y esa es la gráfica que vale la pena mirar. Aquí están los diez coeficientes bajo los tres regímenes — ridge, elastic net y lasso — ajustados sobre los datos estandarizados completos con el mismo λ=1\lambda = 1:

Mira s1 y s2, el par correlacionado. Ridge conserva ambos, en 0.28 y -1.40, compartiendo la señal entre ellos. Elastic net también conserva ambos, en -0.24 y -2.37 — la parte L2 manteniendo unido al grupo. Lasso vierte todo en s1 con -4.84 y corta s2 a exactamente cero: le entregaron un par correlacionado, tomó uno y soltó el otro. Ese es el efecto de agrupamiento en una sola gráfica. Ridge conserva las diez features pero no te da dispersión; lasso te da dispersión pero rompe el grupo correlacionado; elastic net se sienta entre ellos, disperso donde puede serlo y agrupado donde debe estarlo.

Conclusiones

Recurre a elastic net cuando tengas muchas features correlacionadas y quieras un modelo disperso en el que puedas confiar. Esa es la situación específica en la que gana, y es común — los conjuntos de features reales están llenos de casi-duplicados, y el hábito de lasso de quedarse con uno de ellos al azar y poner en cero el resto vuelve inestable su conjunto seleccionado de una forma que silenciosamente socava lo que sea que construyas encima. La parte L2 de elastic net amarra esas features correlacionadas entre sí para que la selección se quede quieta entre remuestreos, mientras la parte L1 igual te entrega los ceros. Es el default seguro entre los modelos lineales penalizados exactamente para los datos que rompen a lasso.

Sé honesto sobre lo que compra y lo que no. En este ajuste de diabetes el R² de prueba fue 0.37 con la penalización y sería más o menos igual sin ella — la regularización no rescató a un modelo que se ahogaba en overfitting, porque no había overfitting que rescatar. Lo que entregó fue estructura: un vector de coeficientes donde las features correlacionadas de colesterol suben y bajan juntas en lugar de que una sobreviva arbitrariamente. Cuando tu modelo tiene muchas más features de las que las muestras pueden fijar, ese mismo mecanismo empieza a pagar también en exactitud, pero la razón para usarlo es la estabilidad, y la exactitud es un bono cuando llega.

Y mantén los dos extremos en la cabeza, porque son los dos capítulos a los lados de este. Desliza la proporción de mezcla a cero y tienes ridge, puro encogimiento y nada de selección; deslízala a uno y tienes lasso, pura selección y nada de agrupamiento. Elastic net no es una tercera cosa — es la perilla entre ambos, y el ajuste correcto es una búsqueda con validación cruzada, no una adivinanza. Empiézala cerca del medio cuando tus features estén correlacionadas, empújala hacia lasso cuando quieras un modelo más ligero y puedas pagar la inestabilidad, y hacia ridge cuando no puedas pagar ninguna. La animación a la mitad de este capítulo es esa perilla; la habilidad completa es saber en qué punto de ella quieren sentarse tus datos.