From one tree to XGBoost: gradient boosting from scratch

The previous article built a decision tree in pure Python and then spent its last section explaining why you should not ship one.

Two things were wrong with it. It overfits: pointed at the breast-cancer dataset with nothing holding it back, it grew to depth 7, scored a perfect 1.000 on the rows it was trained on, and lost on the rows it had not seen to a tree one third its depth. And it is unstable: on the full Cleveland heart study, ten runs over identical patients produced four distinct models, because tied candidate questions were resolved by the order a Python set happened to iterate.

Both problems have the same root. A tree’s structure is a hostage to the particular sample it was grown from, and the cure is not a better tree. It is more trees. There are two ways to combine them:

• Grow many trees independently, each on a different bootstrap sample of the rows and a random subset of the columns, and average their answers. That is a random forest.
• Grow them in sequence, each one fitted to what the ones before it got wrong, and add them up. That is boosting.

This article builds the second, in about forty lines on top of the tree we already have. Then it takes the finished loop apart to show that it was gradient descent all along — performed not on a vector of parameters but on a function — and that is the door XGBoost walks through. The two ensembles are compared properly, on measured numbers, once we have something to compare with.

The regression stump we need

The previous article left a promissory note in its opening section:

point the same code at a numeric target, replace the impurity measure with the variance of that target and the leaf’s label counts with their average, and you have a regression model that predicts a number — nothing else changes.

We are going to need that, so let’s collect it — in compressed form here, since the regression-tree article works through it properly, including the one limitation boosting inherits and never escapes. The change really is two functions. Where the classifier counted labels and measured how mixed the counts were, the regressor averages numbers and measures how spread out they are:

def mean(rows):
    return sum(row[-1] for row in rows) / float(len(rows))

def variance(rows):
    targets = [row[-1] for row in rows]
    m = sum(targets) / len(targets)
    return sum((t - m) ** 2 for t in targets) / len(targets)

variance replaces gini inside info_gain, mean replaces class_counts inside the leaf, and every other line — Question, partition, find_best_split, the recursion in build_tree — is untouched. The greedy search still tries every feature paired with every value it takes and keeps whichever split removes the most impurity; only the meaning of “impurity” moved from mixed labels to spread-out numbers.

Three facts about that tree are all the rest of the article needs: leaves predict the mean of the rows that reach them, splits are chosen to reduce squared error, and — the one thing we still have to add — the tree is capped, so each round can only make a limited correction.

The cap is one line in the recursion:

def build_tree(rows, max_depth, depth=0):
    """CART regressor, capped at max_depth. depth=1 gives a stump."""
    if depth == max_depth:
        return RegressionLeaf(rows)
    gain, question = find_best_split(rows)
    if gain == 0:
        return RegressionLeaf(rows)
    true_rows, false_rows = partition(rows, question)
    return Decision_Node(
        question,
        build_tree(true_rows, max_depth, depth + 1),
        build_tree(false_rows, max_depth, depth + 1),
    )

At max_depth=1 the tree asks exactly one question and stops — one question, two leaves, two possible answers for the entire dataset. That is a stump, and it is the smallest useful weak learner.

The instinct is that this must be a compromise, and that we would use an unconstrained learner if we could afford one. It is not. Boosting’s next tree is trained on what the previous ones left behind, so a learner powerful enough to fit the training data outright leaves nothing behind for it to work on, and the ensemble collapses toward a single overfitted model. Each member has to be limited enough that it corrects only part of the error and passes the rest along.

“Limited” is relative to the problem, not an absolute setting. A stump is the extreme case and it is what we will use here, because on five rows you can check it by hand. On real data the sweet spot is usually deeper — we measure it later, and depth 3 beats stumps on the breast-cancer set. What matters is that each round’s correction is constrained, not that it is tiny.

The data

Five players from the Hitters study, two numeric predictors and a numeric target:

years is seasons played in the major leagues, hits is hits made in the previous season, and salary is the 1987 annual salary in thousands of dollars. Five rows again, so every number below can be checked by hand.

Rows 1 and 2 are worth noting now because they set a floor we will hit later. BillyJo Robidoux and Jack Howell have identical features — two years, forty-one hits — and different salaries, 67.5 and 95.0. No model reading only these two columns can tell them apart, exactly like the two heart patients who produced the 50/50 leaf in the previous article.

Build boosting by hand

Round zero needs a prediction before any tree exists. Take the simplest constant available — the mean salary:

