Capítulo 23 de 37 · intermedio

Gradient boosting

Qué cubre este capítulo

El random forest construía árboles en paralelo y los promediaba: cada árbol una conjetura independiente, y la multitud más sabia que cualquiera de sus miembros. El boosting hace lo contrario. Construye árboles en secuencia, y cada árbol nuevo existe por una sola razón: corregir los errores que los árboles anteriores siguen cometiendo. AdaBoost, el capítulo pasado, lo lograba re-ponderando los ejemplos que el ensamble seguía fallando. El gradient boosting hace algo más limpio y más general: mira los errores como un gradiente y ajusta el siguiente árbol para que apunte directo cuesta abajo.

Este es el algoritmo que calladito se lleva las competencias de datos tabulares. Lo vamos a construir desde cero en NumPy para el caso más limpio —regresión con squared-error loss—, donde el método entero se reduce a algo casi vergonzosamente simple: adivina la media, mira lo que sobra, ajusta un árbol poco profundo a esa sobra, súmale una rebanada encogida y repite. Lo que sobra tiene nombre, los residuales, y la única idea que convierte esto en un método de gradiente es que para el squared error el residual es el gradiente negativo de la loss. Ajusta el residual y ya hiciste gradient descent, solo que lo que vas moviendo no es un vector de pesos: es la propia función de predicción.

La pieza central es una animación: una ronda de boosting por frame sobre una curva en 1-D, con la predicción acumulada trepando desde una línea plana hasta convertirse en una escalera que abraza una onda senoidal mientras los residuales colapsan hacia cero por debajo. Después le damos el mismo trabajo al GradientBoostingRegressor de scikit-learn y verificamos que los números cuadren, y vemos cómo la única perilla que define el método —el learning rate— cambia velocidad por overfitting.

Un poco de historia

El boosting empezó como una pregunta teórica. En 1990 Robert Schapire demostró que un aprendiz "débil" —uno que apenas tiene que ganarle por un pelo a adivinar al azar— podía convertirse mediante boosting en uno "fuerte" y arbitrariamente preciso, un resultado sorprendente sobre lo que era posible en principio. Yoav Freund y Schapire lo volvieron un algoritmo práctico con AdaBoost en 1995, y durante unos años AdaBoost fue medio un misterio: funcionaba mucho mejor de lo que la teoría decía que debería, y nadie estaba muy seguro de por qué.

La aclaración vino desde la estadística. En 1998 Leo Breiman notó que AdaBoost estaba haciendo una especie de gradient descent, y en 2000 Jerome Friedman, Trevor Hastie y Robert Tibshirani escribieron el paper que replanteó todo el asunto: AdaBoost está ajustando un modelo aditivo mediante una optimización por etapas de una loss exponencial. Una vez que lo ves así, la loss específica deja de ser algo especial. Friedman llevó la idea hasta sus últimas consecuencias en "Greedy Function Approximation: A Gradient Boosting Machine" —un reporte técnico de 1999 publicado en los Annals of Statistics en 2001— y nos dio la receta general: escoge cualquier loss diferenciable, calcula su gradiente negativo en las predicciones actuales, ajusta un árbol de regresión a ese gradiente y da un paso. AdaBoost sale como el caso especial de una loss en particular. El squared error, el caso que construimos aquí, sale como el más simple, donde el gradiente es nada más el residual. Ese único paper es el ancestro directo de todas las librerías de gradient boosting —XGBoost, LightGBM, CatBoost— que hoy están en la cima del leaderboard de datos tabulares.

La intuición

Aquí está el truco, en términos llanos, antes de cualquier cálculo. Supón que tu modelo actual predice algún número para cada punto, y le atina mal. Los errores —por cuánto fallaste cada punto— son en sí mismos un dataset: una entrada x y un "target" igual a la falla. Entonces entrena un modelo que prediga las fallas. Si logra predecir aunque sea un poquito del patrón en tus errores, puedes sumarle sus predicciones a tu modelo y encoger los errores. Ahora tienes un modelo ligeramente mejor, con errores ligeramente más chicos, y lo vuelves a hacer. Cada ronda persigue lo que sobró.

