XGBoost

Qué cubre este capítulo

El capítulo anterior construyó gradient boosting: haces crecer árboles uno por uno, cada uno ajustado al gradiente negativo de la pérdida, y los sumas con un learning rate pequeño. Funciona, y durante quince años fue lo más fuerte que podías apuntarle a una tabla de números. Este capítulo es lo que pasó después. En 2016 XGBoost tomó esa misma idea aditiva y cambió dos cosas de cómo se construye cada árbol, y el resultado ganó tantas competencias de Kaggle que "nomás usa XGBoost" se volvió un chiste recurrente que además era cierto.

Los dos cambios son todo el capítulo. Primero, XGBoost usa una visión de segundo orden de la pérdida: cada muestra carga un gradiente y una hessiana, y el árbol minimiza una aproximación cuadrática de la pérdida en lugar de solo perseguir el gradiente. Eso es un paso de Newton, no un paso de gradiente. Segundo, el árbol viene regularizado desde adentro: una penalización L2 sobre los valores de las hojas y un castigo de complejidad por hoja están escritos directamente en la matemática del split, así que el árbol que sale ya viene podado por su propio objetivo. Junta ambas cosas y obtienes tres fórmulas limpias — el peso de la hoja, la ganancia del split y la actualización encogida — que construimos desde cero en NumPy y verificamos contra el paquete real de xgboost.

La pieza central es una animación: una escalera en 1-D que se va pegando a una curva una ronda de boosting a la vez, con un panel de pérdida cayendo debajo. Después corremos nuestro booster hecho a mano, el xgboost real y el gradient boosting simple de sklearn sobre los mismos datos de diabetes, y vemos cómo el modelo regularizado de segundo orden se adelanta.

Un poco de historia

Gradient boosting es de Jerome Friedman. En dos papers cerca del cambio de milenio — "Greedy Function Approximation: A Gradient Boosting Machine" (1999, publicado en 2001) — planteó el boosting como descenso de gradiente en el espacio de funciones: en cada paso, ajustas un weak learner al gradiente negativo de la pérdida y das un paso en esa dirección. Ese mismo año, Friedman, Trevor Hastie y Robert Tibshirani escribieron LogitBoost, que ya usaba un paso de Newton de segundo orden para la pérdida de clasificación. Todas las piezas estaban sobre la mesa para 2001.

XGBoost es de Tianqi Chen. Lo arrancó alrededor de 2014 como proyecto de investigación en la Universidad de Washington, bajo la Distributed Machine Learning Community, y en 2016 él y Carlos Guestrin publicaron "XGBoost: A Scalable Tree Boosting System" en KDD. La contribución del paper es en parte el algoritmo — el objetivo regularizado y la ganancia de split de segundo orden que construimos abajo — y en parte una hazaña de ingeniería: un buscador de splits consciente de la dispersión, patrones de acceso amigables al cache, cómputo out-of-core y una forma de proponer cortes candidatos sin escanear cada valor. La combinación importó. El paper reporta que de 29 soluciones ganadoras publicadas en Kaggle en 2015, 17 usaron XGBoost, y que todos los equipos del top-10 del KDD Cup 2015 lo usaron. Durante varios años, si los datos eran una tabla, el modelo ganador era este. LightGBM y CatBoost después le limaron algunas asperezas, pero son variaciones sobre el mismo tema que fijó Chen: tree boosting regularizado, de segundo orden, hecho para escalar.

La intuición

Empecemos con una imagen de lo que hace el boosting. Aquí hay un conjunto pequeño en 1-D — una entrada x, una salida y, noventa puntos, una curva suave enterrada bajo ruido. El cian son los datos de entrenamiento que el modelo ve; los diamantes ámbar son un set de validación apartado que no ve.

El boosting ajusta esta curva sumando. Arranca con una adivinanza plana — el promedio de y. Está mal en todos lados, y lo que le falta, punto por punto, es el residual. Ajusta un árbol poco profundo a esos residuales, súmale una fracción pequeña a la adivinanza, y la adivinanza se dobla un poco hacia los datos. Ahora hay un residual nuevo, más chico. Ajusta otro árbol a ese. Sigue. Cada árbol es una corrección pequeña a la suma de todos los árboles anteriores, y la predicción acumulada avanza en escalera hacia la curva.

El gradient boosting simple corta la historia ahí: el residual es el gradiente negativo de la pérdida cuadrática, así que "ajusta el residual" es "ajusta el gradiente". XGBoost hace una pregunta más filosa. No solo quiere saber en qué dirección dar el paso: quiere saber qué tan lejos, y para eso necesita la curvatura de la pérdida, la hessiana. Una muestra donde la pérdida es empinada y muy curva debería jalar la hoja hacia ella con más fuerza que una muestra sobre una parte plana y suave de la pérdida. La hessiana es lo que carga esa información, y es justo lo que el gradient boosting simple tira a la basura.

La matemática

El boosting es aditivo. Después de $t-1$ rondas el modelo predice $\hat{y}^{(t-1)}$, y la ronda $t$ agrega un árbol $f_t$. XGBoost escribe la pérdida que quiere que ese árbol minimice y luego la aproxima con una expansión de Taylor de segundo orden alrededor de la predicción actual:

