Gradient boosting: XGBoost, LightGBM y CatBoost.

Boosting: árboles pequeños que corrigen a los anteriores

Un bosque (lección 5.5.2) promedia árboles independientes para bajar su varianza. El gradient boosting hace otra cosa: construye los árboles uno detrás de otro y cada uno corrige lo que los anteriores hacen mal. Con la pérdida $L(y, F)$ y un modelo inicial $F_0$ (la media de la etiqueta), en cada paso $m$:

1. se calcula el gradiente negativo de la pérdida en cada fila, $r_i = -\partial L(y_i, F_{m-1}(x_i)) / \partial F$ (con error cuadrático, el residuo $y_i - F_{m-1}(x_i)$);
2. se ajusta un árbol pequeño $T_m$ a esos $r_i$;
3. se suma con un peso pequeño, la tasa de aprendizaje $\eta$: $F_m = F_{m-1} + \eta\,T_m$, y se repite hasta $M$ árboles.

flowchart TB
  F0["F0 = media"] -->|residuos| T1[Árbol 1]
  T1 --> F1["F1 = F0 + η·T1"]
  F1 -->|residuos| T2[Árbol 2]
  T2 --> F2["F2 = F1 + η·T2"]

La formulación es de Friedman (2001). XGBoost (Chen y Guestrin, 2016) añadió una aproximación de segundo orden de la pérdida con regularización L2 en los valores de las hojas y un algoritmo escalable; LightGBM (Ke et al., 2017), el crecimiento por hojas y el muestreo de los gradientes; CatBoost (Prokhorenkova et al., 2018), árboles simétricos y un tratamiento de las categóricas pensado para evitar una fuga del objetivo. Las tres usan histogramas: discretizan cada feature en unos cientos de cubos y buscan cortes solo entre cubos, lo que hace el ajuste rápido con cientos de miles de filas.

En datos tabulares de muchos dominios, los árboles con boosting siguen ganando a las redes neuronales en los bancos de pruebas de tamaño medio (Grinsztajn, Oyallon y Varoquaux, 2022). En finanzas la pregunta relevante es otra: si ganan a un modelo lineal con buenas features. En datos reales de acciones de EE. UU., Gu, Kelly y Xiu (2020) encuentran que árboles y redes superan a las regresiones penalizadas y atribuyen la ventaja a interacciones no lineales entre predictores. Si la señal es casi lineal, como en el panel de este nivel, no hay nada que ganar y sí varianza que pagar. Esta lección lo mide.

El mismo problema para los tres

Los datos son el panel sintético del nivel (400 acciones ficticias, 21 features semanales, verdad publicada en public/curso/datos/nivel-5/README.md; nada de lo que sale aquí dice nada de mercados reales), con features y etiqueta en rangos normales de cada semana como en la lección 5.5.2. Hacen falta tres particiones con purga de una semana entre ellas, porque la etiqueta de la última semana de cada tramo acaba en el siguiente (lección 5.4.2):

• Entrenamiento 2012-2016 (259 semanas): para ajustar con distintos números de árboles.
• Validación 2017-2018 (103 semanas): para elegir el número de árboles. Nunca se usa para medir el resultado final.
• Prueba 2019-2025 (365 semanas): una sola vez, con el modelo reentrenado en 2012-2018.

import numpy as np
import pandas as pd
from scipy.stats import norm
from sklearn.linear_model import Ridge

f = pd.read_parquet("panel5_features.parquet")
e = pd.read_parquet("panel5_etiquetas.parquet")
d = (f.merge(e[["fecha", "ticker", "ret_sig_semana"]],
             on=["fecha", "ticker"])
      .dropna(subset=["ret_sig_semana"])
      .reset_index(drop=True))
FE = [c for c in f.columns[3:] if c != "revision_analistas"]

def g(s):
    """Rango normal de cada semana (lección 5.2.4); sin dato, 0."""
    por = s.groupby(d["fecha"])
    u = (por.rank() - 0.5) / por.transform("count")
    return pd.Series(norm.ppf(u), index=s.index).fillna(0.0)

