Capítulo 23 de 37 · intermedio

Gradient boosting

What this chapter covers

The random forest built trees in parallel and averaged them: every tree an independent guess, the crowd wiser than any member. Boosting does the opposite. It builds trees in sequence, and each new tree exists for one reason — to fix the mistakes the trees before it are still making. AdaBoost, the chapter before this, did that by re-weighting the examples the ensemble kept getting wrong. Gradient boosting does something cleaner and more general: it looks at the errors as a gradient and fits the next tree to point straight downhill.

This is the algorithm that quietly wins tabular-data competitions. We'll build it from scratch in NumPy for the cleanest case — regression with squared-error loss — where the whole method collapses to something almost embarrassingly simple: guess the mean, look at what's left over, fit a shallow tree to the leftovers, add a shrunken slice of it, and repeat. The leftovers have a name, the residuals, and the one idea that makes this a gradient method is that for squared error the residual is the negative gradient of the loss. Fit the residual and you've done gradient descent, except the thing you're stepping isn't a vector of weights — it's the prediction function itself.

The centerpiece is an animation: one boosting round per frame on a 1-D curve, the running prediction creeping up from a flat line into a staircase that hugs a sine wave while the residuals collapse toward zero underneath it. Then we hand the same job to scikit-learn's GradientBoostingRegressor and check the numbers line up, and we watch the one knob that defines the method — the learning rate — trade speed against overfitting.

A bit of history

Boosting started as a theoretical question. In 1990 Robert Schapire proved that a "weak" learner — one that only has to beat random guessing by a hair — could be boosted into an arbitrarily accurate "strong" one, which was a surprising result about what was possible in principle. Yoav Freund and Schapire turned it into a practical algorithm with AdaBoost in 1995, and for a few years AdaBoost was a bit of a mystery: it worked far better than the theory said it should, and nobody was quite sure why.

The clarifying move came from statistics. In 1998 Leo Breiman noticed that AdaBoost was doing a kind of gradient descent, and in 2000 Jerome Friedman, Trevor Hastie, and Robert Tibshirani wrote the paper that reframed the whole thing: AdaBoost is fitting an additive model by a stagewise optimization of an exponential loss. Once you see it that way, the specific loss stops being special. Friedman took the idea to its conclusion in "Greedy Function Approximation: A Gradient Boosting Machine" — a 1999 tech report published in the Annals of Statistics in 2001 — and gave us the general recipe: pick any differentiable loss, compute its negative gradient at the current predictions, fit a regression tree to that gradient, and take a step. AdaBoost falls out as the special case for one particular loss. Squared error, the case we build here, falls out as the simplest one, where the gradient is just the residual. That single paper is the direct ancestor of every gradient-boosting library — XGBoost, LightGBM, CatBoost — that now sits at the top of the tabular-data leaderboard.

The intuition

Here's the trick, in plain terms, before any calculus. Suppose your current model predicts some number for every point, and it's wrong. The errors — how much you missed each point by — are themselves a dataset: an input x and a "target" equal to the miss. So train a model to predict the misses. If it can predict even a little of the pattern in your errors, you can add its predictions to your model and shrink the errors. Now you have a slightly better model, with slightly smaller errors, and you do it again. Each round chases the leftover.

That leftover is the residual, y - F(x): the truth minus what the model currently says. The very first model is the laziest one possible — predict the mean of the target for everybody, which is the single best constant guess. Its residuals are big. Fit a shallow tree to those residuals, add a fraction of it to the prediction, and the residuals shrink. Fit another tree to the new, smaller residuals. The model is the running sum of that starting constant and every shrunken tree.

Why shallow trees, and why a fraction? Because each tree is supposed to be a small correction, not a full answer. A deep tree fit to the residuals in one shot would wipe them out on the training set and memorize the noise along with the signal. A stubby tree — two or three levels — can only capture a coarse piece of the leftover pattern, which is exactly what you want when you're going to stack fifty or a hundred of them. The fraction is the learning rate, and it's the same idea as the step size in gradient descent: take small steps and you approach carefully; take big ones and you overshoot. The whole method is gradient descent, done in the space of prediction functions instead of the space of weights.

The math

Boosting builds an additive model — a running prediction that starts at a constant and gets one tree added at a time. After $M$ rounds it is

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

