XGBoost

What this chapter covers

The previous chapter built gradient boosting: grow trees one at a time, each one fit to the negative gradient of the loss, add them up with a small learning rate. It works, and for fifteen years it was the strongest thing you could point at a table of numbers. This chapter is what happened next. In 2016 XGBoost took that same additive idea and changed two things about how each tree gets built, and the result won so many Kaggle competitions that "just use XGBoost" became a running joke that was also correct.

The two changes are the whole chapter. First, XGBoost uses a second-order view of the loss — every sample carries a gradient and a hessian, and the tree minimizes a quadratic approximation of the loss instead of just chasing the gradient. That's a Newton step, not a gradient step. Second, the tree is regularized from the inside: an L2 penalty on the leaf values and a per-leaf complexity penalty are written directly into the split math, so the tree that comes out is already pruned by its own objective. Put those together and you get three clean formulas — the leaf weight, the split gain, and the shrunk update — that we build from scratch in NumPy and check against the real xgboost package.

The centerpiece is an animation: a 1-D staircase that snaps to a curve one boosting round at a time, with a loss panel falling underneath it. Then we run our hand-built booster, the real xgboost, and sklearn's plain gradient boosting on the same diabetes data and watch the regularized second-order model pull ahead.

A bit of history

Gradient boosting is Jerome Friedman's. In two papers around the turn of the millennium — "Greedy Function Approximation: A Gradient Boosting Machine" (1999, published 2001) — he framed boosting as gradient descent in function space: at each step, fit a weak learner to the negative gradient of the loss and take a step in that direction. The same year, Friedman, Trevor Hastie, and Robert Tibshirani wrote LogitBoost, which already used a second-order Newton step for the classification loss. The pieces were all on the table by 2001.

XGBoost is Tianqi Chen's. He started it around 2014 as a research project at the University of Washington, under the Distributed Machine Learning Community, and in 2016 he and Carlos Guestrin published "XGBoost: A Scalable Tree Boosting System" at KDD. The paper's contribution is partly the algorithm — the regularized objective and the second-order split gain we build below — and partly an engineering feat: a sparsity-aware split finder, cache-aware access patterns, out-of-core computation, and a way to propose split candidates without scanning every value. The combination mattered. The paper reports that of 29 winning solutions posted on Kaggle in 2015, 17 used XGBoost, and every top-10 team in the KDD Cup 2015 used it. For a stretch of years, if the data was a table, the winning model was this one. LightGBM and CatBoost later filed off some of its edges, but they're variations on the same theme Chen set: regularized, second-order tree boosting, built to scale.

The intuition

Start with a picture of what boosting is doing. Here's a small 1-D set — one input x, one output y, ninety points, a smooth curve buried under noise. Cyan is the training data the model sees; amber diamonds are a held-out validation set it doesn't.

Boosting fits this curve by addition. Start with a flat guess — the average of y. It's wrong everywhere, and the amount it's wrong by, point by point, is the residual. Fit a shallow tree to those residuals, add a small fraction of it to the guess, and the guess bends a little toward the data. Now there's a new, smaller residual. Fit another tree to that. Keep going. Each tree is a small correction to the sum of all the trees before it, and the running prediction staircases toward the curve.

Plain gradient boosting stops the story there: the residual is the negative gradient of squared-error loss, so "fit the residual" is "fit the gradient." XGBoost asks a sharper question. It doesn't just want to know which direction to step — it wants to know how far, and for that it needs the curvature of the loss, the hessian. A sample where the loss is steep and sharply curved should pull the leaf toward it harder than a sample on a shallow, flat part of the loss. The hessian is what carries that information, and it's the thing plain gradient boosting throws away.

The math

Boosting is additive. After $t-1$ rounds the model predicts $\hat{y}^{(t-1)}$, and round $t$ adds one tree $f_t$. XGBoost writes the loss it wants that tree to minimize and then approximates it with a second-order Taylor expansion around the current prediction:

Here $g_i = \partial_{\hat{y}} \ell(y_i, \hat{y}^{(t-1)})$ is the gradient of the loss at sample $i$, and $h_i = \partial^2_{\hat{y}} \ell(y_i, \hat{y}^{(t-1)})$ is the hessian. That $h_i$ in the quadratic term is the entire difference from plain gradient boosting, which keeps only the linear $g_i$ piece. The loss enters only through $g$ and $h$ — change the loss and nothing downstream changes but those two numbers.