def ic(score, mask):
    """IC de rangos de cada semana entre score y la etiqueta."""
    t = pd.DataFrame({"s": score, "f": sem[mask],
                      "y": d.loc[mask, "ret_sig_semana"].to_numpy()})
    t[["s", "y"]] = t.groupby("f")[["s", "y"]].rank()
    return t.groupby("f")[["s", "y"]].corr().xs("s", level=1)["y"]

X = pd.DataFrame({k: g(d[k]) for k in FE}).astype("float32")
X["sector"] = d["sector"].astype("category")
yr = g(d["ret_sig_semana"]).to_numpy()
sem = d["fecha"].to_numpy()
fe = lambda t: np.datetime64(t)
ent = sem < fe("2016-12-26")                          # 2012-2016
val = (sem >= fe("2017-01-06")) & (sem < fe("2018-12-24"))
tr_todo = sem < fe("2018-12-24")                      # 2012-2018
te = sem >= fe("2019-01-01")                          # 2019-2025
print({k: len(np.unique(sem[m])) for k, m in
       [("ent", ent), ("val", val), ("prueba", te)]})

ridge = Ridge(alpha=10.0)
ic_val_r = ic(ridge.fit(X.loc[ent, FE], yr[ent])
              .predict(X.loc[val, FE]), val).mean()
ic_r = ic(ridge.fit(X.loc[tr_todo, FE], yr[tr_todo])
          .predict(X.loc[te, FE]), te)
print(f"Ridge: IC validación {ic_val_r:.4f}  prueba {ic_r.mean():.4f} "
      f"± {ic_r.std() / np.sqrt(len(ic_r)):.4f}")

{'ent': 259, 'val': 103, 'prueba': 365}
Ridge: IC validación 0.0155  prueba 0.0180 ± 0.0045

La baseline: Ridge con α = 10 (con 145 000 filas es, en la práctica, MCO; lección 5.5.1) da un IC de 0,0155 en la validación de 2017-2018 entrenando con 2012-2016, y de 0,0180 ± 0,0045 en la prueba de 2019-2025 entrenando con 2012-2018. El sector se guarda como category de pandas: así lo reconocen XGBoost y LightGBM como categórica.

Una configuración, tres dialectos

Cada librería llama de una forma distinta a los mismos controles, y algunos no tienen equivalente exacto. Esta es la traducción de la configuración de la lección:

Tres diferencias de fondo. Crecimiento: XGBoost parte todos los nodos de un nivel antes de pasar al siguiente; LightGBM parte la hoja con mayor reducción de pérdida, esté donde esté, y por eso su tamaño se controla con num_leaves (con max_depth como freno); CatBoost usa el mismo corte en todos los nodos de un nivel (árboles simétricos u oblivious), muy regularizados y rápidos de evaluar. Con profundidad 2 los tres árboles tienen como mucho cuatro hojas, pero no son iguales: el de CatBoost usa una sola pregunta por nivel. Hoja mínima: CatBoost no la admite con árboles simétricos; su regularización viene de la simetría y de l2_leaf_reg. Categóricas: cada una las trata a su manera (sección «El sector como variable categórica»).

¿Por qué profundidad 2, hojas de 3 000 filas y tasa 0,02? Por lo que enseñó la lección 5.5.2: en este panel, los árboles útiles son diminutos y promedian muchas filas. Un booster con los valores por defecto (profundidad 6 en XGBoost y CatBoost, 31 hojas en LightGBM, tasa 0,1 a 0,3) sobreajusta en pocas iteraciones.

import catboost as cb
import lightgbm as lgb
import xgboost as xgb

