58  Gradient Boosting y XGBoost

XGBoost (eXtreme Gradient Boosting) fue propuesto por Chen y Guestrin en 2016 (Chen y Guestrin 2016) y desde entonces ha dominado competencias de machine learning en Kaggle y otros benchmarks de referencia. Se trata de una implementación regularizada y optimizada a segundo orden del algoritmo de gradient boosting (Friedman 2001), con innovaciones clave en la búsqueda aproximada de particiones, el manejo eficiente de datos dispersos, y un diseño de sistema orientado a la escalabilidad. Este capitulo cubre los fundamentos matemáticos del algoritmo, sus hiperparámetros principales, la interpretación de importancia de variables mediante SHAP, y la implementación práctica en Python.

1 Intuición: cada árbol corrige al anterior

Imagina que estimas el precio de una casa. Un primer tasador da una cifra; un segundo no tasa la casa desde cero, sino que estudia cuánto se equivocó el primero y propone una corrección; un tercero corrige lo que aún falta, y así sucesivamente. Si cada corrección es pequeña y prudente, la suma converge a una buena estimación. Eso es el gradient boosting: una suma de árboles pequeños en la que cada uno aprende los residuos que dejaron los anteriores. El capítulo anterior presentó AdaBoost, que repesa ejemplos; aquí el mecanismo es más general: cada árbol nuevo se ajusta al gradiente de la pérdida, lo que permite usar cualquier pérdida diferenciable, desde el error cuadrático hasta la log-loss de la clasificación.

2 Construcción paso a paso

2.1 Gradient boosting a mano

Cinco casas con un tamaño \(x\) y un precio \(y\) (en millones):

\(x\)	1	2	3	4	5
\(y\)	2	4	9	10	15

Tabla 1: Cinco observaciones para construir gradient boosting a mano.

Usaremos tocones de regresión (árboles de un solo corte) y una tasa de aprendizaje \(\eta = 0.5\).

Paso 1. Predicción inicial. Con pérdida cuadrática, la mejor constante es el promedio: \(F_0 = 40/5 = 8\). Su suma de cuadrados es \(\text{SSE}_0 = 36 + 16 + 1 + 4 + 49 = 106\).

Paso 2. Residuos y primer árbol. Los residuos \(r_i = y_i - F_0\) son \((-6, -4, 1, 2, 7)\). El mejor tocón para esos residuos corta entre \(x = 2\) y \(x = 3\): a la izquierda predice el promedio \(-5\) y a la derecha \((1 + 2 + 7)/3 = 3.33\).

Paso 3. Encoger y sumar. No sumamos la corrección completa, sino la mitad: \[F_1(x) = F_0 + \eta\,h_1(x) = \begin{cases} 8 + 0.5(-5) = 5.5 & x \le 2.5,\\ 8 + 0.5(3.33) = 9.67 & x > 2.5.\end{cases}\] La suma de cuadrados baja a \(\text{SSE}_1 = 43.5\).

Paso 4. Repetir. Los nuevos residuos son \((-3.5, -1.5, -0.67, 0.33, 5.33)\): el gran error que queda es la última casa. El segundo tocón corta entre \(x = 4\) y \(x = 5\), con valores \(-1.33\) y \(5.33\), y \(F_2 = (4.83, 4.83, 9, 9, 12.33)\), con \(\text{SSE}_2 = 16.83\).

Ronda	Predicciones \(F_m\)	Corte del tocón	SSE
0	8, 8, 8, 8, 8	—	106.00
1	5.50, 5.50, 9.67, 9.67, 9.67	\(x \le 2.5\)	43.50
2	4.83, 4.83, 9.00, 9.00, 12.33	\(x \le 4.5\)	16.83
3	3.92, 3.92, 9.61, 9.61, 12.94	\(x \le 2.5\)	8.43

Tabla 2: Tres rondas de gradient boosting con η = 0.5.

Dos observaciones importantes. Primero, un mismo corte puede reaparecer (rondas 1 y 3): cada ronda solo ve los residuos actuales. Segundo, con \(\eta = 1\) el primer tocón corregiría por completo su parte del error, pero también aprendería por completo su ruido; con \(\eta\) pequeño hacen falta más rondas, a cambio de un modelo más estable. Esa es la contracción (shrinkage).

2.2 Lo que agrega XGBoost: gradiente, hessiano y regularización

Con pérdida cuadrática \(\frac{1}{2}(y - F)^2\), el gradiente respecto a la predicción es \(g_i = F - y_i = -r_i\) y la segunda derivada (hessiano) es \(h_i = 1\). En la primera ronda, \(g = (6, 4, -1, -2, -7)\). XGBoost calcula el valor de cada hoja con la Ec. 58.12, \(w = -G/(H + \lambda)\), donde \(G\) y \(H\) suman \(g_i\) y \(h_i\) de los ejemplos de la hoja. Con \(\lambda = 1\) y el corte \(x \le 2.5\):

\[ w_L = -\frac{6 + 4}{2 + 1} = -3.33, \qquad w_R = -\frac{-1 - 2 - 7}{3 + 1} = 2.5. \]

Comparados con los promedios del tocón clásico (\(-5\) y \(3.33\)), los valores se encogen hacia cero, más en la hoja pequeña: \(\lambda\) actúa como una muestra ficticia que tira cada hoja hacia la predicción actual. La calidad del corte se mide con la ganancia de la Ec. 58.14:

\[ \text{Ganancia} = \frac{1}{2}\left[\frac{10^2}{2 + 1} + \frac{(-10)^2}{3 + 1} - \frac{0^2}{5 + 1}\right] - \gamma = \frac{1}{2}(33.3 + 25) - \gamma = 29.2 - \gamma. \]

Corte	\(G_L\)	\(G_R\)	Ganancia (\(\lambda = 1\), \(\gamma = 0\))
\(x \le 1.5\)	6	−6	12.60
\(x \le 2.5\)	10	−10	29.17
\(x \le 3.5\)	9	−9	23.63
\(x \le 4.5\)	7	−7	17.15

Tabla 3: Ganancia de cada corte en la primera ronda de XGBoost.

Si \(\gamma\) fuera mayor que 29.17, ningún corte compensaría y el árbol se quedaría en una sola hoja: así poda XGBoost.

2.3 Clasificación: el mismo mecanismo con la log-loss

En clasificación binaria, el modelo suma árboles en la escala de los log-odds y la probabilidad es \(p = \sigma(F)\). Para la log-loss, \(g_i = p_i - y_i\) y \(h_i = p_i(1 - p_i)\). Si una hoja recibe tres ejemplos con \(y = (1, 1, 0)\) y la predicción actual es \(p = 0.5\) para todos, entonces \(G = (-0.5) + (-0.5) + 0.5 = -0.5\), \(H = 3 \times 0.25 = 0.75\) y, con \(\lambda = 1\), \(w = 0.5/1.75 = 0.286\): los log-odds de esa hoja suben de 0 a \(0.286\,\eta\). El hessiano hace que las hojas donde el modelo ya está seguro (\(p\) cercano a 0 o 1, \(h\) pequeño) reciban correcciones más cautelosas en relación con su gradiente.

El siguiente código reproduce la Tabla 58.2 con una implementación desde cero y verifica con xgboost los pesos de hoja de la primera ronda:

import numpy as np
import pandas as pd
import xgboost as xgb

x5 = np.array([1, 2, 3, 4, 5], dtype=float)
y5 = np.array([2, 4, 9, 10, 15], dtype=float)

def mejor_tocon(x, r):
    """Tocón de regresión: corte que minimiza la suma de cuadrados de los residuos."""
    mejor = None
    for k in range(1, len(x)):
        izq, der = r[:k], r[k:]
        sse = ((izq - izq.mean()) ** 2).sum() + ((der - der.mean()) ** 2).sum()
        if mejor is None or sse < mejor[0]:
            mejor = (sse, (x[k - 1] + x[k]) / 2, izq.mean(), der.mean())
    return mejor[1:]

F, eta, filas = np.full(5, y5.mean()), 0.5, []
filas.append({"ronda": 0, "predicciones": np.round(F, 2), "corte": None, "SSE": ((y5 - F) ** 2).sum()})
for m in range(1, 4):
    corte, a, b = mejor_tocon(x5, y5 - F)
    F = F + eta * np.where(x5 <= corte, a, b)
    filas.append({"ronda": m, "predicciones": np.round(F, 2), "corte": corte, "SSE": round(((y5 - F) ** 2).sum(), 2)})
print(pd.DataFrame(filas).to_string(index=False))