where $F_0$ is the initial constant, each $h_m$ is a shallow regression tree, and $\nu$ is the learning rate — a small positive number, typically 0.1, that shrinks every tree's contribution. The symbol $x$ is one input row and $y$ its target; there are $n$ training rows.

We fit the model one term at a time to drive down a loss. For regression the loss is squared error, written here with a one-half that makes the derivative clean:

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

Where do we start? With the constant $F_0$ that minimizes the total loss on its own. For squared error that constant is the mean of the target:

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

Now the key step. We want each new tree to move the prediction in the direction that reduces the loss fastest — the negative gradient of the loss with respect to the current prediction, evaluated at each training point. Differentiate the squared-error loss with respect to $F$ :

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

Flip the sign and the negative gradient at row $i$ , evaluated at the current model $F_{m-1}$ , is

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

That is the residual, exactly. This is the whole reason squared error is the clean case to learn on: the negative gradient you're supposed to chase is just "truth minus current guess," no calculus needed to compute it in practice. We fit the new tree $h_m$ to those residuals by ordinary least squares — the regression tree's native objective — and then take a shrunken step:

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

Repeat for $M$ rounds. Swap squared error for a different differentiable loss — absolute error, Huber, the logistic loss for classification — and only one line changes: the formula for $r_i$ . Everything else, the tree fit and the shrunken step, stays put. That generality is the point of the word "gradient" in the name.

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

Gradient boosting is, for my money, the strongest thing you can throw at a tabular dataset without a neural network, and usually with one too. It inherits the good manners of trees — no feature scaling, mixed numeric and categorical inputs, nonlinear boundaries and interactions found on their own — and then it does something the random forest can't: because each tree corrects the last one's errors instead of just averaging an independent guess, it drives the bias down round after round and reaches accuracy a forest of the same trees can't touch. Give it a differentiable loss and it optimizes that loss directly, so you can fit it to squared error, absolute error, quantiles, or a classification objective by changing one formula. When there's a real prize on a tabular problem, the winning model is almost always a gradient-boosted ensemble.

The costs are real and they're the mirror image of the forest's. It's sequential, so it doesn't parallelize the way bagging does — tree m needs tree m-1's predictions to compute its residuals. It has more knobs that interact: the learning rate and the number of trees trade off against each other, and get them wrong and you either underfit or memorize. And unlike a forest, more trees is not always safer — past a point the ensemble starts fitting noise and the test error turns around and climbs, which we'll watch happen. It's a model you tune, not one you fit and forget. The payoff is worth the babysitting more often than not, but you do have to babysit.

The data

Two datasets, two jobs. The concept animation uses a 1-D curve I can draw every point of: 120 points with x spread over [0, 10] and the target a sine wave on a gentle upward slope, y = sin(x) + 0.25x, plus Gaussian noise. One input, one output, a shape no single shallow tree can capture — perfect for watching a staircase of trees build up toward it.

For the honest accuracy numbers we switch to a real regression set: scikit-learn's diabetes data, 442 patients with ten standardized physiological measurements each, and a target that's a quantitative measure of disease progression one year after baseline. It's a small, noisy, genuinely hard regression — nobody gets a great R² on it — which makes it the right place to see what boosting actually buys over a single tree and where the learning rate starts to hurt. Same style of split as the earlier chapters: 309 patients to train on, 133 held out to test.

Build it, one function at a time

Two pieces stack here. First the weak learner — a shallow regression tree — then the boosting loop that grows a sequence of them. The tree is week 3's CART with one substitution: it splits to reduce squared error instead of Gini impurity, and a leaf predicts a number instead of a class.

Purity for a regression node is just how spread out the target is inside it, which is the sum of squared errors around the node's own mean:

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

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

This is the exact analogue of Gini. A node full of identical values has sse zero; the more the values scatter, the bigger it gets, and a split earns its keep by leaving its two children with less total sse than the parent had. A leaf, when the splitting stops, predicts the one number that minimizes its own sse — the mean:

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

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

The best-split search is week 3's, line for line, with sse swapped in for the Gini score. Loop over every feature, every candidate threshold, and keep the split that most reduces the total squared error of the two halves:

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

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