It is a terrible model. It predicts 420 for a player earning 67.5 and 420 for one earning 787.5, and its total squared error is 431,037.5. But it is a starting point, and boosting only ever needs a starting point.

Now ask what each player’s prediction is missing — the residual, the actual salary minus what we currently predict:

Read that last column as a to-do list. It says: this model is 352.5 too high for Robidoux and 367.5 too low for Moseby. Fixing it means adding −352.5 to the first prediction and +367.5 to the last.

And here is the move the whole method rests on. We fit a tree to that column. Same two features, years and hits — a tree always splits on features, and those never change — but the target is no longer salary. The target is what salary is still missing:

residuals = [row[-1] - p for row, p in zip(rows, predictions)]
# same features, new target: the residual
residual_rows = [row[:-1] + [r] for row, r in zip(rows, residuals)]
tree = build_tree(residual_rows, max_depth)

Run the stump-builder on those five residual rows and it comes back with:

Is hits >= 77?
--> True:  Predict +225.83
--> False: Predict -338.75

Which is the greedy split we already know how to compute. Players with 77 or more hits are the three whose residuals are positive — +60, +250, +367.5 — and the leaf holds their mean, $677.5 / 3 = 225.83$. The other two have residuals of −352.5 and −325, and their leaf holds $-677.5 / 2 = -338.75$.

(That split is not as inevitable as it looks: Is years >= 3? scores exactly the same, because it cuts these five players into the same two groups, and the >= comparison in find_best_split breaks the tie in favour of whichever column is scanned last. The regression-tree article works through that table in full.)

Add that correction to every prediction and round one is done:

Squared error falls from 431,037 to 48,532 — one question, and 89% of the error is gone.

Look at what happened to the two collided players. Both sit at exactly 81.25, which is the mean of 67.5 and 95.0 — the best any model can do for a pair it cannot distinguish. One stump found the floor for them on the first try. Their residuals are now ±13.75 and will never usefully shrink again.

Then you do it again. The residual column is recomputed against $F_1$, a second stump is fitted to that, and so on:

Drag through the rounds and watch two things. The residual column drains — the numbers get smaller and change sign as the model overshoots and gets corrected. And the questions change: round one asked hits >= 77, round two asks hits >= 149, round three asks hits >= 130. Each stump is answering a different question because each is looking at a different target. The features are fixed; what is being predicted is not.

The other thing worth catching is that progress is not monotone per player. At round 1 the two collided players sit at their optimal 81.25; by round 2 they have been dragged to 45.83, because that round’s stump was chasing Alvin Davis’s −165.83 and they were on the wrong side of the split. Greedy correction is a global bargain, not a promise to each row.

Now that the operation has a shape, it can be named. Every round appends one more term, so the finished model is a sum:

Each $f_m$ is one weak learner, and the prediction for a row is what you get by dropping it down all $M$ of them and adding up the answers. Nothing is ever revised — trees are only ever appended.

The whole algorithm

That is the entire method, so we can write it down. Assembled, it is short enough to read in one sitting:

def fit_gbm(rows, n_trees, learning_rate, max_depth=1):
    """Returns (F0, [tree, tree, ...]). Squared error, so residuals = -gradient."""
    F0 = mean(rows)
    predictions = [F0] * len(rows)
    trees = []

    for _ in range(n_trees):
        residuals = [row[-1] - p for row, p in zip(rows, predictions)]
        # same features, new target: the residual
        residual_rows = [row[:-1] + [r] for row, r in zip(rows, residuals)]
        tree = build_tree(residual_rows, max_depth)
        step = [predict(row, tree) for row in residual_rows]
        predictions = [p + learning_rate * s for p, s in zip(predictions, step)]
        trees.append(tree)

    return F0, trees


def predict_gbm(row, F0, trees, learning_rate):
    total = F0
    for tree in trees:
        total += learning_rate * predict(row, tree)
    return total

Eleven lines of loop on top of the previous article’s tree, and predict_gbm is the additive model made literal — start at the base constant, walk every tree, add up what they say. The one piece that has not been explained yet is learning_rate; set it to 1.0 and this is exactly the arithmetic worked through above, which is what the next section is about.

We now have a complete gradient boosting regressor for squared error. Everything that follows explains why this same loop works for losses that are not squared error at all.