Esa sobra es el residual, y - F(x): la verdad menos lo que el modelo dice actualmente. El primer modelo es el más flojo posible: predecir la media del target para todos, que es la mejor conjetura constante que existe. Sus residuales son grandes. Ajusta un árbol poco profundo a esos residuales, súmale una fracción a la predicción, y los residuales se encogen. Ajusta otro árbol a los nuevos residuales, ya más chicos. El modelo es la suma acumulada de esa constante inicial y de cada árbol encogido.

¿Por qué árboles poco profundos, y por qué una fracción? Porque cada árbol debe ser una corrección chiquita, no una respuesta completa. Un árbol profundo ajustado a los residuales de un jalón los borraría del training set y memorizaría el ruido junto con la señal. Un árbol chaparrito —de dos o tres niveles— solo puede capturar una pieza burda del patrón que sobró, que es justo lo que quieres cuando vas a apilar cincuenta o cien de ellos. La fracción es el learning rate, y es la misma idea que el tamaño de paso en gradient descent: da pasos chicos y te acercas con cuidado; da pasos grandes y te pasas de rosca. El método entero es gradient descent, hecho en el espacio de funciones de predicción en lugar del espacio de pesos.

Las matemáticas

El boosting construye un modelo aditivo: una predicción acumulada que arranca en una constante y a la que se le suma un árbol a la vez. Después de $M$ rondas queda

F_M(x) = F_0 + \nu \sum_{m=1}^{M} h_m(x)

donde $F_0$ es la constante inicial, cada $h_m$ es un árbol de regresión poco profundo y $\nu$ es el learning rate: un número positivo chico, típicamente 0.1, que encoge la contribución de cada árbol. El símbolo $x$ es una fila de entrada y $y$ su target; hay $n$ filas de entrenamiento.

Ajustamos el modelo un término a la vez para bajar una loss. Para regresión la loss es el squared error, escrita aquí con un un medio que deja limpia la derivada:

L(y, F) = \tfrac{1}{2}\,(y - F)^2

¿Dónde empezamos? Con la constante $F_0$ que por sí sola minimiza la loss total. Para el squared error esa constante es la media del target:

F_0 = \arg\min_{c} \sum_{i=1}^{n} L(y_i, c) = \frac{1}{n}\sum_{i=1}^{n} y_i

Ahora el paso clave. Queremos que cada árbol nuevo mueva la predicción en la dirección que reduce la loss más rápido: el gradiente negativo de la loss respecto a la predicción actual, evaluado en cada punto de entrenamiento. Deriva la loss de squared error respecto a $F$ :

\frac{\partial L(y, F)}{\partial F} = -(y - F)

Cámbiale el signo y el gradiente negativo en la fila $i$ , evaluado en el modelo actual $F_{m-1}$ , es

r_i = -\left.\frac{\partial L(y_i, F)}{\partial F}\right|_{F = F_{m-1}(x_i)} = y_i - F_{m-1}(x_i)

Eso es el residual, exactamente. Esta es toda la razón por la que el squared error es el caso limpio para aprender: el gradiente negativo que se supone debes perseguir es nada más "la verdad menos la conjetura actual", sin necesidad de cálculo para computarlo en la práctica. Ajustamos el árbol nuevo $h_m$ a esos residuales por mínimos cuadrados ordinarios —el objetivo nativo del árbol de regresión— y luego damos un paso encogido:

F_m(x) = F_{m-1}(x) + \nu\, h_m(x)

Repite durante $M$ rondas. Cambia el squared error por otra loss diferenciable —error absoluto, Huber, la loss logística para clasificación— y solo cambia una línea: la fórmula de $r_i$ . Todo lo demás, el ajuste del árbol y el paso encogido, se queda igual. Esa generalidad es el punto de la palabra "gradient" en el nombre.

En qué es bueno y en qué no

El gradient boosting es, para mi gusto, lo más fuerte que le puedes aventar a un dataset tabular sin usar una red neuronal, y normalmente también con una. Hereda los buenos modales de los árboles —nada de escalar features, entradas numéricas y categóricas mezcladas, fronteras no lineales e interacciones que encuentra solo— y después hace algo que el random forest no puede: como cada árbol corrige los errores del anterior en vez de solo promediar una conjetura independiente, baja el sesgo ronda tras ronda y alcanza una accuracy que un bosque de los mismos árboles ni de broma toca. Dale una loss diferenciable y la optimiza directo, así que puedes ajustarlo a squared error, error absoluto, cuantiles o un objetivo de clasificación cambiando una fórmula. Cuando hay un premio de verdad en un problema tabular, el modelo ganador casi siempre es un ensamble con gradient boosting.