The regularizer $\Omega$ is the second new idea. A tree with $T$ leaves and leaf values $w_1, \dots, w_T$ is penalized for both:

$\gamma$ charges a fixed price per leaf, and $\lambda$ is an L2 penalty pulling every leaf value toward zero. Now fix the shape of the tree and ask: what leaf value $w_j$ minimizes the objective? Each leaf $j$ owns a set of samples $I_j$; write $G_j = \sum_{i \in I_j} g_i$ and $H_j = \sum_{i \in I_j} h_i$ for their gradient and hessian sums. The objective in $w_j$ is a simple upward parabola, and its minimum is closed-form:

That is the leaf weight, the single most important formula in the chapter. Read it: the leaf points against the gradient (the minus sign), scaled by the inverse curvature (the hessian in the denominator), and shrunk toward zero by $\lambda$. Drop $\lambda$ and it's the pure Newton step $-G/H$. Substitute $w_j^{*}$ back into the objective and you get the best value a fixed tree can reach, its structure score:

Lower is better, so each leaf wants its $G_j^2/(H_j+\lambda)$ as large as possible. That gives the split rule. When we split one leaf into a left and right child, the change in the structure score is the gain:

The bracket is how much lower the objective goes by giving the two children their own leaf weights instead of one shared weight, and $\gamma$ is the toll for the extra leaf. A split only happens when the structure it buys beats $\gamma$. That subtraction is the whole of XGBoost's pruning — no separate pruning pass, the objective refuses bad splits on its own.

Two losses cover most work. For squared-error regression, $g_i = \hat{y}_i - y_i$ and $h_i = 1$. For binary log-loss with probability $p_i$, $g_i = p_i - y_i$ and $h_i = p_i(1 - p_i)$. That's the only thing that changes between regression and classification. Everything above — leaf weight, gain, pruning — is identical. Finally, the round's tree isn't added whole. Shrinkage scales it by a learning rate $\eta$:

so every tree only moves the prediction part of the way and later trees still have work left. That's a third regularizer, stacked on top of $\lambda$ and $\gamma$.

What it's good at, what it isn't

XGBoost is the model to beat on tabular data, and it earns that by not asking much of you. It inherits everything trees are good at — no feature scaling, mixed numeric and categorical inputs, interactions found for free by nesting splits — and adds the two things a single tree lacks: the boosting ensemble drives down bias by stacking hundreds of corrections, and the regularized objective keeps that ensemble from memorizing the training set. It handles missing values natively by learning a default direction at each split, it trains fast because the split finder was engineered for it, and out of the box, with almost no tuning, it lands near the top of most tabular leaderboards. When the data is a spreadsheet, this is the first thing I reach for.

The costs are real, though. It has a lot of knobs — tree count, depth, learning rate, $\lambda$, $\gamma$, subsampling, column sampling, minimum child weight — and they interact, so getting the last few percent means a tuning budget that a random forest doesn't ask for. It is not interpretable the way one tree is; a pile of three hundred trees is a black box you probe with feature importances and SHAP values, not a flowchart you hand to a domain expert. And it wants tabular data specifically — on images, audio, or raw text, where the structure is spatial or sequential, a neural network eats it alive. XGBoost is a specialist, and its specialty happens to be the most common data shape in the working world.

The data

Two datasets, two jobs. The concept animation uses the 1-D toy from the intuition section — ninety points from a noisy curve, split into sixty-five training points and twenty-five held out — small enough that you can watch the staircase land on points you can see.

For the honest numbers we use the diabetes regression set that ships with scikit-learn: 442 patients, ten baseline measurements each (age, sex, BMI, blood pressure, and six blood serum values, all mean-centered and scaled), and a continuous target measuring disease progression a year later. It's a small, noisy regression problem where no model does spectacularly — which makes it a good place to see regularization actually change the answer. An 80/20 split gives 353 patients to train on and 89 held out, seeded so every run is identical.

Build it, one function at a time

The build follows the math in the order the math was derived. The loss enters only through gradients and hessians, so that's first — and it's two lines per loss:

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

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


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

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