Worth pausing first on what training looked like. The previous article observed that a tree “never got gradually better; it got gradually built” — one greedy pass, no epochs, no convergence. Boosting puts the iteration back, but one level up: each individual tree is still built in one greedy pass and never revisited, while the ensemble genuinely converges, round after round, toward lower loss.

Controlling each correction

Three knobs decide how much the model changes per round and how many rounds there are. They are easiest to understand together, because all three constrain the same thing.

The learning rate

The version above adds each tree at full strength. Real implementations add a fraction:

where $\eta$, the learning rate (or shrinkage), is typically 0.1 or 0.01. Every tree still gets built the same way, but only a fraction of its correction is applied.

This looks like pure waste — you fit a tree and then deliberately ignore 90% of it. What it buys is that no single tree gets to dictate the answer. A stump fitted to a residual column is fitted to noise as well as signal, and at full strength its mistakes go straight into the model, where later rounds have to spend themselves undoing it. At a tenth strength, any one tree’s contribution is a suggestion the next hundred trees can outvote.

The cost is rounds. Fitting our five players down to their irreducible floor takes:

So $\eta$ and the number of trees $M$ are tightly coupled: lower the rate and you need more trees to travel the same distance. It is tempting to treat them as a single budget traded against each other, but they are not interchangeable. Shrinkage is a regularizer, not just a step size — a model that reaches a given training loss in 500 small steps generalizes differently from one that got there in 50 large ones, and usually better. The measured sweep later in this article shows the two settings landing on genuinely different test scores rather than the same score at different speeds.

The widget below is the same loop on a smooth one-dimensional curve, which makes the mechanism visible in a way five rows cannot. The grey whiskers are the residuals — the exact quantity the next stump will be trained on — and the blue staircase is $F_m$:

At learning rate 1.0 the staircase lunges at the points and is essentially converged by round 20. At 0.1 it creeps up from the flat line, and after 60 rounds it is still visibly short.

Notice also what the staircase is: every riser is a threshold some stump chose, and the height of each tread is the sum of every stump’s answer for that strip of the axis. A hundred stumps produce a hundred-odd risers, which is how a model made of one-question trees ends up approximating a smooth curve. It never stops being piecewise-constant — it just uses very small pieces.

Depth is the interaction order

Depth in a boosted model does not mean what it means in a single tree. There, depth was raw capacity — how finely the tree could carve up the data before it started memorising rows. Here each tree is one small correction, and depth controls how many features a single correction is allowed to combine:

• max_depth=1 — every tree asks about one feature. The model is a sum of one-feature functions, an additive model with no interactions at all. It cannot express “high cholesterol and low heart rate”.
• max_depth=2 — each tree can chain two questions, so the ensemble captures pairwise interactions.
• max_depth=d — up to $d$-way interactions.

So depth is a statement about the structure of your problem rather than a dial to be maximised. Two to eight is the usual range, and this is why the stumps we have been using are a teaching choice rather than the default you would reach for.

The number of rounds

The third knob is $M$ itself, and unlike the other two it has no safe default. Every round adds capacity, so a long enough run will eventually start fitting noise. This is the one place a boosted model behaves worse than a forest, and it is measured on real data below.

Why residuals are gradients

Everything so far was arithmetic on residuals, and the word “gradient” has not appeared. Here is where it does.

Write the loss for a single observation as squared error, with the customary one-half in front:

Now differentiate it with respect to the prediction — not with respect to a parameter, which is what we are used to differentiating, but with respect to the number the model currently outputs:

The residual is the negative gradient of the loss with respect to the prediction. We have been computing gradients since round zero without calling them that.

That reframing costs nothing and changes everything, because nothing in the loop required the loss to be squared error. Look back at the code: we computed a column of numbers, fitted a tree to it, and added the result. Only the first of those three steps mentioned the loss. So replace it:

These are the pseudo-residuals. For squared error they are the ordinary residuals we have been using. For any other differentiable loss they are some other column of numbers — and every remaining line of the algorithm runs unchanged, because a tree fitted to a column of numbers does not care where the numbers came from.

This is the flexibility the method is famous for. One loop serves regression, classification, ranking, Poisson counts, and survival times; what changes between them is one derivative. The only real requirement is the one that gives the method its name: the loss must be differentiable with respect to the prediction, and cheaply so, because you evaluate it once per training row per round.

Why fit a tree to the gradient at all?

The gradient column tells you exactly how to improve the prediction for each of the $n$ training rows. Why not apply it directly and skip the tree?