Aquí $g_i = \partial_{\hat{y}} \ell(y_i, \hat{y}^{(t-1)})$ es el gradiente de la pérdida en la muestra $i$, y $h_i = \partial^2_{\hat{y}} \ell(y_i, \hat{y}^{(t-1)})$ es la hessiana. Ese $h_i$ en el término cuadrático es toda la diferencia contra el gradient boosting simple, que se queda solo con la parte lineal $g_i$. La pérdida entra únicamente a través de $g$ y $h$: cambia la pérdida y nada río abajo cambia salvo esos dos números.

El regularizador $\Omega$ es la segunda idea nueva. Un árbol con $T$ hojas y valores de hoja $w_1, \dots, w_T$ se penaliza por ambas cosas:

$\gamma$ cobra un precio fijo por hoja, y $\lambda$ es una penalización L2 que jala todo valor de hoja hacia cero. Ahora fija la forma del árbol y pregunta: ¿qué valor de hoja $w_j$ minimiza el objetivo? Cada hoja $j$ es dueña de un conjunto de muestras $I_j$; escribe $G_j = \sum_{i \in I_j} g_i$ y $H_j = \sum_{i \in I_j} h_i$ para sus sumas de gradiente y hessiana. El objetivo en $w_j$ es una parábola simple hacia arriba, y su mínimo tiene forma cerrada:

Ese es el peso de la hoja, la fórmula más importante del capítulo. Léela: la hoja apunta en contra del gradiente (el signo menos), escalada por la curvatura inversa (la hessiana en el denominador), y encogida hacia cero por $\lambda$. Quita $\lambda$ y queda el paso de Newton puro $-G/H$. Sustituye $w_j^{*}$ de regreso en el objetivo y obtienes el mejor valor que puede alcanzar un árbol de forma fija, su structure score:

Más bajo es mejor, así que cada hoja quiere su $G_j^2/(H_j+\lambda)$ lo más grande posible. De ahí sale la regla de split. Cuando partimos una hoja en un hijo izquierdo y uno derecho, el cambio en el structure score es la ganancia:

El corchete es cuánto baja el objetivo al darle a los dos hijos sus propios pesos de hoja en lugar de un peso compartido, y $\gamma$ es la cuota por la hoja extra. Un split solo ocurre cuando la estructura que compra le gana a $\gamma$. Esa resta es toda la poda de XGBoost: no hay pase de poda aparte, el objetivo rechaza los splits malos por su cuenta.

Dos pérdidas cubren la mayoría del trabajo. Para regresión con error cuadrático, $g_i = \hat{y}_i - y_i$ y $h_i = 1$. Para log-loss binaria con probabilidad $p_i$, $g_i = p_i - y_i$ y $h_i = p_i(1 - p_i)$. Eso es lo único que cambia entre regresión y clasificación. Todo lo de arriba — peso de hoja, ganancia, poda — es idéntico. Por último, el árbol de la ronda no se suma completo. El shrinkage lo escala por un learning rate $\eta$:

así cada árbol solo mueve la predicción una parte del camino y a los árboles posteriores todavía les queda trabajo. Ese es un tercer regularizador, apilado encima de $\lambda$ y $\gamma$.

En qué es bueno y en qué no

XGBoost es el modelo a vencer en datos tabulares, y se lo gana sin pedirte gran cosa. Hereda todo lo bueno de los árboles — sin escalar features, entradas numéricas y categóricas mezcladas, interacciones descubiertas gratis al anidar splits — y le suma las dos cosas que le faltan a un solo árbol: el ensamble de boosting baja el sesgo apilando cientos de correcciones, y el objetivo regularizado evita que ese ensamble se memorice el set de entrenamiento. Maneja valores faltantes de forma nativa aprendiendo una dirección por defecto en cada split, entrena rápido porque el buscador de splits fue diseñado para eso, y de fábrica, casi sin tuning, se ubica cerca de la punta en la mayoría de los leaderboards tabulares. Cuando los datos son una hoja de cálculo, esto es lo primero que agarro.

Los costos son reales, eso sí. Tiene un montón de perillas — número de árboles, profundidad, learning rate, $\lambda$, $\gamma$, subsampling, muestreo de columnas, peso mínimo por hijo — y todas interactúan, así que sacarle el último par de puntos porcentuales implica un presupuesto de tuning que un random forest no te pide. No es interpretable como lo es un árbol solo; un montón de trescientos árboles es una caja negra que sondeas con feature importances y valores SHAP, no un diagrama de flujo que le entregas a un experto del dominio. Y quiere datos tabulares específicamente: en imágenes, audio o texto crudo, donde la estructura es espacial o secuencial, una red neuronal se lo come vivo. XGBoost es un especialista, y su especialidad resulta ser la forma de datos más común en el mundo laboral.

Los datos

Dos datasets, dos trabajos. La animación conceptual usa el juguete en 1-D de la sección de intuición — noventa puntos de una curva ruidosa, divididos en sesenta y cinco de entrenamiento y veinticinco apartados — suficientemente chico para que veas la escalera aterrizar sobre puntos que alcanzas a distinguir.