If you read the decision-trees chapter this will look familiar to the character — the only change is what "purity" means. Growing the tree is the same recursion, too, stopped early by a shallow depth cap because in boosting we want a weak learner, not a perfect one:

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

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

max_depth of two or three is the whole discipline here. Each tree is meant to capture a coarse slice of the leftover pattern and nothing more. Prediction is a walk from the root to a leaf, returning the leaf's mean:

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


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

That's the weak learner. Now the part that makes it boosting. Start the prediction at the mean, and then round after round: compute the residuals, fit a tree to them, and add a shrunken step of that tree to the running prediction:

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

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

Read the three lines in the loop against the math. residual = y - F is the negative gradient of squared-error loss — the thing we just derived. build_tree on that residual is the least-squares fit of $h_m$ . F = F + learning_rate * tree_predict(...) is the shrunken step $F_m = F_{m-1} + \nu h_m$ . Three lines, and they are the entire algorithm; everything above them is the tree they call. The record flag just logs the state of the prediction after each round so the animation has something to replay.

Prediction with the finished model is the additive formula spelled out — the initial constant plus every shrunken tree:

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

And two metrics to score it, mean squared error and the $R^2$ that tells you how much better than guessing the mean you're doing:

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


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

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

Watch it work

Here's the payoff, and it's the clearest picture of boosting I know how to draw. This is the real fit_gradient_boost running on the 1-D curve, one round per frame, at a deliberately small learning rate of 0.1 so the staircase builds slowly enough to watch. The top panel shows the current prediction F(x) — an orange step function, because a sum of shallow trees is piecewise constant — creeping up toward the true curve (the dashed gray line) through the cyan data. The bottom panel is the residuals, y - F(x), the leftover each round is trying to kill. The caption names the round, the learning rate, and the training MSE.

Press play. Round zero is the lazy model: F(x) is a flat line at the mean, 1.35, and the residuals in the bottom panel are a wide cloud — the model knows nothing yet, so every point is "wrong" by its full distance from the average. Then the trees start landing. The first few carve the biggest, coarsest steps — they grab the overall up-slope and the largest hump of the sine — and you can see the residual cloud contract hardest in those first rounds. As the rounds pile on, the steps get finer and more numerous, the staircase bends around the wave, and the residuals squeeze toward the zero line. Training MSE falls from 1.007 at round zero to 0.069 by round fifty, monotonically, every single round smaller than the last.

Watch two things in particular. First, no single frame's tree does much — each one nudges. That's the shrinkage: at a learning rate of 0.1 every tree contributes a tenth of what it found, so the fit approaches carefully instead of snapping into place. Second, look at where it is around the middle of the run, near round twenty-five. The staircase already tracks the curve well and the residuals are roughly the size of the noise I baked into the data — that's about as good as any model can honestly do here. The rounds after that are refining, and toward the very end they start chasing wiggles that are noise, not signal. That creeping overfit is invisible on this clean 1-D toy, but on real data it's exactly what turns the test error around, which is the next thing to look at.

The learning rate and the number of trees are the two knobs that define the method, and they pull against each other. Here's what they do on the real diabetes data: test MSE as trees are added, for three learning rates.

This is the picture to sit with. The big learning rate, 0.5 (amber), drops fastest — it reaches its best test MSE of about 3244 in only 18 trees — and then turns and climbs steeply, overfitting hard, all the way past 5300 by 300 trees, which is no better than predicting the mean. It sprinted to a decent answer and then ran off a cliff. The small learning rate, 0.05 (cyan), is the tortoise: it descends slowly, takes about 105 trees to reach its best of roughly 3107, and overfits so gently that even at 300 trees it's only crept back to 3241. The middle rate, 0.1 (purple), bottoms out around 3105 near 40 trees, then overfits at a moderate pace. The lesson every gradient-boosting practitioner learns: a smaller learning rate reaches a lower floor but needs more trees to get there, and it's more forgiving if you overshoot the tree count. In practice you set the learning rate small, set the tree count generously, and let a validation set or early stopping tell you where to quit. You never leave the learning rate at 0.5.

The full implementation

The whole thing, no library, top to bottom — the shallow regression tree and the boosting loop that stacks it. This is the file the animation ran:

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

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

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

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