# Primera ronda de XGBoost con λ = 1: pesos de hoja −G/(H+λ)
m1 = xgb.XGBRegressor(n_estimators=1, max_depth=1, learning_rate=1.0, reg_lambda=1.0, base_score=8.0,
                      tree_method="exact", min_child_weight=0).fit(x5.reshape(-1, 1), y5)
print(m1.get_booster().get_dump(with_stats=True)[0])

 ronda                    predicciones  corte    SSE
     0       [8.0, 8.0, 8.0, 8.0, 8.0]    NaN 106.00
     1    [5.5, 5.5, 9.67, 9.67, 9.67]    2.5  43.50
     2   [4.83, 4.83, 9.0, 9.0, 12.33]    4.5  16.83
     3 [3.92, 3.92, 9.61, 9.61, 12.94]    2.5   8.43
0:[f0<2.5] yes=1,no=2,missing=1,gain=58.3333321,cover=5
    1:leaf=-3.33333325,cover=2
    2:leaf=2.5,cover=3

Las hojas valen −3.33 y 2.5, como en el cálculo a mano. xgboost reporta la ganancia sin el factor ½: \(58.33 = 2 \times 29.17\).

3 De los árboles de decisión al boosting

3.1 El problema del sesgo-varianza y los modelos individuales

Un árbol de decisión profundo tiene la capacidad de memorizar el conjunto de entrenamiento: al aumentar la profundidad, el modelo se ajusta cada vez más a los datos observados y su varianza aumenta. Por el contrario, un árbol poco profundo (o un “tronco” de un solo nivel) introduce un sesgo alto porque su capacidad expresiva es limitada. Este dilema entre sesgo y varianza es central en el aprendizaje automático y motiva el uso de métodos de ensamble.

Los métodos de ensamble combinan múltiples modelos débiles para obtener un modelo fuerte. En el caso del bagging (e.g., Random Forest (Breiman 2001)), los árboles se entrenan de forma independiente sobre muestras con reemplazo, y su predicción se promedia para reducir la varianza. En el boosting, la estrategia es diferente: los modelos se construyen de manera secuencial, donde cada nuevo árbol aprende de los errores del ensamble anterior, reduciendo progresivamente el sesgo.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(42)
n = 80
X_raw = np.sort(rng.uniform(0, 2 * np.pi, n))
y_raw = np.sin(X_raw) + rng.normal(0, 0.35, n)
X_fit = X_raw.reshape(-1, 1)

X_plot = np.linspace(0, 2 * np.pi, 300).reshape(-1, 1)

profundidades = [1, 3, 10]
titulos = [
    "Arbol profundidad 1 (alto sesgo)",
    "Arbol profundidad 3 (equilibrado)",
    "Arbol profundidad 10 (alta varianza)"
]

fig = make_subplots(rows=1, cols=3, subplot_titles=titulos)

for idx, (depth, titulo) in enumerate(zip(profundidades, titulos), start=1):
    tree = DecisionTreeRegressor(max_depth=depth, random_state=0)
    tree.fit(X_fit, y_raw)
    y_pred = tree.predict(X_plot)

    fig.add_trace(
        go.Scatter(
            x=X_raw, y=y_raw,
            mode="markers",
            marker=dict(color="#94a3b8", size=5),
            name="Datos",
            showlegend=(idx == 1)
        ),
        row=1, col=idx
    )
    fig.add_trace(
        go.Scatter(
            x=X_plot.flatten(), y=y_pred,
            mode="lines",
            line=dict(color="#2563eb", width=2.5),
            name="Prediccion",
            showlegend=(idx == 1)
        ),
        row=1, col=idx
    )
    fig.add_trace(
        go.Scatter(
            x=X_plot.flatten(), y=np.sin(X_plot.flatten()),
            mode="lines",
            line=dict(color="black", dash="dash", width=1.5),
            name="Funcion verdadera",
            showlegend=(idx == 1)
        ),
        row=1, col=idx
    )

fig.update_layout(
    height=380,
    template="plotly_white",
    legend=dict(orientation="h", yanchor="bottom", y=-0.25, xanchor="center", x=0.5)
)
fig.update_xaxes(title_text="x")
fig.update_yaxes(title_text="y", col=1)
fig.show()

Figura 1: Ilustración del compromiso sesgo-varianza para árboles de decisión de distinta profundidad. Un árbol superficial (profundidad 1) presenta alto sesgo, mientras que uno muy profundo (profundidad 10) presenta alta varianza.

La Figura 58.1 muestra claramente como la profundidad del árbol determina el equilibrio entre sesgo y varianza. El árbol de profundidad 3 logra un ajuste razonable sin sobreajustar el ruido.

3.2 Boosting: aprendizaje secuencial de residuos

El boosting construye un modelo aditivo de forma secuencial. Partiendo de una predicción inicial constante:

\[F_0(\mathbf{x}) = \underset{\gamma}{\mathrm{argmin}} \sum_{i=1}^{n} L(y_i, \gamma) \tag{1}\]

En cada iteración \(m = 1, \ldots, M\), se calcula el gradiente negativo (pseudoresiduos) de la función de pérdida evaluada en la predicción actual:

\[r_i^{(m)} = -\left[\frac{\partial L(y_i, F(\mathbf{x}_i))}{\partial F(\mathbf{x}_i)}\right]_{F = F_{m-1}} \tag{2}\]

Luego se ajusta un modelo débil \(h_m\) a estos pseudoresiduos, y se actualiza el ensamble:

\[F_m(\mathbf{x}) = F_{m-1}(\mathbf{x}) + \eta \cdot h_m(\mathbf{x}) \tag{3}\]

donde \(\eta \in (0, 1]\) es la tasa de aprendizaje o shrinkage. Para la pérdida cuadrática \(L(y, F) = \tfrac{1}{2}(y - F)^2\), el gradiente negativo es simplemente el residuo ordinario \(r_i = y_i - F_{m-1}(\mathbf{x}_i)\). Esto significa que el algoritmo ajusta iterativamente los errores del ensamble: cada nuevo árbol aprende lo que el ensamble anterior no pudo capturar.

Gradient Boosting (versión simplificada)

1. Inicializar \(F_0(\mathbf{x}) = \bar{y}\) (media de la variable objetivo).
2. Para \(m = 1, \ldots, M\):
   • Calcular residuos \(r_i^{(m)} = y_i - F_{m-1}(\mathbf{x}_i)\).
   • Ajustar un árbol débil \(h_m\) a \(\{(\mathbf{x}_i, r_i^{(m)})\}_{i=1}^n\).
   • Actualizar \(F_m(\mathbf{x}) = F_{m-1}(\mathbf{x}) + \eta \cdot h_m(\mathbf{x})\).
3. Devolver \(F_M(\mathbf{x})\).

La siguiente figura ilustra la convergencia del algoritmo de gradient boosting implementado desde cero, sin usar XGBoost.

import numpy as np
import plotly.graph_objects as go

rng2 = np.random.default_rng(7)
n2 = 60
X2 = np.sort(rng2.uniform(0, 2 * np.pi, n2)).reshape(-1, 1)
y2 = np.sin(X2.flatten()) + rng2.normal(0, 0.3, n2)
X_line = np.linspace(0, 2 * np.pi, 300).reshape(-1, 1)
y_true_line = np.sin(X_line.flatten())


class Stump:
    def fit(self, X, r):
        best_feat, best_thr, best_loss = 0, 0, float("inf")
        for feat in range(X.shape[1]):
            thresholds = np.unique(X[:, feat])
            for thr in thresholds:
                left = r[X[:, feat] <= thr]
                right = r[X[:, feat] > thr]
                if len(left) == 0 or len(right) == 0:
                    continue
                loss = np.sum((left - left.mean()) ** 2) + np.sum((right - right.mean()) ** 2)
                if loss < best_loss:
                    best_loss = loss
                    best_feat, best_thr = feat, thr
                    self.val_left, self.val_right = left.mean(), right.mean()
        self.feat, self.thr = best_feat, best_thr

    def predict(self, X):
        return np.where(X[:, self.feat] <= self.thr, self.val_left, self.val_right)


eta = 0.5
M_max = 20
stumps = []
F_current = np.full(n2, y2.mean())

snapshots = {}
rondas = [1, 2, 3, 5, 10, 20]

for m in range(1, M_max + 1):
    residuos = y2 - F_current
    stump = Stump()
    stump.fit(X2, residuos)
    stumps.append(stump)
    F_current = F_current + eta * stump.predict(X2)
    if m in rondas:
        F_line = np.full(300, y2.mean())
        for s in stumps:
            F_line = F_line + eta * s.predict(X_line)
        snapshots[m] = {
            "F_line": F_line.copy(),
            "residuos": (y2 - F_current).copy()
        }