N, LR = 800, 0.02
def modelos(semilla=0):
    return {
        "xgboost": xgb.XGBRegressor(
            n_estimators=N, learning_rate=LR, max_depth=2,
            min_child_weight=3000, subsample=0.5, colsample_bytree=0.8,
            reg_lambda=10, tree_method="hist", enable_categorical=True,
            random_state=semilla, n_jobs=4),
        "lightgbm": lgb.LGBMRegressor(
            n_estimators=N, learning_rate=LR, max_depth=2, num_leaves=4,
            min_child_samples=3000, subsample=0.5, subsample_freq=1,
            colsample_bytree=0.8, reg_lambda=10, random_state=semilla,
            n_jobs=4, deterministic=True, force_row_wise=True,
            verbose=-1),
        "catboost": cb.CatBoostRegressor(
            iterations=N, learning_rate=LR, depth=2, l2_leaf_reg=10,
            subsample=0.5, bootstrap_type="Bernoulli", rsm=0.8,
            random_seed=semilla, thread_count=4, verbose=0,
            allow_writing_files=False),
    }

def predice(nombre, m, Xm, k):
    """Predicción con los k primeros árboles."""
    if nombre == "xgboost":
        return m.predict(Xm, iteration_range=(0, k))
    if nombre == "lightgbm":
        return m.predict(Xm, num_iteration=k)
    return m.predict(Xm, ntree_end=k)

K = [25, 50, 100, 200, 400, 800]
ajustados, curva = {}, {}
for nombre, m in modelos().items():
    ajustados[nombre] = m.fit(X.loc[ent, FE], yr[ent])
    curva[nombre] = [ic(predice(nombre, m, X.loc[val, FE], k),
                        val).mean() for k in K]
pd.DataFrame(curva, index=K).round(4)

     xgboost  lightgbm  catboost
25    0.0098    0.0077    0.0058
50    0.0089    0.0073    0.0049
100   0.0083    0.0075    0.0069
200   0.0077    0.0070    0.0088
400   0.0055    0.0063    0.0095
800   0.0035    0.0041    0.0086

verbose=-1 en LightGBM, verbose=0 y allow_writing_files=False en CatBoost (si no, crea una carpeta catboost_info/ con los registros del entrenamiento) mantienen la salida limpia. Cada librería tiene su forma de predecir con los primeros $k$ árboles de un modelo ya entrenado (iteration_range, num_iteration, ntree_end): así se obtiene la curva entera con un solo ajuste.

La figura dibuja la tabla. XGBoost y LightGBM dan su mejor IC de validación con 25 árboles, el borde inferior de la rejilla (el óptimo podría estar aún más abajo), y empeoran con cada árbol añadido: de 0,0098 a 0,0035 y de 0,0077 a 0,0041. CatBoost, con árboles simétricos más rígidos, sube hasta 400 árboles (0,0095) y baja poco después. Los tres quedan por debajo de Ridge (0,0155) en todo el recorrido. A diferencia del bosque de la lección 5.5.2, el boosting sí sobreajusta al añadir árboles: cada uno se ajusta a los residuos de los anteriores, y con esta etiqueta los residuos son casi todo ruido.

La parada temprana mira la pérdida, no el IC

Las tres librerías tienen parada temprana: se les da un conjunto de validación, evalúan la pérdida tras cada árbol y paran cuando lleva $n$ iteraciones sin mejorar. Es cómodo, y hay que saber qué mide.

lg = lgb.LGBMRegressor(**modelos()["lightgbm"].get_params())
lg.fit(X.loc[ent, FE], yr[ent], eval_X=(X.loc[val, FE],),
       eval_y=(yr[val],),
       callbacks=[lgb.early_stopping(100, verbose=False)])
print("mejor iteración según el error cuadrático:", lg.best_iteration_)

mejor iteración según el error cuadrático: 14

Con el error cuadrático de validación, LightGBM se queda en 14 árboles. No es lo mismo que el máximo del IC: el error cuadrático castiga la escala de las predicciones, y con una señal tan débil cualquier predicción con varianza lo empeora aunque ordene algo mejor; una cartera por deciles solo usa el orden. Por eso la celda anterior calcula el IC de validación, semana a semana, para cada número de árboles. Se podría pasar una métrica propia a la parada temprana (un eval_metric que agrupe por semanas), pero la curva explícita enseña además lo plana y ruidosa que es.