Because a column of $n$ numbers is not a model. It says what to do at the training rows and says nothing at all about any other input. Apply it directly and you have memorised the training set perfectly and learned nothing that transfers — the pathological case the previous article’s fully grown tree walked into.

Fitting a tree to the gradient is what makes the step generalize. The tree is forced to explain the gradient using the features, and because it is capped it can only explain it in coarse terms: everything with hits >= 77 needs about +226. That is a statement about the feature space, not about five specific rows, so it applies to a player nobody has seen. The weak learner is a deliberately low-resolution summary of where the loss wants to move, and the low resolution is the point.

Gradient descent in function space

We now have everything needed to say what the method is.

Training a neural network means holding the architecture fixed and repeatedly nudging a vector of parameters against the gradient:

Boosting does the same thing one level up. The object being optimized is not a parameter vector — it is the prediction function itself:

Each round computes $-\nabla_F L$, and each round adds a scaled copy of it to the model. Every tree in the ensemble is one gradient descent step, and the finished model is the trajectory of that descent, stored as the sum of every step it took.

The one place the analogy needs care is exactly where the tree lives. A gradient in function space is not a vector you can store — it is defined at your $n$ training rows and undefined everywhere else. The tree is what converts those $n$ numbers into an actual function that can be evaluated at unseen inputs, and the price of that conversion is that the step is approximate. You do not move exactly along the gradient; you move along the closest thing a capped tree can express.

That is also the honest answer to why small trees and small learning rates work well together. Each step is a crude approximation of the true descent direction, so you take a small one, and then you take twenty thousand more.

Changing the loss: classification

The previous article was about classification, so let’s prove the generalization rather than just claiming it. We predict whether a patient has heart disease, on the same five patients it used:

The substitution is completely controlled — four components change and nothing else does:

That last row is the one to hold on to. Even for a classification problem, boosting fits regression trees, because what they are fitted to is a column of real-valued gradients.

Predictions live in log-odds. A sum of trees produces an unbounded real number, which cannot be a probability. So we let the sum be a log-odds score $F(x)$ and map it to a probability only at the end, with $p = \sigma(F) = 1/(1 + e^{-F})$. Three Yes against two No gives $F_0 = \log(3/2) = 0.4055$, which is $\sigma^{-1}(0.6)$.

The pseudo-residual is $y - p$. Differentiating $L = -[y\log p + (1-y)\log(1-p)]$ with respect to $F$ — through the sigmoid, where the algebra collapses beautifully — gives $\partial L/\partial F = p - y$. So the negative gradient is $y - p$: the label minus the currently predicted probability. It is a residual again, just measured in probability rather than in dollars.

Round one on the five patients:

Is vessels >= 1?
--> True:  Predict +1.6667
--> False: Predict -1.1111

stress       vessels    y        p   residual      F_1      p_1
normal       0          0   0.6000    -0.6000  -0.7056   0.3306
fixed        0          1   0.6000     0.4000  -0.7056   0.3306
reversable   2          1   0.6000     0.4000   2.0721   0.8882
reversable   1          1   0.6000     0.4000   2.0721   0.8882
fixed        0          0   0.6000    -0.6000  -0.7056   0.3306

The first stump asks Is vessels >= 1? — the same question the previous article’s tree chose as its root, and for the same reason: it is the split that best separates these patients. The tie that made that choice arbitrary back then is still there; it is simply less consequential now, because this stump is one of many rather than the root that determines everything below it.

How much of the new tree to add

Now look at that leaf value, because it is not what the tree computed. The two reversable patients have residuals of +0.4 each, whose mean is +0.4 — but the leaf holds +1.6667, four times that.

This is the step we have been quietly skipping. Having fitted $f_m$ to the pseudo-residuals, the classical algorithm works out how far to travel along it, by solving a one-dimensional minimization called a line search:

In words: the tree said which direction to move each prediction, and this asks how big a step to take, by keeping whichever step size leaves the total loss lowest.

Under squared error you never see it, because the answer is always exactly 1:

round 1: gamma = 382505.2083 / 382505.2083 = 1.000000000000
round 2: gamma =  25086.8056 /  25086.8056 = 1.000000000000
round 3: gamma =  14173.7558 /  14173.7558 = 1.000000000000

The tree had already done the line search for us. Its leaves hold the mean residual of the rows that reach them, and the mean is precisely the constant that minimizes squared error for those rows. Nothing is left to scale, so $\gamma = 1$ and the step is a no-op — which is why the section on the algorithm could skip it without lying.