Para los números honestos usamos el set de regresión de diabetes que viene con scikit-learn: 442 pacientes, diez mediciones base cada uno (edad, sexo, IMC, presión arterial y seis valores de suero sanguíneo, todos centrados en la media y escalados), y un target continuo que mide la progresión de la enfermedad un año después. Es un problema de regresión chico y ruidoso donde ningún modelo brilla espectacularmente — lo cual lo vuelve un buen lugar para ver a la regularización cambiar la respuesta de verdad. Un split 80/20 deja 353 pacientes para entrenar y 89 apartados, con semilla fija para que cada corrida sea idéntica.

Constrúyelo, una función a la vez

La construcción sigue la matemática en el orden en que se derivó. La pérdida entra solo a través de gradientes y hessianas, así que eso va primero — y son dos líneas por pérdida:

def squared_error_grad_hess(y, pred):
    """Gradient and hessian of 1/2 (pred - y)^2 with respect to pred.

    For squared error the derivatives are as simple as it gets: the gradient
    is the residual, and the hessian is a constant 1. This is the loss we use
    for the regression task in this chapter.
    """
    grad = pred - y            # d/dpred of 1/2 (pred - y)^2
    hess = np.ones_like(y)     # second derivative is 1 everywhere
    return grad, hess


def logistic_grad_hess(y, pred):
    """Gradient and hessian of binary log-loss, for classification.

    `pred` is the raw margin (log-odds); p is the sigmoid. The gradient is the
    familiar p - y, and the hessian is p(1-p). Swap this in for
    squared_error_grad_hess and everything downstream — the tree, the leaf
    weights, the boosting loop — is unchanged. That is the whole point of the
    second-order view: the loss only enters through g and h.
    """
    p = 1.0 / (1.0 + np.exp(-pred))
    grad = p - y
    hess = p * (1.0 - p)
    return grad, hess

Para error cuadrático el gradiente es el residual y la hessiana es una constante uno; para logística es p - y y p(1-p). En este capítulo hacemos regresión, pero la versión logística está aquí para dejar el punto concreto: cambia esta única función y todo el resto del archivo se convierte en un clasificador. Nada más se mueve.

Sigue el peso de la hoja — la fórmula a la que todo lo demás le sirve:

def leaf_weight(G, H, lam):
    """Optimal weight for a leaf holding gradient sum G and hessian sum H.

    Minimizing the regularized quadratic objective over the leaf's constant
    output w gives a closed form: w* = -G / (H + lambda). The lambda in the
    denominator is L2 regularization — it shrinks every leaf toward zero, and
    the more it dominates H the harder it pulls. Set lambda = 0 and you recover
    the plain Newton step -G/H.
    """
    return -G / (H + lam)

Una línea: $-G/(H+\lambda)$. El structure score que ordena los splits es la parte positiva de la misma expresión, $G^2/(H+\lambda)$:

def leaf_objective(G, H, lam):
    """The best (lowest) objective a leaf can reach, its structure score.

    Plug w* back into the leaf's quadratic and you get -1/2 * G^2 / (H + lambda).
    We work with the positive quantity G^2 / (H + lambda): the larger it is, the
    more this group of samples wants a leaf of its own. Split gain is just this
    score for the children minus the parent.
    """
    return (G * G) / (H + lam)

Con el score en la mano, la ganancia del split es el score de los hijos menos el del padre, partido a la mitad, menos $\gamma$ — la fórmula de la sección de matemática, tal cual:

def split_gain(G_L, H_L, G_R, H_R, lam, gamma):
    """Gain from splitting a node into a left and right child.

    Gain = 1/2 [ score(L) + score(R) - score(parent) ] - gamma, where
    score(G,H) = G^2 / (H + lambda). The bracket is how much lower the objective
    goes by giving the two children their own leaf weights instead of one shared
    weight. gamma is the price of the extra leaf: a split only survives if the
    structure it buys beats gamma. That single subtraction is XGBoost's built-in
    pruning — no split ever happens unless it pays for itself.
    """
    parent = leaf_objective(G_L + G_R, H_L + H_R, lam)
    children = leaf_objective(G_L, H_L, lam) + leaf_objective(G_R, H_R, lam)
    return 0.5 * (children - parent) - gamma

Ese - gamma del final hace un trabajo silencioso e importante. Un split que apenas mejora el structure score regresa ganancia negativa y se rechaza, así que el árbol deja de crecer exactamente donde la hoja extra deja de pagarse sola. Ahora la búsqueda de splits. Tiene la misma forma que la búsqueda de CART de la semana 3 — cada feature, cada umbral — pero el score es la ganancia de XGBoost, y barremos cada feature en una sola pasada ordenada, acumulando $G_L$ y $H_L$ conforme los puntos cruzan al hijo izquierdo:

def best_split(X, g, h, lam, gamma, min_child_weight):
    """Find the (feature, threshold) with the largest split gain, or None.

    This is the greedy exact split search, the same shape as the CART search
    from week 3 — loop every feature, every candidate threshold — but the score
    is XGBoost's gain built from gradient and hessian sums, not Gini. For each
    feature we sort the rows and sweep the threshold once, accumulating G_L and
    H_L as points move to the left child, so scoring every cut on a feature is
    one linear pass. A split needs positive gain (it must beat gamma) and each
    child needs at least `min_child_weight` total hessian.
    """
    n, d = X.shape
    G, H = float(g.sum()), float(h.sum())
    best = None  # (gain, feature, threshold)
    for f in range(d):
        order = np.argsort(X[:, f], kind="mergesort")
        xf = X[order, f]
        gf, hf = g[order], h[order]
        G_L = H_L = 0.0
        for i in range(n - 1):
            G_L += float(gf[i])
            H_L += float(hf[i])
            # can only cut between two different feature values
            if xf[i] == xf[i + 1]:
                continue
            G_R, H_R = G - G_L, H - H_L
            if H_L < min_child_weight or H_R < min_child_weight:
                continue
            gain = split_gain(G_L, H_L, G_R, H_R, lam, gamma)
            thr = (xf[i] + xf[i + 1]) / 2.0
            if best is None or gain > best[0]:
                best = (gain, f, float(thr))
    if best is None or best[0] <= 0.0:
        return None
    return best[1], best[2], best[0]

El barrido ordenado es el truco que abarata esto: ordenas una feature una vez, y luego mover el umbral de un valor al siguiente solo suma el gradiente y la hessiana de un punto a las sumas izquierdas acumuladas. Cada corte candidato de una feature se evalúa en una pasada lineal, no en un reconteo desde cero. La búsqueda no devuelve nada cuando ningún split tiene ganancia positiva — esa es la regla de paro podada por $\gamma$ haciendo su trabajo. Hacer crecer el árbol es esa búsqueda, aplicada de forma recursiva, con el valor de cada hoja definido por la fórmula del peso:

def build_tree(X, g, h, lam, gamma, max_depth, min_child_weight, depth=0):
    """Grow one regularized tree that maps samples to leaf weights.

    A node is a dict. A leaf carries the weight w* = -G/(H+lambda) computed from
    the gradient and hessian sums of the rows that reached it; an internal node
    carries the feature and threshold plus its two children. We stop and make a
    leaf when depth runs out or no split has positive gain — which is exactly
    the gamma-pruned criterion from best_split. The tree predicts a correction,
    not a label: its output is added to the running prediction.
    """
    G, H = float(g.sum()), float(h.sum())
    node = {"n": int(len(g)), "weight": float(leaf_weight(G, H, lam))}
    if depth >= max_depth:
        node["leaf"] = True
        return node
    split = best_split(X, g, h, lam, gamma, min_child_weight)
    if split is None:
        node["leaf"] = True
        return node
    f, t, gain = split
    left = X[:, f] <= t
    node.update({
        "leaf": False, "feature": int(f), "threshold": t, "gain": float(gain),
        "left": build_tree(X[left], g[left], h[left], lam, gamma,
                            max_depth, min_child_weight, depth + 1),
        "right": build_tree(X[~left], g[~left], h[~left], lam, gamma,
                            max_depth, min_child_weight, depth + 1),
    })
    return node

Una hoja carga $w^{*} = -G/(H+\lambda)$ calculado con las sumas de gradiente y hessiana de los renglones que llegaron a ella. Fíjate en lo que predice el árbol: no una etiqueta, sino una corrección — un número pequeño para sumarse a la predicción acumulada. Rutear una muestra hasta su hoja es el recorrido normal del árbol:

def tree_predict_one(node, x):
    """Route one sample to its leaf and return that leaf's weight."""
    while not node["leaf"]:
        node = node["left"] if x[node["feature"]] <= node["threshold"] \
            else node["right"]
    return node["weight"]


def tree_predict(tree, X):
    """The correction this one tree adds for every row of X."""
    return np.array([tree_predict_one(tree, x) for x in X])

Finalmente el loop de boosting amarra todo. Arranca cada predicción en el base score (la media de y, la mejor constante bajo error cuadrático), y en cada ronda calcula el gradiente y la hessiana en la predicción actual, haz crecer un árbol regularizado sobre ellos, y suma su salida encogida:

class XGBScratch:
    """A from-scratch XGBoost regressor: second-order boosting with shrinkage.

    Fit starts every prediction at a constant base score (the mean of y, which
    is the optimal constant under squared error). Each round it computes the
    gradient and hessian of the loss at the current prediction, grows one
    regularized tree on those, and adds a shrunken version of its output to the
    running prediction. Shrinkage (learning_rate) is the third regularizer:
    each tree only moves the prediction part of the way, so later trees still
    have work to do and no single tree dominates.
    """

    def __init__(self, n_estimators=200, learning_rate=0.3, max_depth=3,
                 reg_lambda=1.0, gamma=0.0, min_child_weight=1.0,
                 grad_hess=squared_error_grad_hess):
        self.n_estimators = n_estimators
        self.learning_rate = learning_rate
        self.max_depth = max_depth
        self.reg_lambda = reg_lambda
        self.gamma = gamma
        self.min_child_weight = min_child_weight
        self.grad_hess = grad_hess
        self.trees = []
        self.base_score = 0.0

    def fit(self, X, y, record=False):
        """Grow n_estimators trees. With record=True, keep the training
        prediction after every round (the animation and the loss curves read
        this)."""
        self.base_score = float(np.mean(y))
        pred = np.full(len(y), self.base_score, dtype=float)
        self.trees = []
        history = []
        for _ in range(self.n_estimators):
            g, h = self.grad_hess(y, pred)
            tree = build_tree(X, g, h, self.reg_lambda, self.gamma,
                              self.max_depth, self.min_child_weight)
            pred = pred + self.learning_rate * tree_predict(tree, X)
            self.trees.append(tree)
            if record:
                history.append(pred.copy())
        return history

    def staged_predict(self, X):
        """Prediction after each round — the additive model unrolled. Used to
        replay the ensemble sharpening one tree at a time."""
        pred = np.full(len(X), self.base_score, dtype=float)
        preds = []
        for tree in self.trees:
            pred = pred + self.learning_rate * tree_predict(tree, X)
            preds.append(pred.copy())
        return preds

    def predict(self, X):
        """Final prediction: base score plus every shrunken tree."""
        pred = np.full(len(X), self.base_score, dtype=float)
        for tree in self.trees:
            pred = pred + self.learning_rate * tree_predict(tree, X)
        return pred