Los costos son reales y son la imagen espejo de los del bosque. Es secuencial, así que no paraleliza como lo hace el bagging: el árbol m necesita las predicciones del árbol m-1 para calcular sus residuales. Tiene más perillas y además interactúan: el learning rate y el número de árboles se compensan entre sí, y si les atinas mal o haces underfitting o memorizas. Y a diferencia de un bosque, más árboles no siempre es más seguro: pasado cierto punto el ensamble empieza a ajustar ruido y el error de prueba da la vuelta y sube, cosa que vamos a ver pasar. Es un modelo que se tunea, no uno que ajustas y olvidas. La recompensa vale la niñera más veces que no, pero sí te toca hacer de niñera.

Los datos

Dos datasets, dos trabajos. La animación conceptual usa una curva en 1-D de la que puedo dibujar cada punto: 120 puntos con x repartida sobre [0, 10] y el target una onda senoidal sobre una pendiente suave hacia arriba, y = sin(x) + 0.25x, más ruido gaussiano. Una entrada, una salida, una forma que ningún árbol poco profundo puede capturar solo: perfecta para ver cómo una escalera de árboles se va armando hacia ella.

Para los números honestos de accuracy nos cambiamos a un set de regresión real: los datos de diabetes de scikit-learn, 442 pacientes con diez mediciones fisiológicas estandarizadas cada uno, y un target que es una medida cuantitativa de la progresión de la enfermedad un año después de la línea base. Es una regresión chica, ruidosa y genuinamente difícil —nadie saca un R² espectacular ahí— lo que la vuelve el lugar correcto para ver qué te compra realmente el boosting sobre un solo árbol y dónde empieza a doler el learning rate. El mismo estilo de split que en los capítulos anteriores: 309 pacientes para entrenar, 133 apartados para probar.

Constrúyelo, una función a la vez

Aquí se apilan dos piezas. Primero el aprendiz débil —un árbol de regresión poco profundo— y luego el ciclo de boosting que hace crecer una secuencia de ellos. El árbol es el CART de la semana 3 con una sustitución: parte para reducir el squared error en lugar de la impureza de Gini, y una hoja predice un número en vez de una clase.

La pureza de un nodo de regresión es simplemente qué tan disperso está el target adentro, que es la suma de errores cuadrados alrededor de la propia media del nodo:

def sse(y):
    """Sum of squared errors of a node around its own mean.

    This is the regression analogue of Gini impurity: how spread out the target
    is inside a node. A node holding identical values has sse 0; the more the
    values scatter, the larger it grows. A split tries to drive the total sse of
    its two children below the sse of the parent — that drop is the split's gain.
    """
    if len(y) == 0:
        return 0.0
    return float(np.sum((y - y.mean()) ** 2))

Este es el análogo exacto de Gini. Un nodo lleno de valores idénticos tiene sse cero; entre más se dispersen los valores, más grande se pone, y un split se gana su lugar dejando a sus dos hijos con menos sse total del que tenía el padre. Una hoja, cuando el particionado se detiene, predice el único número que minimiza su propio sse: la media.

def leaf_value(y):
    """The constant a leaf predicts: the mean of its rows.

    For squared-error loss the mean is the value that minimizes the leaf's own
    error, which is exactly why the tree splits to reduce sse and then predicts
    the average — the two agree on the same objective.
    """
    return float(y.mean())

La búsqueda del mejor split es la de la semana 3, línea por línea, con sse en lugar del score de Gini. Recorre cada feature, cada umbral candidato, y quédate con el split que más reduzca el error cuadrado total de las dos mitades:

def best_split(X, y):
    """Search every feature and threshold for the split that most reduces sse.

    For each feature we take the midpoints between consecutive sorted unique
    values as candidate thresholds. A split sends rows with feature <= t left and
    the rest right; its quality is the parent sse minus the summed sse of the two
    children (the variance reduction). We return the (feature, threshold, gain)
    with the largest reduction, or None if no split helps. This is week 3's
    best-split search with Gini swapped for sse — the shape is identical.
    """
    n, d = X.shape
    parent = sse(y)
    best_gain, best_f, best_t = 0.0, -1, 0.0
    for f in range(d):
        values = np.unique(X[:, f])
        thresholds = (values[:-1] + values[1:]) / 2.0
        for t in thresholds:
            left = X[:, f] <= t
            n_left = int(left.sum())
            if n_left == 0 or n_left == n:
                continue
            child = sse(y[left]) + sse(y[~left])
            gain = parent - child
            if gain > best_gain:
                best_gain, best_f, best_t = gain, f, float(t)
    if best_f < 0:
        return None
    return best_f, best_t, best_gain