Log-loss is the first case where that breaks. The tree’s leaves still hold means, because the tree is grown by variance reduction, which is a squared-error criterion — but the mean is no longer the value that minimizes log-loss over the rows in that leaf. So Friedman’s refinement runs the search per leaf:

For our reversable leaf that is $0.8 / 0.48 = 1.6667$, against a plain mean of $0.4$. The denominator is a sum of curvatures, and $p(1-p)$ never exceeds $0.25$ — whereas the second derivative of squared error is exactly 1. Log-loss curves far more gently than squared error, so dividing by that curvature produces a correspondingly longer stride. A plain mean would have crawled.

Run it 50 rounds at $\eta = 0.3$ and the patients settle:

stress       vessels    y           F           p
normal       0          0     -8.8674     0.0001
fixed        0          1      0.0002     0.5000
reversable   2          1      9.0602     0.9999
reversable   1          1      9.0602     0.9999
fixed        0          0      0.0002     0.5000

Patients 2 and 5 land on exactly 0.5000 — the same 50/50 the previous article’s tree reported at its mixed leaf. They share fixed/0 and disagree on the outcome, and boosting cannot separate them any more than a single tree could. Hundreds of trees, second-order steps, a differentiable loss: none of it manufactures information the columns do not contain. The two models even report the irreducible error identically, one as a leaf count and one as a converged log-odds of zero.

What happens on real data

The pure-Python version above is for reading, not for running at scale — its split search is the brute-force $O(\text{features} \times \text{values} \times \text{rows})$ scan the previous article ended by criticising. For measurements we use the libraries, on the identical split that article used: breast cancer, test_size=0.3, random_state=0, giving 398 training rows, 171 test rows, and 30 features. The single-tree rows below reproduce its numbers exactly, so everything is measured against the same yardstick.

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from xgboost import XGBClassifier

data = load_breast_cancer()
X_tr, X_te, y_tr, y_te = train_test_split(
    data.data, data.target, test_size=0.3, random_state=0
)

for name, model in [
    ("single tree, unpruned",   DecisionTreeClassifier(random_state=0)),
    ("single tree, max_depth=2", DecisionTreeClassifier(max_depth=2, random_state=0)),
    ("random forest, 100 trees", RandomForestClassifier(n_estimators=100, random_state=0)),
    ("gradient boosting",        GradientBoostingClassifier(random_state=0)),
    ("XGBoost, defaults",        XGBClassifier(random_state=0, eval_metric="logloss")),
]:
    model.fit(X_tr, y_tr)
    print(name, model.score(X_tr, y_tr), model.score(X_te, y_te))

The ordering is what we came for: the tree the previous article grew is last, tuning its depth recovers a good deal, a forest does better still, and boosting is at the top — while still scoring 1.000 on training, the same number that meant “memorised” for a single tree and means nothing of the kind here.

Two honest caveats. This is one dataset, one split, one seed, which is enough to illustrate the mechanism and not enough to establish a general ranking; on other data forests and boosted models trade places routinely. And the depth-3 row beating the stump row by 0.024 is the concrete version of the earlier point that “weak” is relative — stumps are the teaching case, not the best setting.

That comparison also fills in the contrast the introduction deferred:

The number of trees, and early stopping

That last row needs measuring rather than asserting, and on this dataset the naive version of the experiment does not show it:

Boosting plateaus; it does not decay. Breast cancer is clean, well-separated data, and there is very little noise available to overfit to.

Flip 15% of the training labels, and the asymmetry appears immediately:

Boosting peaks at 50 trees and then loses ground steadily — it is spending its later rounds fitting the 59 corrupted labels, because those rows carry the largest residuals and the loop is built to chase the largest residuals. The forest climbs and then settles. Under noise, boosting’s strength becomes its failure mode.

Which is why boosted models are trained with early stopping: hold out a validation set, watch its loss each round, and stop when it turns. That is the tuning step a forest does not need.

Learning rate against tree count

Here is the evidence for the earlier claim that $\eta$ and $M$ are coupled but not interchangeable, at max_depth=2:

If the two were a single budget, every row would be reachable from every other by trading one against the other, and the best score would sit on a flat ridge. It does not. The best cell is lr=0.3 with 500 trees; lr=1.0 falls back from 0.977 to 0.971 as trees are added, because full-strength steps overfit however many you take, while lr=0.03 never catches up within this range. More rounds at a smaller rate is not the same journey taken slowly.