Dos detalles de versión. En LightGBM 4.7, eval_set=[(X, y)] emite un aviso de obsolescencia: los conjuntos de validación van en eval_X y eval_y. Y la validación tiene que estar purgada: si las etiquetas fueran de 21 sesiones, una validación pegada al entrenamiento compartiría cuatro semanas de etiqueta con él y la parada temprana elegiría un modelo que memoriza (lecciones 5.4.1 y 5.4.2).

El número de árboles elegido y la prueba

Con el $k$ de mayor IC de validación, cada modelo se reentrena con todo 2012-2018 y se evalúa una sola vez en 2019-2025. La diferencia con Ridge se mide emparejada, semana a semana.

filas = []
for nombre in ajustados:
    k = K[int(np.argmax(curva[nombre]))]
    m = modelos()[nombre]
    m.set_params(**{"iterations" if nombre == "catboost"
                    else "n_estimators": k})
    m.fit(X.loc[tr_todo, FE], yr[tr_todo])
    s = ic(m.predict(X.loc[te, FE]), te)
    dif = s - ic_r
    filas.append({"modelo": nombre, "k": k,
                  "ic_val": max(curva[nombre]), "ic_prueba": s.mean(),
                  "se": s.std() / np.sqrt(len(s)),
                  "menos_ridge": dif.mean(),
                  "se_dif": dif.std() / np.sqrt(len(dif))})
pd.DataFrame(filas).set_index("modelo").round(4)

            k  ic_val  ic_prueba      se  menos_ridge  se_dif
modelo
xgboost    25  0.0098     0.0136  0.0046      -0.0044  0.0022
lightgbm   25  0.0077     0.0139  0.0047      -0.0041  0.0022
catboost  400  0.0095     0.0160  0.0046      -0.0020  0.0019

La figura de la cabecera resume el resultado. XGBoost da 0,0136, LightGBM 0,0139 y CatBoost 0,0160, frente a 0,0180 de Ridge y 0,0256 del oráculo en el mismo periodo (lección 5.5.1). Las diferencias emparejadas con Ridge son negativas en los tres: −0,0044 ± 0,0022 y −0,0041 ± 0,0022 (t ≈ −2) en XGBoost y LightGBM, −0,0020 ± 0,0019 en CatBoost. Es el resultado esperable en este panel: la verdad es casi lineal en los rangos, la única interacción plantada es pequeña (lección 5.5.2) y los boosters la pagan con varianza. Con los valores por defecto de las librerías el resultado sería peor.

El sector como variable categórica

El sector es la típica variable categórica de un panel de acciones. En el panel no tiene efecto plantado: el alfa verdadero no depende de él (los sectores tienen betas distintas, pero eso no se puede predecir sin saber hacia dónde irá el mercado). Lo correcto sería que no ayudara. Las tres librerías lo tratan de forma distinta:

• LightGBM ordena las categorías de cada nodo por su gradiente medio (sum_gradient / sum_hessian, según su documentación) y busca la mejor partición en dos grupos sobre ese orden. Con diez sectores es un corte muy flexible.
• XGBoost con enable_categorical=True hace algo parecido: con menos categorías que max_cat_to_onehot (4 por defecto) prueba una contra el resto; con más, reparte las categorías entre los dos hijos del nodo.
• CatBoost sustituye cada categoría por estadísticas del objetivo calculadas en orden: recorre las filas en una permutación y codifica cada una solo con las etiquetas de las anteriores. Así la codificación de una fila no incluye su propia etiqueta, una fuga del objetivo que sus autores identifican y corrigen.

COL = FE + ["sector"]
k = 200
cat = {
    "xgboost": xgb.XGBRegressor(**{**modelos()["xgboost"]
                                   .get_params(), "n_estimators": k}),
    "lightgbm": lgb.LGBMRegressor(**{**modelos()["lightgbm"]
                                     .get_params(), "n_estimators": k}),
    "catboost": cb.CatBoostRegressor(
        **{**modelos()["catboost"].get_params(), "iterations": k},
        cat_features=["sector"]),
}
filas = []
for nombre, m in cat.items():
    m.fit(X.loc[ent, COL], yr[ent])
    sin = ic(predice(nombre, ajustados[nombre], X.loc[val, FE], k),
             val).mean()
    con = ic(m.predict(X.loc[val, COL]), val).mean()
    filas.append({"modelo": nombre, "ic_val_sin_sector": sin,
                  "ic_val_con_sector": con})