staged_predict es la única pieza que no es estrictamente necesaria para ajustar — desenrolla el ensamble un árbol a la vez, para que podamos reproducir la predicción después de cada ronda. Eso es lo que lee la animación.

Míralo funcionar

Aquí está la recompensa. Este es el booster real del archivo de arriba, ajustado al juguete en 1-D — treinta rondas, árboles de profundidad 2, learning rate 0.3, $\lambda = 1$. El panel de arriba es la predicción: una escalera (naranja) sobre los puntos de entrenamiento (cian) y los puntos de validación apartados (diamantes ámbar). El panel de abajo traza el error cuadrático medio en ambos splits conforme se acumulan las rondas. Cada frame es una ronda real de boosting.

Dale play. La ronda 0 es el base score — una línea plana en la media, mal en todos lados, con el MSE de train y el de val ambos altos. La ronda 1 agrega un árbol de profundidad 2, y la línea plana levanta sus primeros dos escalones donde la curva se dobla más. De ahí en adelante cada ronda agrega otro árbol poco profundo que empuja la escalera más cerca, y el panel de pérdida cae rápido al principio — las correcciones grandes y fáciles llegan temprano — y luego se aplana conforme los árboles empiezan a pelearse por el ruido. Para la ronda 30 la escalera abraza la curva, el MSE de entrenamiento baja a 0.038, y el MSE de validación se asienta alrededor de 0.110. Reinicia y córrela otra vez; sale igual siempre.

Observa la brecha entre las dos líneas de pérdida. Al principio caen juntas — el modelo está aprendiendo estructura real que ayuda en ambas. Ya tarde, la línea de train sigue bajando de a poquito mientras la de val se detiene. Esa separación es el modelo empezando a ajustar ruido, y es exactamente lo que $\lambda$, $\gamma$ y el learning rate pequeño están ahí para frenar. Para ver esa palanca directamente vamos a los datos reales. Este es el RMSE de test en los pacientes de diabetes apartados, ronda por ronda, para tres valores de $\lambda$ con una profundidad fija de 4:

Esta es la gráfica estelar de la regularización. La línea naranja es $\lambda = 0$ — sin penalización L2, el paso de Newton crudo. Cae más rápido y toca su mejor RMSE de test de 53.07 en la ronda 6, luego da la vuelta y sube: de ahí en adelante cada árbol está memorizando ruido de entrenamiento y la generalización empeora. La línea cian, $\lambda = 1$, baja más lento pero llega más abajo, tocando fondo en 52.11 alrededor de la ronda 11 antes de que ella también empiece a subir. La línea morada, $\lambda = 30$, está sobrerregularizada — es tan cautelosa que nunca hace overfitting, pero tampoco baja nunca, atorada alrededor de 55.44. Ese es todo el tradeoff en una sola gráfica: muy poca regularización y el modelo hace overfitting rápido, demasiada y hace underfitting, y el ajuste útil es una penalización moderada que alcanza un mínimo más bajo y lo sostiene más tiempo. El gradient boosting simple solo tiene el learning rate y el número de árboles con qué jugar; XGBoost te entrega $\lambda$ y $\gamma$ como dos más.

La implementación completa

El booster entero, sin librería — las derivadas de la pérdida, las tres fórmulas regularizadas, la búsqueda de splits, el árbol y el loop de boosting. Este es el archivo que corrió la animación:

"""XGBoost from scratch: regularized, second-order gradient boosting.

Plain gradient boosting (week 31) fits each new tree to the negative gradient
of the loss — a first-order step. XGBoost keeps the same additive shape but
changes two things. First, it uses a second-order (Newton) view of the loss:
every sample carries a gradient g AND a hessian h, and the tree is built to
minimize a quadratic approximation of the loss, not just chase the gradient.
Second, the tree itself is regularized — an L2 penalty lambda on the leaf
weights and a per-leaf complexity penalty gamma bake straight into the split
math, so the model that comes out is already pruned by its own objective.

Everything here is pure NumPy (pandas only loads the CSV). The three formulas
that make it XGBoost rather than a plain GBM are the leaf weight
w* = -G/(H+lambda), the split gain with its gamma penalty, and the Newton step
that feeds them gradients and hessians. Each appears in the chapter one region
at a time (the `# region:` markers are what the include directives pull in).

We do regression (squared error) end to end; the logistic gradient/hessian is
included so you can see the only thing that changes for classification is the
two-line loss.
"""

import numpy as np
import pandas as pd