Si leíste el capítulo de árboles de decisión esto te va a resultar familiar de carácter: lo único que cambia es qué significa "pureza". Hacer crecer el árbol es la misma recursión, también, cortada temprano por un tope de profundidad bajo porque en boosting queremos un aprendiz débil, no uno perfecto:

def build_tree(X, y, max_depth, depth=0):
    """Recursively grow a shallow CART regression tree.

    A node is a dict. A leaf carries the mean it predicts; an internal node
    carries the feature index and threshold to split on plus its two children.
    Recursion stops when the depth cap is hit, the node holds one row, or no
    split reduces sse — at which point the node becomes a leaf. Boosting keeps
    these trees deliberately shallow (max_depth 2-3) so each one is a weak
    learner that nudges the fit rather than memorizing it.
    """
    node = {"n": int(len(y)), "value": leaf_value(y)}
    if depth >= max_depth or len(y) <= 1:
        node["leaf"] = True
        return node
    split = best_split(X, y)
    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": gain,
        "left": build_tree(X[left], y[left], max_depth, depth + 1),
        "right": build_tree(X[~left], y[~left], max_depth, depth + 1),
    })
    return node

Un max_depth de dos o tres es toda la disciplina aquí. Cada árbol debe capturar una rebanada burda del patrón sobrante y nada más. Predecir es caminar de la raíz a una hoja y regresar la media de esa hoja:

def tree_predict_one(node, x):
    """Walk one sample from the root to a leaf, following each split."""
    while not node["leaf"]:
        node = node["left"] if x[node["feature"]] <= node["threshold"] \
            else node["right"]
    return node["value"]


def tree_predict(tree, X):
    """Predict every row of X by walking the tree from the root."""
    return np.array([tree_predict_one(tree, x) for x in X])

Ese es el aprendiz débil. Ahora la parte que lo vuelve boosting. Arranca la predicción en la media, y luego ronda tras ronda: calcula los residuales, ajusta un árbol a ellos, y súmale un paso encogido de ese árbol a la predicción acumulada:

def fit_gradient_boost(X, y, n_trees, learning_rate, max_depth, record=False):
    """Fit a gradient-boosting regressor by the residual-fitting loop.

    Start the prediction at F0, the mean of y (the constant that minimizes
    squared error). Then, round after round: compute the residual r = y - F,
    which for squared-error loss is exactly the negative gradient of the loss;
    fit a shallow tree to that residual; and take a shrunken step by adding
    learning_rate * tree(X) to the running prediction F. The model is the
    initial constant plus the list of trees. With `record` on we log the state
    of F after every round so the chapter can animate the fit building up.
    """
    F0 = float(y.mean())
    F = np.full(len(y), F0)
    trees, history = [], []
    if record:
        history.append({"round": 0, "F": F.copy(),
                        "residual": (y - F).copy(), "mse": float(np.mean((y - F) ** 2))})
    for m in range(1, n_trees + 1):
        residual = y - F                      # negative gradient of 1/2 (y-F)^2
        tree = build_tree(X, residual, max_depth)
        F = F + learning_rate * tree_predict(tree, X)
        trees.append(tree)
        if record:
            history.append({"round": m, "F": F.copy(),
                            "residual": (y - F).copy(),
                            "mse": float(np.mean((y - F) ** 2))})
    model = {"F0": F0, "trees": trees, "learning_rate": learning_rate}
    return (model, history) if record else model

Lee las tres líneas del ciclo contra las matemáticas. residual = y - F es el gradiente negativo de la loss de squared error, justo lo que acabamos de derivar. build_tree sobre ese residual es el ajuste por mínimos cuadrados de $h_m$ . F = F + learning_rate * tree_predict(...) es el paso encogido $F_m = F_{m-1} + \nu h_m$ . Tres líneas, y son el algoritmo completo; todo lo que está arriba de ellas es el árbol al que llaman. La bandera record nada más registra el estado de la predicción después de cada ronda para que la animación tenga algo que reproducir.