pd.DataFrame(filas).set_index("modelo").round(4)

          ic_val_sin_sector  ic_val_con_sector
modelo
xgboost              0.0077             0.0043
lightgbm             0.0070             0.0039
catboost             0.0088             0.0077

Con 200 árboles, añadir el sector baja el IC de validación de XGBoost de 0,0077 a 0,0043 y el de LightGBM de 0,0070 a 0,0039: casi la mitad. CatBoost también empeora, pero menos (de 0,0088 a 0,0077). La figura añade el resultado del ejercicio 2: el mismo LightGBM con el sector en diez columnas one-hot apenas cambia (0,0068). La celda siguiente muestra por qué.

cortes = pd.Series(cat["lightgbm"].booster_.feature_importance("split"),
                   index=COL)
print(f"cortes en sector: {cortes['sector']} de {cortes.sum()}")
cortes.sort_values(ascending=False).head(6)

cortes en sector: 104 de 547
sector      104
ret_1s       77
ruido_5      35
mom_12_1     34
sue          31
ret_1m       27
dtype: int32

El sector es la feature con más cortes del modelo, 104 de 547, por delante de ret_1s (77), el efecto más fuerte del panel. Una partición de diez categorías elegida con el propio objetivo encuentra casi siempre algo que mejora el ajuste del entrenamiento: con 511 maneras de partir diez sectores en dos grupos, siempre hay una que coincide con el ruido de esas semanas. Detrás aparece ruido_5, ruido con colas gruesas, con 35 cortes. Una categórica con muchas modalidades y sin efecto es una fuente de sobreajuste; si hay que usarla, mejor con la protección de CatBoost, en one-hot o como estadística calculada fuera de la muestra de entrenamiento.

Determinismo y velocidad

El curso exige que cada salida sea idéntica en cada ejecución. Las tres librerías lo permiten en CPU con los ajustes de la tabla; la celda lo comprueba ajustando dos veces cada modelo.

for nombre in ["xgboost", "lightgbm", "catboost"]:
    a = modelos()[nombre]
    a.set_params(**{"iterations" if nombre == "catboost"
                    else "n_estimators": 100})
    p1 = a.fit(X.loc[ent, FE], yr[ent]).predict(X.loc[val, FE])
    p2 = a.fit(X.loc[ent, FE], yr[ent]).predict(X.loc[val, FE])
    print(f"{nombre:9s} idénticas: {np.array_equal(p1, p2)}")

xgboost   idénticas: True
lightgbm  idénticas: True
catboost  idénticas: True

Idénticas no significa iguales en otra máquina: con otra CPU u otro número de hilos, las sumas en coma flotante pueden cambiar en los últimos decimales y, en un corte que está en el límite, cambiar un árbol. Por eso las salidas del curso se redondean a lo que necesita la conclusión, y en producción se guarda el modelo entrenado en vez de reentrenarlo para reproducir una predicción. En GPU (device="cuda" en XGBoost, task_type="GPU" en CatBoost) hay operaciones que no garantizan el mismo orden de suma; el curso solo usa CPU.

La velocidad es fácil de medir y difícil de comparar, porque depende de la máquina, de los hilos y de la carga. El código para medirla, sin ejecutar en el curso (sus tiempos cambiarían en cada ejecución):

import time

for nombre, m in modelos().items():
    t0 = time.perf_counter()
    m.fit(X.loc[tr_todo, FE], yr[tr_todo])
    print(f"{nombre:9s} {time.perf_counter() - t0:5.1f} s")