# region: losses
def squared_error_grad_hess(y, pred):
    """Gradient and hessian of 1/2 (pred - y)^2 with respect to pred.

    For squared error the derivatives are as simple as it gets: the gradient
    is the residual, and the hessian is a constant 1. This is the loss we use
    for the regression task in this chapter.
    """
    grad = pred - y            # d/dpred of 1/2 (pred - y)^2
    hess = np.ones_like(y)     # second derivative is 1 everywhere
    return grad, hess


def logistic_grad_hess(y, pred):
    """Gradient and hessian of binary log-loss, for classification.

    `pred` is the raw margin (log-odds); p is the sigmoid. The gradient is the
    familiar p - y, and the hessian is p(1-p). Swap this in for
    squared_error_grad_hess and everything downstream — the tree, the leaf
    weights, the boosting loop — is unchanged. That is the whole point of the
    second-order view: the loss only enters through g and h.
    """
    p = 1.0 / (1.0 + np.exp(-pred))
    grad = p - y
    hess = p * (1.0 - p)
    return grad, hess
# endregion


# region: leaf_weight
def leaf_weight(G, H, lam):
    """Optimal weight for a leaf holding gradient sum G and hessian sum H.

    Minimizing the regularized quadratic objective over the leaf's constant
    output w gives a closed form: w* = -G / (H + lambda). The lambda in the
    denominator is L2 regularization — it shrinks every leaf toward zero, and
    the more it dominates H the harder it pulls. Set lambda = 0 and you recover
    the plain Newton step -G/H.
    """
    return -G / (H + lam)
# endregion


# region: leaf_objective
def leaf_objective(G, H, lam):
    """The best (lowest) objective a leaf can reach, its structure score.

    Plug w* back into the leaf's quadratic and you get -1/2 * G^2 / (H + lambda).
    We work with the positive quantity G^2 / (H + lambda): the larger it is, the
    more this group of samples wants a leaf of its own. Split gain is just this
    score for the children minus the parent.
    """
    return (G * G) / (H + lam)
# endregion


# region: split_gain
def split_gain(G_L, H_L, G_R, H_R, lam, gamma):
    """Gain from splitting a node into a left and right child.

    Gain = 1/2 [ score(L) + score(R) - score(parent) ] - gamma, where
    score(G,H) = G^2 / (H + lambda). The bracket is how much lower the objective
    goes by giving the two children their own leaf weights instead of one shared
    weight. gamma is the price of the extra leaf: a split only survives if the
    structure it buys beats gamma. That single subtraction is XGBoost's built-in
    pruning — no split ever happens unless it pays for itself.
    """
    parent = leaf_objective(G_L + G_R, H_L + H_R, lam)
    children = leaf_objective(G_L, H_L, lam) + leaf_objective(G_R, H_R, lam)
    return 0.5 * (children - parent) - gamma
# endregion


# region: best_split
def best_split(X, g, h, lam, gamma, min_child_weight):
    """Find the (feature, threshold) with the largest split gain, or None.

    This is the greedy exact split search, the same shape as the CART search
    from week 3 — loop every feature, every candidate threshold — but the score
    is XGBoost's gain built from gradient and hessian sums, not Gini. For each
    feature we sort the rows and sweep the threshold once, accumulating G_L and
    H_L as points move to the left child, so scoring every cut on a feature is
    one linear pass. A split needs positive gain (it must beat gamma) and each
    child needs at least `min_child_weight` total hessian.
    """
    n, d = X.shape
    G, H = float(g.sum()), float(h.sum())
    best = None  # (gain, feature, threshold)
    for f in range(d):
        order = np.argsort(X[:, f], kind="mergesort")
        xf = X[order, f]
        gf, hf = g[order], h[order]
        G_L = H_L = 0.0
        for i in range(n - 1):
            G_L += float(gf[i])
            H_L += float(hf[i])
            # can only cut between two different feature values
            if xf[i] == xf[i + 1]:
                continue
            G_R, H_R = G - G_L, H - H_L
            if H_L < min_child_weight or H_R < min_child_weight:
                continue
            gain = split_gain(G_L, H_L, G_R, H_R, lam, gamma)
            thr = (xf[i] + xf[i + 1]) / 2.0
            if best is None or gain > best[0]:
                best = (gain, f, float(thr))
    if best is None or best[0] <= 0.0:
        return None
    return best[1], best[2], best[0]
# endregion


# region: build_tree
def build_tree(X, g, h, lam, gamma, max_depth, min_child_weight, depth=0):
    """Grow one regularized tree that maps samples to leaf weights.

    A node is a dict. A leaf carries the weight w* = -G/(H+lambda) computed from
    the gradient and hessian sums of the rows that reached it; an internal node
    carries the feature and threshold plus its two children. We stop and make a
    leaf when depth runs out or no split has positive gain — which is exactly
    the gamma-pruned criterion from best_split. The tree predicts a correction,
    not a label: its output is added to the running prediction.
    """
    G, H = float(g.sum()), float(h.sum())
    node = {"n": int(len(g)), "weight": float(leaf_weight(G, H, lam))}
    if depth >= max_depth:
        node["leaf"] = True
        return node
    split = best_split(X, g, h, lam, gamma, min_child_weight)
    if split is None:
        node["leaf"] = True
        return node
    f, t, gain = split
    left = X[:, f] <= t
    node.update({
        "leaf": False, "feature": int(f), "threshold": t, "gain": float(gain),
        "left": build_tree(X[left], g[left], h[left], lam, gamma,
                            max_depth, min_child_weight, depth + 1),
        "right": build_tree(X[~left], g[~left], h[~left], lam, gamma,
                            max_depth, min_child_weight, depth + 1),
    })
    return node