Predecir con el modelo terminado es la fórmula aditiva escrita tal cual: la constante inicial más cada árbol encogido:

def gb_predict(model, X):
    """Predict with the boosted model: the initial constant plus every shrunken
    tree. F(x) = F0 + nu * sum_m tree_m(x)."""
    F = np.full(X.shape[0], model["F0"])
    for tree in model["trees"]:
        F = F + model["learning_rate"] * tree_predict(tree, X)
    return F

Y dos métricas para calificarlo, el mean squared error y el $R^2$ que te dice qué tanto mejor que adivinar la media lo estás haciendo:

def mse(y_true, y_pred):
    """Mean squared error: the average squared gap between guess and truth."""
    return float(np.mean((y_true - y_pred) ** 2))


def r2_score(y_true, y_pred):
    """Coefficient of determination: the fraction of variance the model explains.

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

Míralo funcionar

Aquí está la recompensa, y es la imagen más clara del boosting que sé dibujar. Este es el fit_gradient_boost real corriendo sobre la curva en 1-D, una ronda por frame, con un learning rate deliberadamente chico de 0.1 para que la escalera se arme lo bastante lento como para verla. El panel de arriba muestra la predicción actual F(x) —una función escalonada naranja, porque una suma de árboles poco profundos es constante por tramos— trepando hacia la curva verdadera (la línea gris punteada) a través de los datos en cian. El panel de abajo son los residuales, y - F(x), la sobra que cada ronda intenta matar. El pie de figura nombra la ronda, el learning rate y el MSE de entrenamiento.

Dale play. La ronda cero es el modelo flojo: F(x) es una línea plana en la media, 1.35, y los residuales en el panel de abajo son una nube ancha —el modelo todavía no sabe nada, así que cada punto está "mal" por toda su distancia al promedio. Luego empiezan a aterrizar los árboles. Los primeros tallan los escalones más grandes y más burdos —agarran la pendiente general y la joroba más grande de la senoidal— y puedes ver la nube de residuales contraerse más fuerte en esas primeras rondas. Conforme se acumulan las rondas, los escalones se vuelven más finos y más numerosos, la escalera se dobla alrededor de la onda, y los residuales se aprietan hacia la línea del cero. El MSE de entrenamiento cae de 1.007 en la ronda cero a 0.069 para la ronda cincuenta, de forma monótona, cada ronda más chico que la anterior.

Fíjate en dos cosas en particular. Primero, el árbol de ningún frame hace mucho por sí solo: cada uno empuja tantito. Eso es el shrinkage: con un learning rate de 0.1 cada árbol aporta un décimo de lo que encontró, así que el ajuste se acerca con cuidado en vez de caer de golpe en su lugar. Segundo, mira dónde va como a la mitad de la corrida, cerca de la ronda veinticinco. La escalera ya sigue bien la curva y los residuales son más o menos del tamaño del ruido que le metí a los datos: eso es tan bueno como cualquier modelo puede hacerlo honestamente aquí. Las rondas siguientes son refinamiento, y hacia el final empiezan a perseguir ondulaciones que son ruido, no señal. Ese overfitting que se va colando es invisible en este juguete limpio en 1-D, pero con datos reales es exactamente lo que voltea el error de prueba, que es lo siguiente que vamos a ver.

El learning rate y el número de árboles son las dos perillas que definen el método, y jalan una contra la otra. Esto es lo que hacen con los datos reales de diabetes: MSE de prueba conforme se agregan árboles, para tres learning rates.

Esta es la imagen para quedarse un rato. El learning rate grande, 0.5 (ámbar), cae más rápido —llega a su mejor MSE de prueba de alrededor de 3244 con apenas 18 árboles— y luego voltea y sube en picada, sobreajustando durísimo, hasta pasar de 5300 con 300 árboles, que no es mejor que predecir la media. Corrió a toda velocidad hacia una respuesta decente y luego se cayó por un barranco. El learning rate chico, 0.05 (cian), es la tortuga: desciende despacio, tarda unos 105 árboles en llegar a su mejor valor de aproximadamente 3107, y sobreajusta tan suavemente que incluso con 300 árboles apenas se regresó a 3241. El rate de en medio, 0.1 (morado), toca fondo alrededor de 3105 cerca de los 40 árboles y luego sobreajusta a ritmo moderado. La lección que todo practicante de gradient boosting aprende: un learning rate más chico llega a un piso más bajo pero necesita más árboles para llegar ahí, y perdona más si te pasas con la cantidad de árboles. En la práctica pones el learning rate chico, la cantidad de árboles generosa, y dejas que un set de validación o el early stopping te digan dónde parar. Nunca dejas el learning rate en 0.5.

La implementación completa

Todo el asunto, sin librerías, de arriba abajo: el árbol de regresión poco profundo y el ciclo de boosting que lo apila. Este es el archivo que corrió la animación:

"""Gradient boosting for regression with squared-error loss, from scratch.