For squared error the gradient is the residual and the hessian is a constant one; for logistic it's p - y and p(1-p). We do regression in this chapter, but the logistic version is here to make the point concrete: swap this one function and the entire rest of the file becomes a classifier. Nothing else moves.

Next the leaf weight — the formula everything else serves:

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

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

One line: $-G/(H+\lambda)$. The structure score that ranks splits is the positive piece of the same expression, $G^2/(H+\lambda)$:

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

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

With the score in hand, split gain is the children's score minus the parent's, halved, minus $\gamma$ — the formula from the math section, verbatim:

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

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

That - gamma at the end is doing quiet, important work. A split that barely improves the structure score returns negative gain and gets rejected, so the tree stops growing exactly where the extra leaf stops paying for itself. Now the split search. It's the same shape as the CART search from week 3 — every feature, every threshold — but the score is XGBoost's gain, and we sweep each feature in a single sorted pass, accumulating $G_L$ and $H_L$ as points cross to the left child:

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

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

The sorted sweep is the trick that makes this cheap: sort a feature once, then moving the threshold from one value to the next just adds one point's gradient and hessian to the running left sums. Every candidate cut on a feature is scored in one linear pass, not a fresh recount. The search returns nothing when no split has positive gain — that is the $\gamma$-pruned stopping rule doing its job. Growing the tree is that search, applied recursively, with each leaf's value set by the weight formula:

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

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

A leaf carries $w^{*} = -G/(H+\lambda)$ computed from the gradient and hessian sums of the rows that reached it. Notice what the tree predicts: not a label, but a correction — a small number to be added to the running prediction. Routing a sample to its leaf is the plain tree walk:

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


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

Finally the boosting loop ties it together. Start every prediction at the base score (the mean of y, the best constant under squared error), and each round compute the gradient and hessian at the current prediction, grow one regularized tree on them, and add its shrunk output:

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

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

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

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

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

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

staged_predict is the only piece that isn't strictly necessary to fit — it unrolls the ensemble one tree at a time, so we can replay the prediction after every round. That's what the animation reads.

Watch it work

Here's the payoff. This is the real booster from the file above, fit to the 1-D toy — thirty rounds, depth-2 trees, learning rate 0.3, $\lambda = 1$. The top panel is the prediction: a staircase (orange) laid over the training points (cyan) and the held-out validation points (amber diamonds). The bottom panel traces mean squared error on both splits as the rounds accumulate. Every frame is one real boosting round.

Press play. Round 0 is the base score — a flat line at the mean, wrong everywhere, with train and val MSE both high. Round 1 adds one depth-2 tree, and the flat line grows its first two risers where the curve bends most. From there each round adds another shallow tree that nudges the staircase closer, and the loss panel falls fast at first — the big, easy corrections come early — then flattens as the trees start fighting over noise. By round 30 the staircase hugs the curve, training MSE is down to 0.038, and validation MSE has settled around 0.110. Reset and run it again; it's the same every time.

Watch the gap between the two loss lines. Early on they fall together — the model is learning real structure that helps on both. Late, the train line keeps inching down while the val line stops. That parting is the model beginning to fit noise, and it's exactly the thing $\lambda$, $\gamma$, and the small learning rate are there to slow down. To see that lever directly we go to the real data. This is test RMSE on the held-out diabetes patients, round by round, for three values of $\lambda$ at a fixed depth of 4:

This is the money shot for regularization. The orange line is $\lambda = 0$ — no L2 penalty, the raw Newton step. It drops fastest and hits its best test RMSE of 53.07 by round 6, then turns and climbs: past that point every tree is memorizing training noise and generalization gets worse. The cyan line, $\lambda = 1$, is slower to descend but goes lower, bottoming out at 52.11 around round 11 before it too begins to drift up. The purple line, $\lambda = 30$, is over-regularized — it's so cautious it never overfits, but it never gets low either, floored around 55.44. That's the whole tradeoff in one chart: too little regularization and the model overfits fast, too much and it underfits, and the useful setting is a moderate penalty that reaches a lower minimum and holds it longer. Plain gradient boosting only has the learning rate and tree count to play with; XGBoost hands you $\lambda$ and $\gamma$ as two more.

The full implementation

The whole booster, no library — the loss derivatives, the three regularized formulas, the split search, the tree, and the boosting loop. This is the file the animation ran:

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

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

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

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