# endregion


# region: tree_predict
def tree_predict_one(node, x):
    """Route one sample to its leaf and return that leaf's weight."""
    while not node["leaf"]:
        node = node["left"] if x[node["feature"]] <= node["threshold"] \
            else node["right"]
    return node["weight"]


def tree_predict(tree, X):
    """The correction this one tree adds for every row of X."""
    return np.array([tree_predict_one(tree, x) for x in X])
# endregion


# region: boost
class XGBScratch:
    """A from-scratch XGBoost regressor: second-order boosting with shrinkage.

    Fit starts every prediction at a constant base score (the mean of y, which
    is the optimal constant under squared error). Each round it computes the
    gradient and hessian of the loss at the current prediction, grows one
    regularized tree on those, and adds a shrunken version of its output to the
    running prediction. Shrinkage (learning_rate) is the third regularizer:
    each tree only moves the prediction part of the way, so later trees still
    have work to do and no single tree dominates.
    """

    def __init__(self, n_estimators=200, learning_rate=0.3, max_depth=3,
                 reg_lambda=1.0, gamma=0.0, min_child_weight=1.0,
                 grad_hess=squared_error_grad_hess):
        self.n_estimators = n_estimators
        self.learning_rate = learning_rate
        self.max_depth = max_depth
        self.reg_lambda = reg_lambda
        self.gamma = gamma
        self.min_child_weight = min_child_weight
        self.grad_hess = grad_hess
        self.trees = []
        self.base_score = 0.0

    def fit(self, X, y, record=False):
        """Grow n_estimators trees. With record=True, keep the training
        prediction after every round (the animation and the loss curves read
        this)."""
        self.base_score = float(np.mean(y))
        pred = np.full(len(y), self.base_score, dtype=float)
        self.trees = []
        history = []
        for _ in range(self.n_estimators):
            g, h = self.grad_hess(y, pred)
            tree = build_tree(X, g, h, self.reg_lambda, self.gamma,
                              self.max_depth, self.min_child_weight)
            pred = pred + self.learning_rate * tree_predict(tree, X)
            self.trees.append(tree)
            if record:
                history.append(pred.copy())
        return history

    def staged_predict(self, X):
        """Prediction after each round — the additive model unrolled. Used to
        replay the ensemble sharpening one tree at a time."""
        pred = np.full(len(X), self.base_score, dtype=float)
        preds = []
        for tree in self.trees:
            pred = pred + self.learning_rate * tree_predict(tree, X)
            preds.append(pred.copy())
        return preds

    def predict(self, X):
        """Final prediction: base score plus every shrunken tree."""
        pred = np.full(len(X), self.base_score, dtype=float)
        for tree in self.trees:
            pred = pred + self.learning_rate * tree_predict(tree, X)
        return pred
# endregion


def rmse(y, pred):
    """Root mean squared error."""
    return float(np.sqrt(np.mean((y - pred) ** 2)))


def r2_score(y, pred):
    """Coefficient of determination."""
    ss_res = float(np.sum((y - pred) ** 2))
    ss_tot = float(np.sum((y - np.mean(y)) ** 2))
    return 1.0 - ss_res / ss_tot


def load_diabetes(path="../data/diabetes-regression.csv"):
    """The sklearn diabetes regression set: 442 patients, 10 features, a
    continuous disease-progression target."""
    df = pd.read_csv(path)
    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 con librería

Nadie manda a producción un booster hecho a mano, y tú tampoco deberías una vez que ya construiste uno. El paquete real de xgboost es el mismo algoritmo que acabamos de escribir — boosting de segundo orden, el peso de hoja $-G/(H+\lambda)$, la ganancia penalizada por $\gamma$ — con una década de ingeniería encima. Para que la comparación sea honesta dejamos las perillas iguales y le pedimos que use el buscador de splits exacto:

def fit_xgboost(X_train, y_train, X_test,
                n_estimators, learning_rate, max_depth,
                reg_lambda, gamma, min_child_weight, base_score):
    """The real XGBoost, configured to match our scratch booster.

    tree_method="exact" makes XGBoost do the same greedy full-scan split search
    we do (the fast default, "hist", buckets features and would diverge). We
    pass the same lambda, gamma, learning rate, depth, and base score, and turn
    off the subsampling and L1 that our version doesn't have, so any remaining
    difference is implementation detail, not a different model.
    """
    model = xgboost.XGBRegressor(
        n_estimators=n_estimators,
        learning_rate=learning_rate,
        max_depth=max_depth,
        reg_lambda=reg_lambda,
        gamma=gamma,
        min_child_weight=min_child_weight,
        reg_alpha=0.0,
        subsample=1.0,
        colsample_bytree=1.0,
        tree_method="exact",
        base_score=base_score,
        objective="reg:squarederror",
        random_state=0,
    )
    model.fit(X_train, y_train)
    return model.predict(X_test), model