The idea in one line: start from a constant prediction (the mean of the target),
then repeatedly fit a shallow regression tree to what's left over — the residuals
— and add a shrunken step of that tree to the running prediction. The residuals
are the negative gradient of squared-error loss, so "fit the residual" and
"take a step downhill on the loss" are the same move. That's the whole method,
and it's why it generalizes: swap in a different loss and you fit its gradient
instead.

Two pieces live in this file. First a shallow CART regression tree — the weak
learner — that splits to reduce squared error instead of Gini impurity (a leaf
predicts the mean of its rows). Then the boosting loop that grows a sequence of
those trees on residuals and sums their shrunken predictions.

Pure NumPy — no ML library anywhere in this file (pandas only loads the CSV).
Every function below appears in the chapter one step at a time; the `# region:`
markers are what the book's include directives pull in.
"""

import numpy as np
import pandas as pd


# region: sse
def sse(y):
    """Sum of squared errors of a node around its own mean.

    This is the regression analogue of Gini impurity: how spread out the target
    is inside a node. A node holding identical values has sse 0; the more the
    values scatter, the larger it grows. A split tries to drive the total sse of
    its two children below the sse of the parent — that drop is the split's gain.
    """
    if len(y) == 0:
        return 0.0
    return float(np.sum((y - y.mean()) ** 2))
# endregion


# region: leaf_value
def leaf_value(y):
    """The constant a leaf predicts: the mean of its rows.

    For squared-error loss the mean is the value that minimizes the leaf's own
    error, which is exactly why the tree splits to reduce sse and then predicts
    the average — the two agree on the same objective.
    """
    return float(y.mean())
# endregion


# region: best_split
def best_split(X, y):
    """Search every feature and threshold for the split that most reduces sse.

    For each feature we take the midpoints between consecutive sorted unique
    values as candidate thresholds. A split sends rows with feature <= t left and
    the rest right; its quality is the parent sse minus the summed sse of the two
    children (the variance reduction). We return the (feature, threshold, gain)
    with the largest reduction, or None if no split helps. This is week 3's
    best-split search with Gini swapped for sse — the shape is identical.
    """
    n, d = X.shape
    parent = sse(y)
    best_gain, best_f, best_t = 0.0, -1, 0.0
    for f in range(d):
        values = np.unique(X[:, f])
        thresholds = (values[:-1] + values[1:]) / 2.0
        for t in thresholds:
            left = X[:, f] <= t
            n_left = int(left.sum())
            if n_left == 0 or n_left == n:
                continue
            child = sse(y[left]) + sse(y[~left])
            gain = parent - child
            if gain > best_gain:
                best_gain, best_f, best_t = gain, f, float(t)
    if best_f < 0:
        return None
    return best_f, best_t, best_gain
# endregion


# region: build_tree
def build_tree(X, y, max_depth, depth=0):
    """Recursively grow a shallow CART regression tree.

    A node is a dict. A leaf carries the mean it predicts; an internal node
    carries the feature index and threshold to split on plus its two children.
    Recursion stops when the depth cap is hit, the node holds one row, or no
    split reduces sse — at which point the node becomes a leaf. Boosting keeps
    these trees deliberately shallow (max_depth 2-3) so each one is a weak
    learner that nudges the fit rather than memorizing it.
    """
    node = {"n": int(len(y)), "value": leaf_value(y)}
    if depth >= max_depth or len(y) <= 1:
        node["leaf"] = True
        return node
    split = best_split(X, y)
    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": gain,
        "left": build_tree(X[left], y[left], max_depth, depth + 1),
        "right": build_tree(X[~left], y[~left], max_depth, depth + 1),
    })
    return node
# endregion


# region: tree_predict
def tree_predict_one(node, x):
    """Walk one sample from the root to a leaf, following each split."""
    while not node["leaf"]:
        node = node["left"] if x[node["feature"]] <= node["threshold"] \
            else node["right"]
    return node["value"]


def tree_predict(tree, X):
    """Predict every row of X by walking the tree from the root."""
    return np.array([tree_predict_one(tree, x) for x in X])
# endregion


# region: fit
def fit_gradient_boost(X, y, n_trees, learning_rate, max_depth, record=False):
    """Fit a gradient-boosting regressor by the residual-fitting loop.

    Start the prediction at F0, the mean of y (the constant that minimizes
    squared error). Then, round after round: compute the residual r = y - F,
    which for squared-error loss is exactly the negative gradient of the loss;
    fit a shallow tree to that residual; and take a shrunken step by adding
    learning_rate * tree(X) to the running prediction F. The model is the
    initial constant plus the list of trees. With `record` on we log the state
    of F after every round so the chapter can animate the fit building up.
    """
    F0 = float(y.mean())
    F = np.full(len(y), F0)
    trees, history = [], []
    if record:
        history.append({"round": 0, "F": F.copy(),
                        "residual": (y - F).copy(), "mse": float(np.mean((y - F) ** 2))})
    for m in range(1, n_trees + 1):
        residual = y - F                      # negative gradient of 1/2 (y-F)^2
        tree = build_tree(X, residual, max_depth)
        F = F + learning_rate * tree_predict(tree, X)
        trees.append(tree)
        if record:
            history.append({"round": m, "F": F.copy(),
                            "residual": (y - F).copy(),
                            "mse": float(np.mean((y - F) ** 2))})
    model = {"F0": F0, "trees": trees, "learning_rate": learning_rate}
    return (model, history) if record else model
# endregion


# region: predict
def gb_predict(model, X):
    """Predict with the boosted model: the initial constant plus every shrunken
    tree. F(x) = F0 + nu * sum_m tree_m(x)."""
    F = np.full(X.shape[0], model["F0"])
    for tree in model["trees"]:
        F = F + model["learning_rate"] * tree_predict(tree, X)
    return F
# endregion


# region: metrics
def mse(y_true, y_pred):
    """Mean squared error: the average squared gap between guess and truth."""
    return float(np.mean((y_true - y_pred) ** 2))


def r2_score(y_true, y_pred):
    """Coefficient of determination: the fraction of variance the model explains.

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