frames = []
for m in rondas:
    snap = snapshots[m]
    frame = go.Frame(
        data=[
            go.Scatter(
                x=X2.flatten(), y=y2,
                mode="markers",
                marker=dict(color="#94a3b8", size=6),
                name="Datos"
            ),
            go.Scatter(
                x=X_line.flatten(), y=y_true_line,
                mode="lines",
                line=dict(color="black", dash="dash", width=1.5),
                name="Funcion verdadera"
            ),
            go.Scatter(
                x=X_line.flatten(), y=snap["F_line"],
                mode="lines",
                line=dict(color="#2563eb", width=2.5),
                name="Prediccion F_m"
            ),
            go.Scatter(
                x=X2.flatten(), y=snap["residuos"],
                mode="markers",
                marker=dict(color="#f97316", size=6, symbol="diamond"),
                name="Residuos"
            )
        ],
        name=str(m)
    )
    frames.append(frame)

snap0 = snapshots[1]
fig2 = go.Figure(
    data=[
        go.Scatter(x=X2.flatten(), y=y2, mode="markers",
                   marker=dict(color="#94a3b8", size=6), name="Datos"),
        go.Scatter(x=X_line.flatten(), y=y_true_line, mode="lines",
                   line=dict(color="black", dash="dash", width=1.5), name="Funcion verdadera"),
        go.Scatter(x=X_line.flatten(), y=snap0["F_line"], mode="lines",
                   line=dict(color="#2563eb", width=2.5), name="Prediccion F_m"),
        go.Scatter(x=X2.flatten(), y=snap0["residuos"], mode="markers",
                   marker=dict(color="#f97316", size=6, symbol="diamond"), name="Residuos")
    ],
    frames=frames
)

steps = []
for m in rondas:
    step = dict(
        method="animate",
        args=[[str(m)], {"frame": {"duration": 600, "redraw": True}, "mode": "immediate"}],
        label=str(m)
    )
    steps.append(step)