import numpy as np
import pandas as pd


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

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


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

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


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

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


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

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


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


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


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

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


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


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


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

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


def load_diabetes_csv(path="../data/diabetes.csv"):
    """Diabetes regression set: 442 patients, 10 standardized features, a
    disease-progression score one year on. Source: sklearn's load_diabetes."""
    df = pd.read_csv(path)
    features = [c for c in df.columns if c != "target"]
    X = df[features].to_numpy(float)
    y = df["target"].to_numpy(float)
    return X, y, features

The library version

Nobody hand-rolls a booster in production. scikit-learn's GradientBoostingRegressor is the same algorithm — start at the mean, fit shallow trees to the residuals, add each one back with a learning-rate step — written in C and wired for the knobs that matter:

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

loss="squared_error" and the mean initialization are the defaults, so this is a faithful twin of our fit_gradient_boost: same objective, same residual-fitting loop, same shrinkage. The differences are refinements, not a different idea. sklearn splits on friedman_mse by default — a small improvement on the plain variance-reduction split we use, weighting the split score by the child sizes — and it adds the machinery real problems need: subsample for stochastic gradient boosting (fit each tree on a random fraction of the rows, which de-correlates them a little like a forest does), n_iter_no_change for early stopping, and losses beyond squared error. Same engine, more instrumentation.

Scratch versus library

Fit both at scikit-learn's own defaults — 100 trees, learning rate 0.1, depth-3 trees — on the 309-patient training half of diabetes, score on the 133 held out. And to make the case for boosting itself, put a single depth-3 regression tree, one lonely weak learner, next to them:

The two boosters land together: our from-scratch ensemble scores 0.3766 test R², sklearn's scores 0.3916, a gap of fifteen thousandths that comes down entirely to sklearn's friedman_mse splits choosing slightly different thresholds than our plain squared-error ones. Same algorithm, essentially the same fit, which is the result you want — our hand-written residual loop and sklearn's C implementation agree. In error terms it's a test MSE of 3332.69 for ours against 3252.65 for sklearn's, both a long way below the mean baseline's 5349.11.

Now the number that argues for boosting. A single depth-3 tree — one weak learner, the thing we boosted — scores 0.2564 R² on this test set. Stacking a hundred of them with residual fitting took that to 0.3766, better than half again as much explained variance, from the same shallow trees and no new kind of model. That's the whole pitch: a depth-3 tree is a mediocre model on its own, and a hundred of them correcting each other is a strong one. One honest caveat, the same one that haunts the diabetes set everywhere: an R² of 0.38 is not a triumph, it's just what there is to get on hard, noisy medical data. The point isn't the altitude, it's the climb from 0.26 to 0.38 that boosting bought over the single tree, and the fact that our scratch code made that climb neck and neck with the library.

Takeaways

Gradient boosting is AdaBoost's idea, generalized and made honest. Where AdaBoost re-weighted the hard examples, gradient boosting names what "hard" means precisely — the negative gradient of a loss — and fits the next tree straight at it. For squared error that gradient is the residual, so the whole method reduces to the plainest loop imaginable: guess the mean, fit a shallow tree to what's left over, add a shrunken slice, repeat. Every line of the from-scratch loop maps to one line of the math, and the word "gradient" earns its place because swapping the loss changes exactly one formula and nothing else.

Reach for it when you're serious about a tabular problem. A random forest is the baseline you fit first to see if there's signal; a gradient-boosted ensemble is what you fit when you want to win, because building trees that correct each other drives the error lower than averaging independent ones ever can. The price is the tuning: the learning rate and the tree count are coupled, smaller rates need more trees but reach a lower floor and forgive you more, and unlike a forest it will overfit if you keep adding trees past the point the validation error turns around. Set the learning rate small, the tree count generous, and let early stopping find the wall. Never ship it at learning rate 0.5.

What we built is the 2001 gradient-boosting machine in its simplest clothes — plain squared-error splits, one tree at a time, from scratch. What ships today is that same skeleton with better splits, second-order gradients, regularization on the leaf values, column and row subsampling, and an engineering effort poured into making it fast on millions of rows. That's XGBoost, and it's the next chapter — not a new idea so much as this idea, sharpened into the model that sits at the top of more tabular leaderboards than anything else.