def load_diabetes_csv(path="../data/diabetes.csv"):
    """Diabetes regression set: 442 patients, 10 standardized features, a
    disease-progression score one year on. Source: sklearn's load_diabetes."""
    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 hace un booster a mano en producción. El GradientBoostingRegressor de scikit-learn es el mismo algoritmo —arranca en la media, ajusta árboles poco profundos a los residuales, suma cada uno de regreso con un paso de learning rate— escrito en C y cableado con las perillas que importan:

def sk_gbr(X_train, y_train, X_test, y_test, n_trees, learning_rate, max_depth):
    """Fit a gradient-boosting regressor and return (train_mse, test_mse,
    test_r2, model). loss="squared_error" and init="mean" are the defaults, so
    this is the library twin of fit_gradient_boost — same objective, same
    residual-fitting loop, same shrinkage."""
    model = GradientBoostingRegressor(
        n_estimators=n_trees, learning_rate=learning_rate,
        max_depth=max_depth, random_state=0,
    )
    model.fit(X_train, y_train)
    train_pred = model.predict(X_train)
    test_pred = model.predict(X_test)
    train_mse = float(((train_pred - y_train) ** 2).mean())
    test_mse = float(((test_pred - y_test) ** 2).mean())
    test_r2 = float(model.score(X_test, y_test))
    return train_mse, test_mse, test_r2, model

loss="squared_error" y la inicialización en la media son los valores por defecto, así que este es un gemelo fiel de nuestro fit_gradient_boost: mismo objetivo, mismo ciclo de ajuste al residual, mismo shrinkage. Las diferencias son refinamientos, no otra idea. sklearn parte con friedman_mse por defecto —una pequeña mejora sobre el split simple de reducción de varianza que usamos nosotros, que pondera el score del split por el tamaño de los hijos— y agrega la maquinaria que los problemas reales necesitan: subsample para stochastic gradient boosting (ajustar cada árbol sobre una fracción aleatoria de las filas, lo que las descorrelaciona un poco como hace un bosque), n_iter_no_change para early stopping, y losses más allá del squared error. El mismo motor, más instrumentación.