En la máquina donde se escribió la lección (16 núcleos compartidos con otros procesos, 4 hilos por modelo), 800 árboles con 145 200 filas tardaron unos 3 s en XGBoost, 4 s en LightGBM y 4-4,5 s en CatBoost; Ridge, 0,04 s. Con deterministic=True y force_row_wise=True LightGBM renuncia a parte de su velocidad, y con árboles de profundidad 2 la ventaja del crecimiento por hojas no aparece. Las comparaciones de velocidad publicadas usan otros datos y otros ajustes: lo único seguro es que los tres son lo bastante rápidos para reentrenar cada semana y que Ridge es cien veces más rápido.

Errores típicos

• Usar los valores por defecto. Profundidad 6 o 31 hojas y tasa 0,1-0,3 están pensados para problemas con señal; con rendimientos sobreajustan en pocos árboles.
• Parar con la pérdida en vez de con la métrica del objetivo. El error cuadrático eligió 14 árboles; el IC, 25. Y en cualquier caso, la validación purgada.
• Elegir el número de árboles en la prueba. La curva de 2019-2025 no se mira para decidir; para eso está la validación. Tras elegir, se reentrena con entrenamiento y validación juntos.
• Meter categóricas sin efecto. El sector se lleva más cortes que la mejor feature y baja el IC a la mitad.
• Olvidar la baseline. Un IC de 0,0160 parece bueno hasta que se ve que Ridge, sobre las mismas semanas y en una fracción del tiempo, da 0,0180.
• Reproducibilidad a medias. Sin deterministic=True y force_row_wise=True en LightGBM, o sin allow_writing_files=False en CatBoost, el resultado o el directorio de trabajo cambian en cada ejecución.

Ejercicios

1. La comparación con Ridge se ha hecho con una sola partición. Repítela con un walk-forward anual de 2015 a 2025 (reentrenando cada enero con todo lo anterior y una semana de purga) para LightGBM con 25 árboles. ¿Se mantiene?

anio = pd.DatetimeIndex(sem).year.to_numpy()
p_lg = np.full(len(d), np.nan)
p_r = np.full(len(d), np.nan)
for a in range(2015, 2026):
    tr = sem < fe(f"{a}-01-01") - np.timedelta64(7, "D")
    t_ = anio == a
    m = modelos()["lightgbm"].set_params(n_estimators=25)
    p_lg[t_] = m.fit(X.loc[tr, FE], yr[tr]).predict(X.loc[t_, FE])
    p_r[t_] = (Ridge(alpha=10.0).fit(X.loc[tr, FE], yr[tr])
               .predict(X.loc[t_, FE]))
wf = anio >= 2015
i_lg, i_r = ic(p_lg[wf], wf), ic(p_r[wf], wf)
dif = i_lg - i_r
print(f"LightGBM {i_lg.mean():.4f}  Ridge {i_r.mean():.4f}  "
      f"diferencia {dif.mean():+.4f} ± "
      f"{dif.std() / np.sqrt(len(dif)):.4f}")

LightGBM 0.0137  Ridge 0.0178  diferencia -0.0041 ± 0.0018

2. Codifica el sector en diez columnas one-hot (pd.get_dummies) y repite la medida de LightGBM con 200 árboles. ¿Por qué el one-hot hace menos daño que la categórica nativa?

uno = pd.get_dummies(d["sector"], prefix="s", dtype="float32")
Xo = pd.concat([X[FE], uno], axis=1)
m = modelos()["lightgbm"].set_params(n_estimators=200)
m.fit(Xo[ent], yr[ent])
print(f"IC validación con sector en one-hot: "
      f"{ic(m.predict(Xo[val]), val).mean():.4f}")

IC validación con sector en one-hot: 0.0068

3. Sin código: un compañero sustituye el Ridge de producción por un LightGBM porque «en el backtest el IC sube de 0,018 a 0,031». El LightGBM tiene 31 hojas, 1 000 árboles, tasa 0,1, el sector y la industria como categóricas, y el número de árboles se eligió con la parada temprana sobre los dos últimos años del backtest. ¿Qué revisarías?