Depth

The measured version of the interaction-order argument, at 200 trees and lr=0.1:

The collapse from 0.982 to 0.918 between depth 3 and depth 5 is the cost of letting each correction express more than the data supports.

What XGBoost changes

Everything above is Friedman’s 1999–2001 gradient boosting machine. XGBoost is a 2016 reimplementation whose paper is titled “A Scalable Tree Boosting System”, and what it changes comes in three steps.

1. Classical boosting fits a tree to first-order gradients. That is the loop we built: compute $g_i = \partial L/\partial F$, fit a tree to it by variance reduction, then work out leaf values afterwards. Note the seam that leaves — the tree is grown by one criterion and filled by another, and the two have nothing to do with each other.

2. XGBoost derives both from one regularized second-order objective. It also computes the Hessian $h_i = \partial^2 L/\partial F^2$, and rather than fitting a tree to the gradient and patching up the leaves, it writes the new tree directly into the objective, Taylor-expands to second order, and minimizes. Two things fall out. The leaf value comes in closed form,

and the split criterion becomes the actual change in the objective, replacing variance reduction and closing the seam. The regularization lives inside the thing being minimized — $\lambda$ on leaf values, $\gamma$ as a fixed price per leaf — so pruning is arithmetic rather than a depth cap imposed from outside: a split that does not pay for itself is not taken.

You have already seen this idea at small scale. The $\sum p(1-p)$ in Friedman’s log-loss leaf value is a sum of Hessians, and $w_j^*$ is that same Newton step with a regularizer added to the denominator. First order says which way to go; second order says how far.

3. The systems work makes it scale. Approximate split finding from feature quantiles instead of scanning every distinct value; a quantile sketch weighted by the Hessians, so candidate thresholds cluster where the loss actually curves; pre-sorted compressed column blocks so the sort happens once rather than per split; cache prefetching and out-of-core sharding.

One of those deserves a specific callback. Sparsity-aware split finding gives every split a learned default direction: it enumerates only the present values, tries sending all the missing ones left, tries sending them all right, and keeps whichever scores better. The previous article ended on patient 267 — a real Cleveland patient with no recorded stress test, whose NA failed an equality test, slid down the False branch, and came back {'No': '100%'} while actually having the disease. Missing-ness stops being a preprocessing decision made in the dark and becomes a parameter learned per split:

# 21% of every feature blanked out, no imputation, no dropped rows
XGBClassifier(random_state=0, eval_metric="logloss").fit(X_train_with_nans, y_tr)
# test accuracy 0.942

The full derivation — the Taylor expansion, the structure score, and where $\gamma$ enters the split gain — is the subject of Inside XGBoost.

Costs, and when to reach for it

Three costs are real.

Interpretability. A single decision tree is a flowchart you can read and hand to a domain expert. Five hundred trees summed together is not, even though every individual piece is trivial. You get it back only through generic post-hoc tooling — SHAP values, partial dependence plots — which explain the model’s behavior rather than showing you its logic. If the model has to be defensible line by line to a regulator or a clinician, a shallow single tree may still win despite being worse.

Training is sequential. Tree $m$ cannot start until tree $m-1$ has updated the predictions. A random forest is embarrassingly parallel across trees; boosting is not, and its parallelism has to be found within each split search instead. This is a large part of what XGBoost’s systems work is for.

It must be tuned. Measured above: 0.918 down to 0.883 under label noise. There is no safe default number of rounds, and a validation set with early stopping is not optional.

Against those, the case for reaching for it: on tabular data — rows and columns where each column has a name and a meaning of its own — it remains extremely hard to beat, which is why it still runs fraud detection, delivery estimation, and search ranking in production in an era of large language models. It is not the tool when the features are pixels or tokens or a graph, because there the hard part is learning what the features are, and a sum of threshold questions has no mechanism for that. The previous article’s closing point stands: the right model is the one whose assumptions match your data.

Where this leaves us

Gradient boosting repeatedly asks a small tree to predict the direction in which the current model’s predictions should move. Squared error makes that direction look like an ordinary residual; any differentiable loss turns it into a general optimization method, run on functions instead of parameters; and XGBoost makes the step second-order, regularized, and scalable.

The code for everything here — the from-scratch regressor, the log-odds classifier, and the library comparisons — is in gbm.py, gbm_classifier.py and real_data.py alongside this article.