El único ajuste que importa para que empaten es tree_method="exact". El default rápido de XGBoost, "hist", mete cada feature en bins y escanea los bins en lugar de los valores crudos — una aceleración enorme en datos grandes, pero encontraría splits ligeramente distintos a los de nuestro barrido completo, y los números no cuadrarían. Con "exact" hace la misma búsqueda greedy que hacemos nosotros. También apagamos las cosas que nuestra versión no tiene — muestreo de columnas y de renglones, la penalización L1 — para que cualquier diferencia que quede sea detalle de implementación, no un modelo distinto. Como contraste traemos el gradient boosting simple de sklearn, la máquina de primer orden del capítulo anterior, sin hessiana y sin regularización explícita de hojas:

def fit_sklearn_gbm(X_train, y_train, X_test,
                    n_estimators, learning_rate, max_depth):
    """sklearn's GradientBoostingRegressor — plain, first-order boosting.

    No hessian, no reg_lambda, no gamma; it fits each tree to the residual and
    scales by the learning rate. Same tree count, depth, and shrinkage as the
    others so the only difference on the face-off is the regularized
    second-order objective.
    """
    model = GradientBoostingRegressor(
        n_estimators=n_estimators,
        learning_rate=learning_rate,
        max_depth=max_depth,
        random_state=0,
    )
    model.fit(X_train, y_train)
    return model.predict(X_test), model

Desde cero contra la librería

Tres boosters, los mismos 200 árboles a profundidad 3, el mismo learning rate de 0.3, ajustados sobre los 353 pacientes de entrenamiento y evaluados sobre los 89 apartados. Nuestro booster de segundo orden hecho desde cero, el xgboost real y el gradient boosting simple de sklearn. RMSE de test, más bajo es mejor:

Aquí hay dos cosas que leer. Primero, nuestro booster y el xgboost real quedan uno encima del otro — 59.50 de RMSE de test para el nuestro, 60.01 para el paquete, una diferencia de menos del uno por ciento, y sus predicciones sobre el set apartado correlacionan en 0.9978. Ese es el resultado que quieres de una construcción desde cero: el mismo modelo, confirmado por la implementación de referencia. La brecha pequeña que queda son las partes que dejamos fuera a propósito — el manejo propio del base score por defecto del paquete, el orden de suma en punto flotante, el desempate en la búsqueda de splits — no un desacuerdo sobre el algoritmo. Las cifras exactas viven en results.json, regeneradas cada que cambia el código.

Segundo, y este es el punto del capítulo: ambos boosters de segundo orden le ganan al gradient boosting simple. El GBM de sklearn, con el mismo número de árboles, profundidad y learning rate pero sin hessiana y sin regularización de hojas, llega a 64.29 de RMSE — una R cuadrada de 0.168 contra 0.288 del nuestro y 0.275 de xgboost. En este set chico y ruidoso el objetivo de Newton regularizado vale como doce puntos de RMSE sobre la máquina de primer orden, gratis, sin tuning extra. Esos son los doce puntos que ganaron todas esas competencias.

Puntos clave

XGBoost es gradient boosting con las dos mejoras que convirtieron a un modelo fuerte en el modelo por defecto. El objetivo de segundo orden usa la curvatura de la pérdida, no solo su pendiente, así que cada árbol da un paso de Newton en lugar de un paso de gradiente y cae más cerca al primer intento. El score regularizado del árbol hornea una penalización L2 de hoja y un cargo de complejidad por hoja directo en la matemática del split, así que el ensamble se poda solo mientras crece en lugar de memorizarse el set de entrenamiento. Agrega el shrinkage y tienes tres frenos independientes contra el overfitting, y los números de diabetes los muestran ganándose el sueldo — doce puntos de RMSE sobre el GBM simple de la misma forma.

Agárralo cuando los datos sean una tabla y quieras el mejor número con la menor ceremonia. Maneja tipos mezclados y valores faltantes, no necesita escalado, entrena rápido y se ubica cerca de la punta de los leaderboards tabulares antes de que hayas tuneado nada. Las perillas que de verdad mueven la aguja, en el orden en que las toco: el learning rate y el número de árboles juntos (tasa más baja, más árboles, hasta que la curva de validación deje de mejorar), luego max_depth para cuánta interacción puede capturar cada árbol, y luego $\lambda$ y $\gamma$ para meterle riendas al overfitting una vez fijada la profundidad. La gráfica de regularización es el modelo mental — hay un ajuste moderado que alcanza un mínimo más bajo y lo sostiene, y tu chamba es encontrarlo con un set de validación, no subir todas las penalizaciones hasta el cielo.

¿Cuándo es exagerado? Cuando un random forest ya hace el trabajo — si no necesitas ese último par de puntos porcentuales y prefieres no andar cuidando un calendario de learning rate, los árboles con bagging te dan casi todo el accuracy con casi nada de tuning, y son mucho más difíciles de sobreajustar. Y un solo árbol profundo sigue siendo la respuesta correcta cuando necesitas algo que un humano pueda leer y discutir, cosa que un ensamble con boosting nunca será. XGBoost es a lo que gradúas cuando el accuracy tabular es lo que más importa y estás dispuesto a gastar un presupuesto de tuning para conseguirlo. En el tipo de datos con el que la mayoría de nosotros trabaja de verdad, eso pasa más seguido que no.