import numpy as np
import pandas as pd


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

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


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

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


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

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


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

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


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

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


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

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


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

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


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


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


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

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

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

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

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

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


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


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


def load_diabetes(path="../data/diabetes-regression.csv"):
    """The sklearn diabetes regression set: 442 patients, 10 features, a
    continuous disease-progression target."""
    df = pd.read_csv(path)
    features = [c for c in df.columns if c != "target"]
    X = df[features].to_numpy(float)
    y = df["target"].to_numpy(float)
    return X, y, features

The library version

Nobody ships a hand-rolled booster, and you shouldn't either once you've built one. The real xgboost package is the same algorithm we just wrote — second-order boosting, the $-G/(H+\lambda)$ leaf weight, the $\gamma$-penalized gain — with a decade of engineering under it. To make the comparison honest we hold the knobs equal and ask it to use the exact split finder:

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

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

The one setting that matters for matching is tree_method="exact". XGBoost's fast default, "hist", buckets each feature into bins and scans the bins instead of the raw values — a huge speedup on big data, but it would find slightly different splits than our full scan, and the numbers wouldn't line up. With "exact" it does the same greedy search we do. We also turn off the things our version doesn't have — column and row subsampling, the L1 penalty — so any remaining difference is implementation detail, not a different model. For contrast we bring in sklearn's plain gradient boosting, the first-order machine from the previous chapter, with no hessian and no explicit leaf regularization:

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

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

Scratch versus library

Three boosters, same 200 trees at depth 3, same 0.3 learning rate, fit on the 353 training patients and scored on the 89 held out. Our from-scratch second-order booster, the real xgboost, and sklearn's plain gradient boosting. Test RMSE, lower is better:

Two things to read here. First, our booster and the real xgboost land on top of each other — 59.50 test RMSE for ours, 60.01 for the package, a difference of less than one percent, and their predictions on the held-out set correlate at 0.9978. That's the result you want from a from-scratch build: the same model, confirmed by the reference implementation. The small gap that remains is the parts we deliberately left out — the package's own default base score handling, floating-point order of summation, tie breaking in the split search — not a disagreement about the algorithm. The exact figures live in results.json, regenerated whenever the code changes.

Second, and this is the point of the chapter: both second-order boosters beat plain gradient boosting. sklearn's GBM, same tree count and depth and learning rate but no hessian and no leaf regularization, comes in at 64.29 RMSE — an R-squared of 0.168 against 0.288 for ours and 0.275 for xgboost. On this small, noisy set the regularized Newton objective is worth about twelve points of RMSE over the first-order machine, for free, with no extra tuning. That's the twelve points that won all those competitions.

Takeaways

XGBoost is gradient boosting with the two upgrades that turned a strong model into the default one. The second-order objective uses the curvature of the loss, not just its slope, so each tree takes a Newton step instead of a gradient step and lands closer on the first try. The regularized tree score bakes an L2 leaf penalty and a per-leaf complexity charge straight into the split math, so the ensemble prunes itself as it grows rather than memorizing the training set. Add shrinkage and you have three independent brakes on overfitting, and the diabetes numbers show them earning their keep — twelve points of RMSE over the same-shaped plain GBM.

Reach for it when the data is a table and you want the best number with the least ceremony. It handles mixed types and missing values, needs no scaling, trains fast, and lands near the top of tabular leaderboards before you've tuned anything. The knobs that actually move the needle, in the order I touch them: the learning rate and tree count together (lower rate, more trees, until the validation curve stops improving), then max_depth for how much interaction each tree can capture, then $\lambda$ and $\gamma$ to rein in overfitting once depth is set. The regularization chart is the mental model — there's a moderate setting that reaches a lower minimum and holds it, and your job is to find it with a validation set, not to crank every penalty to the sky.

When is it overkill? When a random forest already does the job — if you don't need the last few percent and you'd rather not babysit a learning-rate schedule, bagged trees give you most of the accuracy with almost none of the tuning, and they're far harder to overfit. And a single deep tree is still the right answer when you need something a human can read and argue with, which a boosted ensemble never will be. XGBoost is what you graduate to when the tabular accuracy is the thing that matters most and you're willing to spend a tuning budget to get it. On the kind of data most of us actually work with, that's more often than not.