Desde cero contra la librería

Ajusta ambos con los propios defaults de scikit-learn —100 árboles, learning rate 0.1, árboles de profundidad 3— sobre la mitad de entrenamiento de 309 pacientes de diabetes, y califica sobre los 133 apartados. Y para argumentar a favor del boosting en sí, pon junto a ellos un solo árbol de regresión de profundidad 3, un aprendiz débil solitario:

Los dos boosters caen juntos: nuestro ensamble desde cero saca 0.3766 de R² de prueba y el de sklearn 0.3916, una diferencia de quince milésimas que se debe por completo a que los splits friedman_mse de sklearn eligen umbrales ligeramente distintos a nuestros splits de squared error puro. El mismo algoritmo, esencialmente el mismo ajuste, que es el resultado que quieres: nuestro ciclo de residuales escrito a mano y la implementación en C de sklearn coinciden. En términos de error es un MSE de prueba de 3332.69 para el nuestro contra 3252.65 para el de sklearn, ambos muy por debajo del 5349.11 de la baseline de la media.

Ahora el número que argumenta a favor del boosting. Un solo árbol de profundidad 3 —un aprendiz débil, la cosa que impulsamos— saca 0.2564 de R² en este test set. Apilar cien de ellos ajustando residuales lo llevó a 0.3766, más de la mitad extra de varianza explicada, con los mismos árboles poco profundos y sin ningún tipo de modelo nuevo. Ese es todo el argumento: un árbol de profundidad 3 es un modelo mediocre por su cuenta, y cien de ellos corrigiéndose entre sí son uno fuerte. Una advertencia honesta, la misma que persigue al set de diabetes en todas partes: un R² de 0.38 no es un triunfo, es simplemente lo que hay que sacar con datos médicos difíciles y ruidosos. El punto no es la altitud, es la subida de 0.26 a 0.38 que el boosting compró sobre el árbol solo, y el hecho de que nuestro código desde cero hizo esa subida a la par de la librería.

Conclusiones

El gradient boosting es la idea de AdaBoost, generalizada y hecha honesta. Donde AdaBoost re-ponderaba los ejemplos difíciles, el gradient boosting nombra con precisión qué significa "difícil" —el gradiente negativo de una loss— y ajusta el siguiente árbol directo hacia él. Para el squared error ese gradiente es el residual, así que el método entero se reduce al ciclo más simple imaginable: adivina la media, ajusta un árbol poco profundo a lo que sobra, súmale una rebanada encogida, repite. Cada línea del ciclo desde cero mapea a una línea de las matemáticas, y la palabra "gradient" se gana su lugar porque cambiar la loss cambia exactamente una fórmula y nada más.

Échale mano cuando vayas en serio con un problema tabular. Un random forest es la baseline que ajustas primero para ver si hay señal; un ensamble con gradient boosting es lo que ajustas cuando quieres ganar, porque construir árboles que se corrigen entre sí baja el error más de lo que promediar árboles independientes podrá jamás. El precio es el tuneo: el learning rate y la cantidad de árboles están acoplados, los rates más chicos necesitan más árboles pero llegan a un piso más bajo y te perdonan más, y a diferencia de un bosque va a sobreajustar si sigues agregando árboles pasado el punto donde el error de validación da la vuelta. Pon el learning rate chico, la cantidad de árboles generosa, y deja que el early stopping encuentre la pared. Nunca lo mandes a producción con learning rate 0.5.

Lo que construimos es la gradient boosting machine de 2001 en su ropa más simple: splits de squared error puro, un árbol a la vez, desde cero. Lo que se usa hoy es ese mismo esqueleto con mejores splits, gradientes de segundo orden, regularización sobre los valores de las hojas, submuestreo de columnas y de filas, y un esfuerzo de ingeniería vaciado en hacerlo rápido con millones de filas. Eso es XGBoost, y es el siguiente capítulo: no tanto una idea nueva, sino esta misma idea, afilada hasta convertirse en el modelo que encabeza más leaderboards tabulares que ningún otro.