fig2.update_layout(
    template="plotly_white",
    height=440,
    xaxis_title="x",
    yaxis_title="y",
    legend=dict(orientation="h", yanchor="bottom", y=-0.3, xanchor="center", x=0.5),
    sliders=[dict(
        active=0,
        currentvalue=dict(prefix="M = ", font=dict(size=14)),
        pad=dict(t=50),
        steps=steps
    )],
    updatemenus=[dict(
        type="buttons",
        showactive=False,
        y=1.1, x=0.5, xanchor="center",
        buttons=[
            dict(label="Reproducir", method="animate",
                 args=[None, {"frame": {"duration": 800, "redraw": True}, "fromcurrent": True}]),
            dict(label="Pausar", method="animate",
                 args=[[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate"}])
        ]
    )]
)
fig2.show()

Figura 2: Convergencia iterativa del gradient boosting sobre un problema de regresión 1D. Cada frame muestra la predicción acumulada del ensamble (azul) y los residuos actuales (naranja) luego de agregar M árboles.

Observando la Figura 58.2 se aprecia como los residuos (puntos naranjas) se reducen progresivamente a medida que se agregan más árboles al ensamble, y la curva de predicción (azul) se acerca a la función verdadera (línea discontinua negra).

Un aspecto fundamental del boosting es que la tasa de aprendizaje \(\eta\) en la Ec. 58.3 actua como un regulador de la velocidad de convergencia. Con \(\eta = 1\) el algoritmo intenta corregir completamente los errores en cada paso, lo que puede llevar a sobreajuste rápido. Con \(\eta\) pequeño se necesitan más iteraciones pero el modelo final suele generalizar mejor. Este efecto se conoce como shrinkage (Friedman 2002) y es uno de los mecanismos de regularización implicitos del boosting.

El gradient boosting en su forma general (Ec. 58.2) puede verse como un descenso de gradiente funcional: en lugar de actualizar los parámetros de un modelo parametrico en el espacio euclidiano, se actualiza la función de predicción en el espacio funcional. Cada nuevo árbol débil \(h_m\) es el “paso de gradiente” que más reduce la pérdida esperada. Esta perspectiva, formalizada por Friedman (2001), unifica una amplia familia de algoritmos de aprendizaje bajo un único marco matemático.

4 El objetivo de XGBoost

4.1 Modelo aditivo y función objetivo regularizada

XGBoost construye el modelo de forma aditiva:

\[\hat{y}_i^{(t)} = \sum_{k=1}^{t} f_k(\mathbf{x}_i), \quad f_k \in \mathcal{F} \tag{4}\]

donde \(\mathcal{F}\) es el espacio de árboles de regresión. En cada ronda \(t\), se busca el árbol \(f_t\) que minimiza el objetivo regularizado:

\[\text{Obj}^{(t)} = \sum_{i=1}^{n} l\!\left(y_i,\, \hat{y}_i^{(t-1)} + f_t(\mathbf{x}_i)\right) + \Omega(f_t) + \text{constante} \tag{5}\]

La regularización \(\Omega\) penaliza tanto el número de hojas como la magnitud de los pesos:

\[\Omega(f) = \gamma T + \frac{1}{2}\lambda\|\mathbf{w}\|^2 \tag{6}\]

donde \(T\) es el número de hojas del árbol \(f_t\) y \(\mathbf{w}\) son los pesos de las hojas. El parámetro \(\gamma\) controla la poda (un árbol solo se expande si la ganancia supera \(\gamma\)) y \(\lambda\) regula el tamaño de los pesos mediante regularización L2.

Innovación clave de XGBoost: A diferencia del gradient boosting clásico (Friedman 2001), XGBoost incluye explícitamente la regularización en el objetivo de optimización y utiliza información de segundo orden (el hessiano) para encontrar pesos óptimos de forma analítica, lo que conduce a una convergencia más rápida y modelos más robustos.

4.2 Aproximación de Taylor de segundo orden

Para poder optimizar el objetivo de la Ec. 58.5 eficientemente, XGBoost aplica una expansión de Taylor de segundo orden alrededor de \(\hat{y}_i^{(t-1)}\):

\[l\!\left(y_i, \hat{y}_i^{(t-1)} + f_t\right) \approx l\!\left(y_i, \hat{y}_i^{(t-1)}\right) + g_i f_t(\mathbf{x}_i) + \frac{1}{2} h_i f_t^2(\mathbf{x}_i) \tag{7}\]

donde se definen:

\[g_i = \frac{\partial\, l(y_i, \hat{y}^{(t-1)})}{\partial\, \hat{y}^{(t-1)}} \quad \text{(gradiente de primer orden)} \tag{8}\]

\[h_i = \frac{\partial^2\, l(y_i, \hat{y}^{(t-1)})}{\partial\, (\hat{y}^{(t-1)})^2} \quad \text{(hessiano — curvatura local)} \tag{9}\]

Descartando los términos constantes respecto a \(f_t\), el objetivo simplificado es:

\[\tilde{\text{Obj}}^{(t)} = \sum_{i=1}^{n}\left[g_i f_t(\mathbf{x}_i) + \frac{1}{2}h_i f_t^2(\mathbf{x}_i)\right] + \gamma T + \frac{1}{2}\lambda\sum_{j=1}^{T}w_j^2 \tag{10}\]

Esta aproximación cuadrática es exacta para la pérdida cuadrática y una buena aproximación local para cualquier pérdida dos veces diferenciable, lo que hace a XGBoost aplicable a una amplia variedad de problemas (regresión, clasificación, ranking, etc.).

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

y_obs = 1.0
F0 = -0.3

F_vals = np.linspace(-2.0, 2.5, 400)
loss_true = 0.5 * (y_obs - F_vals) ** 2

g_i = -(y_obs - F0)
h_i = 1.0

loss_linear = 0.5 * (y_obs - F0) ** 2 + g_i * (F_vals - F0)
loss_quadratic = 0.5 * (y_obs - F0) ** 2 + g_i * (F_vals - F0) + 0.5 * h_i * (F_vals - F0) ** 2

lam = 1.0
delta_vals = np.linspace(-3.0, 3.0, 400)
obj_delta = g_i * delta_vals + 0.5 * (h_i + lam) * delta_vals ** 2
delta_opt = -g_i / (h_i + lam)
obj_opt = g_i * delta_opt + 0.5 * (h_i + lam) * delta_opt ** 2

fig3 = make_subplots(
    rows=1, cols=2,
    subplot_titles=[
        "Curva de perdida y aproximacion de Taylor",
        "Contribucion cuadratica al objetivo (con lambda)"
    ]
)

fig3.add_trace(
    go.Scatter(x=F_vals, y=loss_true, mode="lines",
               line=dict(color="black", width=2.5),
               name="Perdida real L(y, F)"),
    row=1, col=1
)
fig3.add_trace(
    go.Scatter(x=F_vals, y=loss_linear, mode="lines",
               line=dict(color="#f97316", width=1.8, dash="dot"),
               name="Aprox. lineal (1er orden)"),
    row=1, col=1
)
fig3.add_trace(
    go.Scatter(x=F_vals, y=loss_quadratic, mode="lines",
               line=dict(color="#2563eb", width=2.0, dash="dash"),
               name="Aprox. cuadratica (2do orden)"),
    row=1, col=1
)
fig3.add_trace(
    go.Scatter(x=[F0], y=[0.5 * (y_obs - F0) ** 2],
               mode="markers",
               marker=dict(color="red", size=10, symbol="circle"),
               name="Punto de expansion F_0"),
    row=1, col=1
)

fig3.add_trace(
    go.Scatter(x=delta_vals, y=obj_delta, mode="lines",
               line=dict(color="#2563eb", width=2.5),
               name="g*Delta + (1/2)(h+lambda)*Delta^2",
               showlegend=True),
    row=1, col=2
)
fig3.add_trace(
    go.Scatter(
        x=[delta_opt, delta_opt],
        y=[obj_delta.min() - 0.3, obj_delta.max()],
        mode="lines",
        line=dict(color="red", dash="dash", width=1.8),
        name="Delta* optimo = -g/(h+lambda)"
    ),
    row=1, col=2
)
fig3.add_trace(
    go.Scatter(x=[delta_opt], y=[obj_opt],
               mode="markers",
               marker=dict(color="red", size=10),
               name="Minimo analitico",
               showlegend=False),
    row=1, col=2
)
fig3.add_annotation(
    x=delta_opt, y=obj_opt - 0.5,
    text=f"Delta* = {delta_opt:.2f}",
    showarrow=True, arrowhead=2,
    font=dict(size=12, color="red"),
    row=1, col=2
)

fig3.update_xaxes(title_text="F (prediccion)", row=1, col=1)
fig3.update_yaxes(title_text="Perdida L(y, F)", row=1, col=1)
fig3.update_xaxes(title_text="Delta (cambio en prediccion)", row=1, col=2)
fig3.update_yaxes(title_text="Contribucion al objetivo", row=1, col=2)
fig3.update_layout(
    height=420,
    template="plotly_white",
    legend=dict(orientation="h", yanchor="bottom", y=-0.35, xanchor="center", x=0.5)
)
fig3.show()

Figura 3: Aproximación de Taylor de segundo orden de la función de pérdida. Panel izquierdo: curva de pérdida real (negro), aproximación lineal (naranja) y aproximación cuadrática (azul) alrededor del punto de expansión F=-0.3. Panel derecho: contribución cuadrática al objetivo en función del cambio Delta, con el óptimo analítico senalado.

La Figura 58.3 ilustra por que la aproximación cuadrática es útil: permite encontrar el cambio óptimo \(\Delta^* = -g_i / (h_i + \lambda)\) de forma analítica, sin necesidad de búsqueda iterativa. XGBoost utiliza exactamente esta idea para calcular el peso óptimo de cada hoja.

4.3 Pesos óptimos de las hojas

Agrupando las muestras por hoja \(j\) (donde \(I_j = \{i : q(\mathbf{x}_i) = j\}\) es el conjunto de índices de la hoja \(j\)), el objetivo simplificado se puede reescribir como:

\[\tilde{\text{Obj}}^{(t)} = \sum_{j=1}^{T}\left[\left(\sum_{i\in I_j}g_i\right)w_j + \frac{1}{2}\left(\sum_{i\in I_j}h_i + \lambda\right)w_j^2\right] + \gamma T \tag{11}\]

Definiendo \(G_j = \sum_{i\in I_j} g_i\) y \(H_j = \sum_{i\in I_j} h_i\), y minimizando respecto a \(w_j\), se obtiene el peso óptimo:

\[w_j^* = -\frac{G_j}{H_j + \lambda} \tag{12}\]

Sustituyendo en el objetivo:

\[\text{Obj}^* = -\frac{1}{2}\sum_{j=1}^{T}\frac{G_j^2}{H_j + \lambda} + \gamma T \tag{13}\]

Este valor sirve como puntuación de la estructura del árbol: cuanto más negativo sea, mejor es el árbol. Nótese que \(\lambda > 0\) reduce la magnitud de los pesos y tiene un efecto regularizador directo, ya que \(|w_j^*|\) decrece con \(\lambda\).

4.4 Ganancia de una partición

Para decidir si realizar una división, XGBoost compara el objetivo antes y después de la partición. Si un nodo padre \(P\) se divide en hijos izquierdo \(L\) y derecho \(R\), la ganancia de la división es:

\[\text{Gain} = \frac{1}{2}\left[\frac{G_L^2}{H_L+\lambda} + \frac{G_R^2}{H_R+\lambda} - \frac{G_P^2}{H_P+\lambda}\right] - \gamma \tag{14}\]

Si \(\text{Gain} \leq 0\), la división se rechaza (poda). El término \(-\gamma\) penaliza la complejidad del árbol: un \(\gamma\) mayor hace que solo se acepten divisiones con mejoras sustanciales.

import numpy as np
import plotly.graph_objects as go

rng3 = np.random.default_rng(21)
n3 = 150
X3 = np.sort(rng3.uniform(0, 10, n3))
y3 = 2.0 * np.sin(X3 * 0.8) + rng3.normal(0, 0.5, n3)

F_prev = np.full(n3, y3.mean())
g_vals = -(y3 - F_prev)
h_vals = np.ones(n3)

lam3 = 1.0
G_total = g_vals.sum()
H_total = h_vals.sum()
score_parent = G_total ** 2 / (H_total + lam3)

thresholds = np.unique(X3)[1:-1]
gains = []

for thr in thresholds:
    left_mask = X3 <= thr
    right_mask = ~left_mask
    G_L = g_vals[left_mask].sum()
    H_L = h_vals[left_mask].sum()
    G_R = g_vals[right_mask].sum()
    H_R = h_vals[right_mask].sum()
    if H_L + lam3 <= 0 or H_R + lam3 <= 0:
        gains.append(float("nan"))
        continue
    gain = 0.5 * (G_L ** 2 / (H_L + lam3) + G_R ** 2 / (H_R + lam3) - score_parent)
    gains.append(gain)

gains = np.array(gains)
best_idx = np.nanargmax(gains)
best_thr = thresholds[best_idx]
best_gain = gains[best_idx]

fig4 = go.Figure()

fig4.add_trace(go.Scatter(
    x=thresholds, y=gains,
    mode="lines",
    fill="tozeroy",
    fillcolor="rgba(37, 99, 235, 0.15)",
    line=dict(color="#2563eb", width=2.0),
    name="Ganancia"
))
fig4.add_trace(go.Scatter(
    x=[best_thr, best_thr],
    y=[gains.min() - 0.5, best_gain],
    mode="lines",
    line=dict(color="red", dash="dash", width=2.0),
    name=f"Umbral optimo = {best_thr:.2f}"
))
fig4.add_trace(go.Scatter(
    x=[thresholds.min(), thresholds.max()],
    y=[0, 0],
    mode="lines",
    line=dict(color="gray", dash="dash", width=1.5),
    name="Gain = 0 (poda)"
))
fig4.add_annotation(
    x=best_thr, y=best_gain + 0.3,
    text=f"Optimo: thr={best_thr:.2f}, Gain={best_gain:.2f}",
    showarrow=True, arrowhead=2,
    font=dict(size=12, color="red")
)
fig4.update_layout(
    template="plotly_white",
    height=400,
    xaxis_title="Umbral de particion",
    yaxis_title="Ganancia",
    legend=dict(orientation="h", yanchor="bottom", y=-0.3, xanchor="center", x=0.5)
)
fig4.show()

Figura 4: Ganancia de partición en función del umbral para un conjunto de datos sintético 1D. La línea discontinua roja marca el umbral óptimo (máxima ganancia). La línea horizontal gris indica Gain=0, umbral de poda cuando gamma=0.

La Figura 58.4 muestra como la ganancia varia en función del umbral de corte. El algoritmo busca el umbral que maximiza esta función en cada nivel del árbol. Cuando \(\gamma > 0\), cualquier ganancia por debajo de ese valor lleva a rechazar la partición.

4.5 Búsqueda aproximada de particiones (Approximate Split Finding)

Un aspecto crucial del diseño de XGBoost es que la búsqueda exacta del mejor umbral (recorriendo todos los valores únicos de cada variable) puede ser prohibitivamente costosa para conjuntos de datos grandes. Para abordar esto, Chen y Guestrin (2016) proponen un algoritmo de búsqueda aproximada basado en quantiles ponderados.

La idea es proponer un conjunto reducido de candidatos de partición \(\mathcal{S}_k = \{s_{k1}, s_{k2}, \ldots, s_{kl}\}\) para cada variable \(k\), utilizando los quantiles de la distribución de la variable ponderada por los hessianos \(h_i\):

\[\text{rank}_k(z) = \frac{\sum_{(x,h): x_k < z} h}{\sum_{(x,h)} h} \tag{15}\]

Los candidatos se seleccionan tal que \(|\text{rank}_k(s_{k,j}) - \text{rank}_k(s_{k,j+1})| \leq \varepsilon\), donde \(\varepsilon\) controla el número aproximado de candidatos: \(\approx 1/\varepsilon\) por variable. Esto garantiza que el error de aproximación respecto a la solución exacta esta acotado.

El uso de los hessianos \(h_i\) como pesos en la distribución de quantiles no es arbitrario: las observaciones con mayor hessiano (mayor curvatura de la pérdida) tienen más influencia en la selección de los puntos de corte, lo que dirige la búsqueda hacia las regiones del espacio de características donde la ganancia potencial es mayor.

Sparsity-awareness (manejo de datos dispersos)

XGBoost aprende automáticamente la dirección por defecto para los valores faltantes o cero en cada nodo. Durante el entrenamiento, solo se evaluan los valores no faltantes para cada variable; los valores faltantes se asignan a la rama que minimiza el objetivo. Esto hace que XGBoost sea nativamente compatible con matrices dispersas sin necesidad de imputación previa, una ventaja práctica significativa en aplicaciones reales (Chen y Guestrin 2016).

5 Laboratorio: gradient boosting paso a paso

El laboratorio usa los mismos 40 puntos de entrenamiento y 300 de prueba del laboratorio de bagging (\(y = 2\sin x + 0.3x\) más ruido de varianza 0.36). Cada ronda ajusta un árbol de regresión a los gradientes con la ganancia y los pesos de XGBoost (\(\gamma = 0\)); con \(\lambda = 0\) es el gradient boosting clásico. El panel derecho muestra los residuos actuales y el árbol que se agregará en la ronda siguiente.

Laboratorio · Gradient boosting

¿Cuántas rondas y qué tan grande cada paso?

Avanza las rondas y compara tasas de aprendizaje. Observa cuándo el error de prueba deja de bajar.

Rondas M: Tasa de aprendizaje η: Profundidad de cada árbol: Regularización λ:

ECM entrenamiento

ECM prueba

Mejor M en prueba

Error irreducible σ²0.360

Predicción del ensamble F_M

entrenamientofunción realF_M

Residuos y próximo árbol

residuos y − F_Márbol h_M+1 (antes de multiplicar por η)

Error según el número de rondas

entrenamientopruebaσ² irreducibleronda actual

Guía del laboratorio.

1. Con \(M = 0\), el árbol rojo es la primera corrección. ¿Dónde está el mayor residuo y cómo lo atiende el árbol?
2. Avanza hasta \(M = 10\) con \(\eta = 0.3\). Predice antes de mover: ¿qué pasará con \(\eta = 1\)? Compara las curvas de error.
3. Lleva \(M\) a 300 con profundidad 4. ¿En qué ronda alcanza su mínimo el error de prueba? ¿Qué le pasa después?
4. Sube \(\lambda\) a 10. ¿Cómo cambia la altura del árbol rojo y la velocidad de convergencia?
5. Explica por qué en boosting más árboles sí pueden sobreajustar, a diferencia de bagging.

6 Hiperparámetros principales

XGBoost expone numerosos hiperparámetros que permiten controlar la capacidad del modelo, la regularización, el subsampling y la velocidad de convergencia. La siguiente tabla resume los más importantes.

Desde la perspectiva del aprendizaje estadístico, los hiperparámetros pueden agruparse en tres categorías conceptuales:

1. Capacidad del modelo: max_depth, n_estimators, min_child_weight. Determinan cuanta información puede memorizar el modelo. Valores altos de max_depth y n_estimators pueden llevar a sobreajuste.

2. Regularización: gamma, reg_alpha, reg_lambda. Penalizan la complejidad del modelo y reducen la varianza de los pesos de las hojas. Corresponden directamente a los términos de la Ec. 58.6.

3. Aleatoriedad y robustez: subsample, colsample_bytree, learning_rate. Introducen estocasticidad en el proceso de construcción, similar al gradient boosting estocástico (Friedman 2002), y controlan la velocidad de convergencia mediante el shrinkage \(\eta\).

Parámetro	Descripción	Valor típico
n_estimators	Número de árboles en el ensamble	100–1000
max_depth	Profundidad máxima por árbol	3–8
learning_rate (\(\eta\))	Tasa de aprendizaje (shrinkage)	0.01–0.3
subsample	Fracción de filas muestreadas por árbol	0.6–1.0
colsample_bytree	Fracción de columnas por árbol	0.6–1.0
min_child_weight	Mínimo de \(\sum h_i\) en un nodo hoja	1–10
gamma (\(\gamma\))	Ganancia mínima para realizar una partición	0–5
reg_alpha	Regularización L1 sobre pesos	0–1
reg_lambda	Regularización L2 sobre pesos	1 (default)

Tabla 4: Hiperparámetros principales de XGBoost y sus rangos típicos.

El número de árboles (n_estimators) y la tasa de aprendizaje (learning_rate) interactuan directamente: tasas de aprendizaje menores requieren más árboles pero generalmente producen mejores modelos. El subsample y el colsample_bytree introducen aleatoriedad similar al Random Forest, lo que reduce la varianza y el sobreajuste (Friedman 2002). El min_child_weight controla la complejidad de cada árbol al exigir que cada hoja tenga una suma de hessianos mínima.

6.1 Convergencia y early stopping

import numpy as np
import plotly.graph_objects as go
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import xgboost as xgb

data_ca = fetch_california_housing()
X_ca, y_ca = data_ca.data, data_ca.target
X_train, X_val, y_train, y_val = train_test_split(X_ca, y_ca, test_size=0.2, random_state=42)

dtrain = xgb.DMatrix(X_train, label=y_train)
dval = xgb.DMatrix(X_val, label=y_val)

params = {
    "objective": "reg:squarederror",
    "max_depth": 4,
    "learning_rate": 0.1,
    "subsample": 0.8,
    "colsample_bytree": 0.8,
    "reg_lambda": 1.0,
    "seed": 42
}

evals_result = {}
model_conv = xgb.train(
    params,
    dtrain,
    num_boost_round=300,
    evals=[(dtrain, "train"), (dval, "validacion")],
    early_stopping_rounds=20,
    evals_result=evals_result,
    verbose_eval=False
)

train_rmse = evals_result["train"]["rmse"]
val_rmse = evals_result["validacion"]["rmse"]
best_round = model_conv.best_iteration
rounds = list(range(1, len(train_rmse) + 1))

fig5 = go.Figure()

fig5.add_trace(go.Scatter(
    x=list(range(best_round + 1, len(rounds) + 1)),
    y=val_rmse[best_round:],
    fill="tozeroy",
    fillcolor="rgba(220, 38, 38, 0.08)",
    mode="none",
    name="Region de sobreajuste",
    showlegend=True
))
fig5.add_trace(go.Scatter(
    x=rounds, y=train_rmse,
    mode="lines",
    line=dict(color="#2563eb", width=2.0),
    name="RMSE entrenamiento"
))
fig5.add_trace(go.Scatter(
    x=rounds, y=val_rmse,
    mode="lines",
    line=dict(color="#dc2626", width=2.0),
    name="RMSE validacion"
))
fig5.add_vline(
    x=best_round + 1,
    line=dict(color="black", dash="dash", width=1.8),
    annotation_text=f"Early stop: ronda {best_round + 1}",
    annotation_position="top right"
)
fig5.update_layout(
    template="plotly_white",
    height=420,
    xaxis_title="Ronda de boosting",
    yaxis_title="RMSE",
    legend=dict(orientation="h", yanchor="bottom", y=-0.3, xanchor="center", x=0.5)
)
fig5.show()

Figura 5: Curvas de error RMSE en entrenamiento y validación en función del número de rondas de boosting para el conjunto de datos California Housing. La línea vertical indica la ronda de parada anticipada (early stopping). La zona roja sombreada senala la región de sobreajuste potencial.

La Figura 58.5 muestra el comportamiento típico del gradient boosting: el error de entrenamiento decrece monotonamente, mientras el error de validación se estabiliza y puede comenzar a subir si se excede la cantidad óptima de árboles. La técnica de early stopping detiene el entrenamiento cuando el error de validación no mejora durante un número definido de rondas consecutivas.

6.2 Búsqueda de hiperparámetros: profundidad y tasa de aprendizaje

import numpy as np
import plotly.graph_objects as go
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import xgboost as xgb

data_hp = fetch_california_housing()
X_hp, y_hp = data_hp.data, data_hp.target
X_tr_hp, X_va_hp, y_tr_hp, y_va_hp = train_test_split(X_hp, y_hp, test_size=0.2, random_state=42)

depths = [2, 3, 4, 5, 6]
lrs = [0.01, 0.05, 0.1, 0.2, 0.3]

rmse_grid = np.zeros((len(depths), len(lrs)))

for i, depth in enumerate(depths):
    for j, lr in enumerate(lrs):
        model_g = xgb.XGBRegressor(
            n_estimators=100,
            max_depth=depth,
            learning_rate=lr,
            subsample=0.8,
            colsample_bytree=0.8,
            reg_lambda=1.0,
            random_state=42,
            verbosity=0
        )
        model_g.fit(X_tr_hp, y_tr_hp)
        preds = model_g.predict(X_va_hp)
        rmse_grid[i, j] = np.sqrt(np.mean((y_va_hp - preds) ** 2))

best_i, best_j = np.unravel_index(np.argmin(rmse_grid), rmse_grid.shape)

text_vals = [[f"{rmse_grid[i, j]:.3f}" for j in range(len(lrs))] for i in range(len(depths))]

fig6 = go.Figure()

fig6.add_trace(go.Heatmap(
    z=rmse_grid.tolist(),
    x=[str(lr) for lr in lrs],
    y=[str(d) for d in depths],
    colorscale="RdYlGn_r",
    text=text_vals,
    texttemplate="%{text}",
    colorbar=dict(title="RMSE"),
    showscale=True
))

bx = best_j
by = best_i
dx = 0.5
dy = 0.5

fig6.add_trace(go.Scatter(
    x=[bx - dx, bx + dx, bx + dx, bx - dx, bx - dx],
    y=[by - dy, by - dy, by + dy, by + dy, by - dy],
    mode="lines",
    line=dict(color="red", width=3),
    name=f"Mejor: depth={depths[best_i]}, lr={lrs[best_j]}",
    showlegend=True
))

fig6.update_layout(
    template="plotly_white",
    height=440,
    xaxis_title="learning_rate",
    yaxis_title="max_depth",
    legend=dict(orientation="h", yanchor="bottom", y=-0.2, xanchor="center", x=0.5)
)
fig6.show()

Figura 6: RMSE de validación sobre una grilla de max_depth x learning_rate en el conjunto California Housing. Valores más bajos (verde) indican mejor rendimiento. La celda óptima se marca con un borde rojo.

La Figura 58.6 revela la interacción entre la profundidad del árbol y la tasa de aprendizaje. En general, tasas de aprendizaje muy bajas (0.01) requieren más árboles para converger y con solo 100 estimadores pueden tener alto sesgo. Profundidades intermedias (3–5) suelen dar el mejor equilibrio entre capacidad y generalización.

7 Importancia de variables

XGBoost ofrece tres métricas para cuantificar la importancia de cada variable en el modelo entrenado:

• Frecuencia (weight): número de veces que la variable se usa en una división a lo largo de todos los árboles. Es sensible al número de categorías distintas de cada variable.
• Ganancia (gain): reducción promedio de la pérdida aportada por cada división que usa esa variable. Es la métrica más informativa porque captura la contribución real al rendimiento del modelo.
• Cobertura (cover): número promedio de observaciones que pasan por cada división de la variable. Refleja cuantas muestras son “afectadas” por esa variable.

La métrica de ganancia suele ser la más recomendable para interpretar cuales variables son verdaderamente relevantes, ya que la frecuencia puede sobreestimar la importancia de variables con muchos valores posibles.

Es importante destacar que todas estas métricas miden la importancia global del modelo, es decir, el aporte promedio de cada variable a través de todas las observaciones. Para obtener importancias locales (especificas a una predicción individual), se deben usar los valores SHAP, descritos en la Sección 58.9.

Matemáticamente, la importancia de ganancia para la variable \(k\) se calcula como:

\[\text{Gain}_k = \frac{1}{|\mathcal{T}_k|} \sum_{t \in \mathcal{T}_k} \text{Gain}_t \tag{16}\]

donde \(\mathcal{T}_k\) es el conjunto de todas las divisiones que usan la variable \(k\) en todos los árboles del ensamble. La importancia de frecuencia es simplemente \(|\mathcal{T}_k|\), y la de cobertura es el promedio de \(|I_L| + |I_R|\) sobre todas esas divisiones.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import xgboost as xgb

data_imp = fetch_california_housing()
X_imp, y_imp = data_imp.data, data_imp.target
X_tr_imp, X_va_imp, y_tr_imp, y_va_imp = train_test_split(X_imp, y_imp, test_size=0.2, random_state=42)
feat_names = list(data_imp.feature_names)

dtrain_imp = xgb.DMatrix(X_tr_imp, label=y_tr_imp, feature_names=feat_names)
dval_imp = xgb.DMatrix(X_va_imp, label=y_va_imp, feature_names=feat_names)

params_imp = {
    "objective": "reg:squarederror",
    "max_depth": 4,
    "learning_rate": 0.1,
    "subsample": 0.8,
    "colsample_bytree": 0.8,
    "reg_lambda": 1.0,
    "seed": 42
}
evals_r = {}
model_imp = xgb.train(
    params_imp, dtrain_imp,
    num_boost_round=300,
    evals=[(dtrain_imp, "train"), (dval_imp, "validacion")],
    early_stopping_rounds=20,
    evals_result=evals_r,
    verbose_eval=False
)

imp_types = ["weight", "gain", "cover"]
subtitles = ["Frecuencia", "Ganancia", "Cobertura"]

fig7 = make_subplots(rows=1, cols=3, subplot_titles=subtitles)

for col_idx, (itype, subtitle) in enumerate(zip(imp_types, subtitles), start=1):
    scores = model_imp.get_score(importance_type=itype)
    all_feats = feat_names
    vals = [scores.get(f, 0.0) for f in all_feats]
    sorted_idx = np.argsort(vals)
    sorted_feats = [all_feats[i] for i in sorted_idx]
    sorted_vals = [vals[i] for i in sorted_idx]

    fig7.add_trace(
        go.Bar(
            x=sorted_vals,
            y=sorted_feats,
            orientation="h",
            marker=dict(color="#6366f1"),
            name=subtitle,
            showlegend=False
        ),
        row=1, col=col_idx
    )
    fig7.update_xaxes(title_text=subtitle, row=1, col=col_idx)

fig7.update_layout(
    template="plotly_white",
    height=420
)
fig7.show()

Figura 7: Importancia de variables del modelo XGBoost entrenado en California Housing, según tres criterios: Frecuencia (número de divisiones), Ganancia (reducción media de pérdida) y Cobertura (muestras promedio por división).

Según la Figura 58.7, la variable MedInc (ingreso mediano del bloque censal) es consistentemente la más importante en los tres criterios, lo cual es coherente con el conocimiento del dominio: el ingreso es el predictor más fuerte del precio de la vivienda en California. Las variables Latitude y Longitude también aparecen como relevantes, capturando efectos geograficos.

8 Regularización y poda

8.1 Efecto de \(\lambda\) (L2) y \(\gamma\) (ganancia mínima)

La regularización en XGBoost actua por dos mecanismos complementarios, derivados directamente de la Ec. 58.6:

Regularización L2 (\(\lambda\)): El peso óptimo de cada hoja es \(w_j^* = -G_j / (H_j + \lambda)\). Un \(\lambda\) mayor reduce la magnitud de los pesos, produciendo predicciones más conservadoras y previniendo el sobreajuste. Es equivalente a la regularización Ridge en modelos lineales pero aplicada a los pesos de cada hoja del árbol.

Poda por ganancia mínima (\(\gamma\)): Una división solo se realiza si \(\text{Gain} > \gamma\) (ver Ec. 58.14). Un \(\gamma\) mayor lleva a árboles más pequeños y con menos divisiones. A diferencia de la poda post-entrenamiento, esta poda ocurre durante la construcción del árbol, lo que es más eficiente computacionalmente.

La combinación de ambos parámetros es complementaria: \(\lambda\) controla la magnitud de cada predicción de hoja, mientras que \(\gamma\) controla la estructura del árbol. Se recomienda priorizar \(\lambda\) alto sobre \(\gamma\) alto, ya que \(\lambda\) tiene un efecto más suave y continuo.

import numpy as np
import plotly.graph_objects as go
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import xgboost as xgb

data_reg = fetch_california_housing()
X_reg, y_reg = data_reg.data, data_reg.target
X_tr_reg, X_va_reg, y_tr_reg, y_va_reg = train_test_split(X_reg, y_reg, test_size=0.2, random_state=42)

lambdas = [0.01, 0.1, 1, 5, 10]
gammas = [0, 0.5, 1, 2, 5]

rmse_reg = np.zeros((len(lambdas), len(gammas)))

for i, lam in enumerate(lambdas):
    for j, gam in enumerate(gammas):
        m_reg = xgb.XGBRegressor(
            n_estimators=100,
            max_depth=4,
            learning_rate=0.1,
            subsample=0.8,
            colsample_bytree=0.8,
            reg_lambda=lam,
            gamma=gam,
            random_state=42,
            verbosity=0
        )
        m_reg.fit(X_tr_reg, y_tr_reg)
        preds_reg = m_reg.predict(X_va_reg)
        rmse_reg[i, j] = np.sqrt(np.mean((y_va_reg - preds_reg) ** 2))

text_reg = [[f"{rmse_reg[i, j]:.3f}" for j in range(len(gammas))] for i in range(len(lambdas))]

fig8 = go.Figure(go.Heatmap(
    z=rmse_reg.tolist(),
    x=[str(g) for g in gammas],
    y=[str(l) for l in lambdas],
    colorscale="RdYlGn_r",
    text=text_reg,
    texttemplate="%{text}",
    colorbar=dict(title="RMSE")
))

fig8.update_layout(
    title="Validacion RMSE: efecto de lambda y gamma",
    template="plotly_white",
    height=420,
    xaxis_title="gamma (poda minima)",
    yaxis_title="lambda (L2)"
)
fig8.show()

Figura 8: RMSE de validación en función de la regularización L2 (lambda) y la ganancia mínima de poda (gamma) en California Housing. Colores más verdes indican menor error.

La Figura 58.8 confirma que valores muy altos de \(\gamma\) (mayor a 2) pueden degradar el rendimiento al podar demasiado agresivamente, impidiendo que el modelo capture patrones relevantes. Un \(\lambda\) moderado (entre 1 y 5) generalmente ofrece el mejor equilibrio.

Adicionalmente, XGBoost soporta regularización L1 sobre los pesos de las hojas mediante el parámetro reg_alpha. El objetivo con regularización elastica (L1 + L2) toma la forma:

\[\Omega(f) = \gamma T + \alpha \sum_{j=1}^{T}|w_j| + \frac{1}{2}\lambda\sum_{j=1}^{T}w_j^2 \tag{17}\]

La regularización L1 (\(\alpha > 0\)) induce esparsidad en los pesos de las hojas: muchos pesos se vuelven exactamente cero, lo que produce árboles donde algunas hojas no contribuyen a la predicción. Esto puede ser útil cuando se sospecha que muchas divisiones son redundantes o ruidosas. En la práctica, la regularización L2 (\(\lambda\)) es suficiente en la mayoría de los casos y se recomienda comenzar con \(\alpha = 0\).

8.2 Efecto del subsampling sobre la regularización

El parámetro subsample (fracción de filas) y colsample_bytree (fracción de columnas) introducen aleatoriedad en la construcción de cada árbol. Desde una perspectiva estadística, esto puede verse como una forma de regularización implícita: al entrenar cada árbol con un subconjunto distinto de datos, se reduce la correlación entre los árboles del ensamble, lo que disminuye la varianza de las predicciones. Friedman (2002) demostro empíricamente que el subsampling de filas mejora consistentemente la generalización del gradient boosting, especialmente en conjuntos de datos con ruido.

9 Contribuciones individuales (SHAP simplificado)

Los valores SHAP (SHapley Additive exPlanations), propuestos por Lundberg y Lee (2017), permiten descomponer la predicción de cualquier modelo en contribuciones aditivas de cada variable para una observación específica. Para XGBoost, Lundberg et al. (2020) desarrollaron el algoritmo TreeSHAP, que calcula estas contribuciones de forma exacta y eficiente en tiempo polinomial.

Formalmente, para una predicción \(\hat{y}_i\), los valores SHAP \(\phi_j\) satisfacen:

\[\hat{y}_i = \phi_0 + \sum_{j=1}^{p} \phi_j \tag{18}\]

donde \(\phi_0 = \mathbb{E}[\hat{y}]\) es la predicción base (valor esperado del modelo) y \(\phi_j\) es la contribución de la variable \(j\) para esa observación en particular. Las contribuciones pueden ser positivas (la variable empuja la predicción hacia arriba) o negativas (la variable la baja).

import numpy as np
import plotly.graph_objects as go
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import xgboost as xgb

data_wf = fetch_california_housing()
X_wf, y_wf = data_wf.data, data_wf.target
feat_names_wf = list(data_wf.feature_names)

X_tr_wf, X_te_wf, y_tr_wf, y_te_wf = train_test_split(X_wf, y_wf, test_size=0.2, random_state=42)

dtrain_wf = xgb.DMatrix(X_tr_wf, label=y_tr_wf, feature_names=feat_names_wf)

params_wf = {
    "objective": "reg:squarederror",
    "max_depth": 4,
    "learning_rate": 0.1,
    "subsample": 0.8,
    "colsample_bytree": 0.8,
    "reg_lambda": 1.0,
    "seed": 42
}
evals_wf = {}
model_wf = xgb.train(
    params_wf, dtrain_wf,
    num_boost_round=300,
    evals=[(dtrain_wf, "train")],
    evals_result=evals_wf,
    verbose_eval=False
)

obs_idx = 42
X_single = X_te_wf[obs_idx:obs_idx + 1]
d_single = xgb.DMatrix(X_single, feature_names=feat_names_wf)

contribs = model_wf.predict(d_single, pred_contribs=True)[0]
bias = contribs[-1]
feature_contribs = contribs[:-1]

sorted_idx = np.argsort(np.abs(feature_contribs))[::-1]
sorted_feats = [feat_names_wf[i] for i in sorted_idx]
sorted_contribs = feature_contribs[sorted_idx]

cumulative = bias
running = [bias]
for c in sorted_contribs:
    running.append(running[-1] + c)

final_pred = running[-1]

bar_colors = ["#16a34a" if c >= 0 else "#dc2626" for c in sorted_contribs]
labels = ["Base"] + sorted_feats
measures = ["absolute"] + ["relative"] * len(sorted_feats)

fig9 = go.Figure(go.Waterfall(
    orientation="v",
    measure=measures,
    x=labels,
    y=[bias] + sorted_contribs.tolist(),
    connector=dict(line=dict(color="gray", width=1)),
    increasing=dict(marker=dict(color="#16a34a")),
    decreasing=dict(marker=dict(color="#dc2626")),
    totals=dict(marker=dict(color="#6366f1"))
))

fig9.add_annotation(
    x=len(labels) - 1,
    y=final_pred + 0.15,
    text=f"Prediccion final: {final_pred:.3f}",
    showarrow=True,
    arrowhead=2,
    font=dict(size=12, color="black")
)

fig9.update_layout(
    title="Contribucion de cada variable para una prediccion individual",
    template="plotly_white",
    height=480,
    xaxis_title="Variable",
    yaxis_title="Contribucion a la prediccion",
    showlegend=False
)
fig9.show()

Figura 9: Gráfico de cascada (waterfall) que muestra la contribución de cada variable para una predicción individual en el conjunto California Housing. Las barras verdes suman a la predicción y las rojas restan. Se muestran ordenadas por magnitud de contribución.

La Figura 58.9 muestra como la predicción final se construye sumando la contribución de cada variable al valor base. Para esta observación particular, MedInc tiene la mayor influencia positiva, mientras que algunas variables geograficas pueden reducir la predicción estimada. Esta visualización es especialmente útil para comunicar predicciones individuales a usuarios no técnicos o para auditar el comportamiento del modelo.

10 XGBoost para clasificación

XGBoost es igualmente aplicable a problemas de clasificación, con soporte nativo para:

• Clasificación binaria: objective="binary:logistic" — el modelo genera probabilidades mediante la función sigmoide \(p_i = \sigma(\hat{y}_i)\), y la función de pérdida es la log-verosimilitud negativa (log-loss). Los gradientes toman la forma \(g_i = p_i - y_i\) y el hessiano \(h_i = p_i(1-p_i)\).
• Clasificación multiclase: objective="multi:softprob" con num_class=K — se entrena un conjunto de árboles por clase y se aplica softmax para normalizar las probabilidades.

La ventaja de usar el hessiano \(h_i = p_i(1-p_i)\) es que pesa menos las observaciones donde el modelo ya esta seguro (probabilidad cercana a 0 o 1), concentrando el aprendizaje en los casos difíciles (probabilidad cercana a 0.5).

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
import xgboost as xgb

data_clf = load_breast_cancer()
X_c, y_c = data_clf.data, data_clf.target
X_tr_c, X_te_c, y_tr_c, y_te_c = train_test_split(X_c, y_c, test_size=0.2, random_state=42)

clf = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=4,
    learning_rate=0.1,
    subsample=0.8,
    colsample_bytree=0.8,
    eval_metric="logloss",
    random_state=42,
    verbosity=0
)
clf.fit(X_tr_c, y_tr_c)
auc = roc_auc_score(y_te_c, clf.predict_proba(X_te_c)[:, 1])
acc = (clf.predict(X_te_c) == y_te_c).mean()
print(f"AUC-ROC en prueba: {auc:.4f}")
print(f"Exactitud en prueba: {acc:.4f}")

AUC-ROC en prueba: 0.9925
Exactitud en prueba: 0.9561

Tabla 5

El resultado muestra el excelente rendimiento de XGBoost incluso con hiperparámetros predeterminados. En el conjunto de cancer de mama (569 observaciones, 30 variables), se obtiene un AUC-ROC superior a 0.99, lo cual refleja tanto la calidad del clasificador como la separabilidad del conjunto de datos.

10.1 Curva ROC y análisis de umbrales

La curva ROC (Receiver Operating Characteristic) es una herramienta estándar para evaluar clasificadores binarios. Traza la tasa de verdaderos positivos (sensibilidad) frente a la tasa de falsos positivos (1 - especificidad) para todos los umbrales de decisión posibles. El área bajo la curva (AUC) cuantifica la capacidad discriminatoria del modelo en un único número: AUC = 1 indica clasificación perfecta, AUC = 0.5 indica predicción aleatoria.

import numpy as np
import plotly.graph_objects as go
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_curve, roc_auc_score
import xgboost as xgb

data_roc = load_breast_cancer()
X_roc, y_roc = data_roc.data, data_roc.target
X_tr_roc, X_te_roc, y_tr_roc, y_te_roc = train_test_split(X_roc, y_roc, test_size=0.2, random_state=42)

clf_roc = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=4,
    learning_rate=0.1,
    subsample=0.8,
    colsample_bytree=0.8,
    eval_metric="logloss",
    random_state=42,
    verbosity=0
)
clf_roc.fit(X_tr_roc, y_tr_roc)
probs = clf_roc.predict_proba(X_te_roc)[:, 1]
fpr, tpr, thresholds_roc = roc_curve(y_te_roc, probs)
auc_val = roc_auc_score(y_te_roc, probs)

fig_roc = go.Figure()

fig_roc.add_trace(go.Scatter(
    x=[0, 1], y=[0, 1],
    mode="lines",
    line=dict(color="gray", dash="dash", width=1.5),
    name="Clasificador aleatorio (AUC = 0.5)"
))

fig_roc.add_trace(go.Scatter(
    x=fpr.tolist(), y=tpr.tolist(),
    mode="lines",
    fill="tozeroy",
    fillcolor="rgba(37, 99, 235, 0.12)",
    line=dict(color="#2563eb", width=2.5),
    name=f"XGBoost (AUC = {auc_val:.4f})"
))

fig_roc.add_annotation(
    x=0.6, y=0.35,
    text=f"AUC = {auc_val:.4f}",
    showarrow=False,
    font=dict(size=16, color="#2563eb"),
    bgcolor="white",
    bordercolor="#2563eb",
    borderwidth=1
)

fig_roc.update_layout(
    template="plotly_white",
    height=440,
    xaxis_title="Tasa de falsos positivos (1 - Especificidad)",
    yaxis_title="Tasa de verdaderos positivos (Sensibilidad)",
    xaxis=dict(range=[0, 1]),
    yaxis=dict(range=[0, 1.02]),
    legend=dict(orientation="h", yanchor="bottom", y=-0.25, xanchor="center", x=0.5)
)
fig_roc.show()

Figura 10: Curva ROC del clasificador XGBoost en el conjunto de prueba de cancer de mama. El área bajo la curva (AUC) indica la capacidad discriminatoria del modelo. La línea diagonal punteada representa un clasificador aleatorio.

La Figura 58.10 confirma visualmente el alto rendimiento del clasificador: la curva ROC se acerca al punto ideal (0, 1), y el área bajo la curva es prácticamente perfecta para este conjunto de datos. En aplicaciones reales, donde la separabilidad es menor, la curva ROC permite elegir el umbral de decisión adecuado según el costo relativo de falsos positivos y falsos negativos.

11 Buenas prácticas y comparación con alternativas

11.1 Comparación de métodos de ensamble

La siguiente tabla resume el posicionamiento de XGBoost en el contexto de otros métodos de aprendizaje supervisado.

Método	Complejidad	Interpretabilidad	Rendimiento típico
Árbol único	Baja	Alta	Bajo
Random Forest	Media	Media	Bueno
Gradient Boosting (sklearn)	Media	Media	Muy bueno
XGBoost	Media-alta	Media (+ SHAP)	Excelente
LightGBM / CatBoost	Media-alta	Media (+ SHAP)	Excelente

Tabla 6: Comparación de métodos de ensamble según complejidad, interpretabilidad y rendimiento típico.

11.2 Recomendaciones prácticas

Las siguientes guias resumen las mejores prácticas para usar XGBoost de forma efectiva:

1. Estrategia de ajuste de hiperparámetros: Comenzar con learning_rate=0.1 y n_estimators=200–500. Una vez identificado un buen número de árboles mediante early stopping, reducir la tasa de aprendizaje a 0.05 o 0.01 y aumentar proporcionalmente los estimadores para refinar el modelo (Hastie et al. 2009).

2. Early stopping obligatorio: Siempre usar early_stopping_rounds con un conjunto de validación separado. Esto evita el sobreajuste sin necesidad de una búsqueda exhaustiva del número óptimo de árboles.

3. Subsampling para reducir varianza: Activar subsample=0.8 y colsample_bytree=0.8 por defecto. La aleatoriedad adicional que introducen, similar al Random Forest, generalmente mejora la generalización (Friedman 2002).

4. Métrica de importancia: Preferir importance_type="gain" sobre "weight", ya que la frecuencia puede ser enganiosa cuando las variables tienen cardinalidades muy distintas.

5. Regularización: Mantener reg_lambda=1 (valor por defecto) como punto de partida. Aumentar si hay seniales de sobreajuste. Usar gamma con moderación (valores entre 0 y 1) para podar divisiones poco informativas.

6. Validación cruzada: Para conjuntos de datos pequeños, usar xgb.cv() con k-fold para estimar el número óptimo de árboles de forma más robusta que con un único conjunto de validación.

7. Manejo de valores faltantes: XGBoost aprende automáticamente la dirección óptima para los valores faltantes durante el entrenamiento (sparsity-awareness (Chen y Guestrin 2016)). No es necesario imputar previamente.

8. Comparar con LightGBM: Para conjuntos de datos muy grandes (millones de filas o cientos de variables), LightGBM suele ser más rápido gracias a su estrategia de construcción de árboles basada en histogramas (leaf-wise growth). La API es muy similar a XGBoost.

Referencias

Breiman, Leo. 2001. «Random Forests». Machine Learning 45 (1): 5-32. https://doi.org/10.1023/A:1010933404324.

Chen, Tianqi, y Carlos Guestrin. 2016. «XGBoost: A Scalable Tree Boosting System». Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 785-94. https://doi.org/10.1145/2939672.2939785.

Friedman, Jerome H. 2001. «Greedy Function Approximation: A Gradient Boosting Machine». Annals of Statistics 29 (5): 1189-232. https://doi.org/10.1214/aos/1013203451.

Friedman, Jerome H. 2002. «Stochastic Gradient Boosting». Computational Statistics & Data Analysis 38 (4): 367-78. https://doi.org/10.1016/S0167-9473(01)00065-2.

Hastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2.ª ed. Springer Series en Statistics. Springer. https://doi.org/10.1007/978-0-387-84858-7.

Lundberg, Scott M., Gabriel Erion, Hugh Chen, et al. 2020. «From Local Explanations to Global Understanding with Explainable AI for Trees». Nature Machine Intelligence 2: 56-67. https://doi.org/10.1038/s42256-019-0138-9.

Lundberg, Scott M., y Su-In Lee. 2017. «A Unified Approach to Interpreting Model Predictions». Advances in Neural Information Processing Systems 30.