Greedy Function Approximation: A Gradient Boosting Machine, annotated

How to read this page

• Any dotted word explains itself when you hover it, tab to it, or tap it, and so does every symbol in every equation.
• The diagrams are live: hover or tap a part. The playground in §5 boosts small trees in your browser: switch the loss, drag the learning rate and the number of trees, and watch three bad labels pull the fit off course, or not.

Each idea climbs the same ladder: an everyday picture, a tiny example you can check by hand, a diagram, the math, and why it matters today. The running example is the one the trees and boosting lesson uses: four houses of size 1, 2, 3 and 4 that sold for 1, 2, 6 and 7. The lesson's GradientBoosting is this paper's Algorithm 1 in about thirty lines of NumPy.

Abstract · original

“Function estimation/approximation is viewed from the perspective of numerical optimization in function space, rather than parameter space.”Friedman (2001), Abstract

Everyday picture

A sculptor does not carve a statue in one blow. Each tap removes a little of what is still wrong, guided by where the stone most differs from the figure in mind. This paper says boosting is exactly that: every new small tree is a tap aimed at the direction in which the current model's mistakes shrink fastest. That direction is the gradient of the loss, and because every loss has one, the same recipe works for prices, for yes/no labels, for data full of outliers.

What the paper claims

• Adding one small model at a time, each chosen to reduce the loss, is gradient descent carried out on the predictions rather than on weights (§2 and §3).
• A single generic algorithm, Gradient_Boost, works for any differentiable loss: fit a tree by least squares to the negative gradient, then choose the step size for the actual loss (Algorithm 1).
• With trees as the small models, each leaf can get its own best step, giving the “TreeBoost” algorithms for squared error, absolute error, the Huber loss and logistic classification (§4).
• Shrinking every step by a learning rate ν makes the result markedly more accurate, at the price of more trees (§5).
• Boosted trees can be read through relative importance and partial dependence plots (§8), and inherit trees' robustness to messy data (§10).

Why it matters today

Gradient-boosted trees are the most common winning model on tabular data. XGBoost, LightGBM, CatBoost and scikit-learn's boosting all run this paper's loop, with refinements. See the XGBoost companion for the most influential of those refinements.

1 Function estimation · original

Everyday picture

Learning from examples means finding the rule F that turns inputs x (a house's size) into guesses for the output y (its price) so that, on average over all houses the world could produce, a chosen penalty for bad guesses is as small as possible. The penalty is the loss function L(y, F). Pick a different penalty and a different rule is best.

Tiny example

Suppose the rule may only be a single constant guess c for every house. With the squared penalty (y − c)², the best constant for prices 1, 2, 6 and 7 is their average, 4. With the absolute penalty |y − c|, any value between the two middle prices is best, and the median 4 is the usual pick. For these four prices the two agree; §4.2 shows a case where they do not.

In words: “the best rule is the one whose average penalty, over every input and output the world could produce, is smallest.”

With the numbers: with only constants allowed and the four prices as the whole world, the average squared penalty at c = 4 is (9 + 4 + 4 + 9)/4 = 6.5, lower than at c = 3 (7.5) or c = 5 (7.5). The average absolute penalty at c = 4 is (3 + 2 + 2 + 3)/4 = 2.5.

In Python:

y = [1, 2, 6, 7]
# E L(y, c): the average penalty of the constant guess c
def avg_squared(c):
    return sum((y_i - c) ** 2 for y_i in y) / len(y)
[avg_squared(c) for c in (3, 4, 5)]  # → [7.5, 6.5, 7.5]
# the absolute penalty at c = 4
sum(abs(y_i - 4) for y_i in y) / len(y)  # → 2.5

The paper then restricts F to an additive expansion: a sum of M simple pieces h(x; a_m), each with its own parameters a_m and weight β_m. Neural networks, splines and wavelets all have this shape; the paper's case of interest is where each piece is a small regression tree, whose parameters are its splitting variables, split points and leaf values.

In words: “the model is a weighted sum of M simple functions of the input, each shaped by its own parameters.”

With the numbers: one stump, h(x) = −2.5 if size ≤ 2.5 and +2.5 otherwise, with weight β = 1, added to a constant 4: F(4) = 4 + 1 × 2.5 = 6.5 and F(1) = 4 − 2.5 = 1.5.

In Python:

# h(x; a): a stump that splits size at 2.5
def h(x):
    return -2.5 if x <= 2.5 else 2.5
beta = 1.0
# F(x) = 4 + β h(x)
[4 + beta * h(x) for x in (1, 4)]  # → [1.5, 6.5]

1.1 and 1.2: optimizing parameters, one step at a time

With the shape fixed, learning becomes choosing parameters. Most methods do that by building the answer as a sum of steps: start from a guess, then repeatedly move the parameters a little in the direction that lowers the loss fastest, which is the negative gradient, by an amount found by a line search. That is steepest descent, the ancestor of every optimizer in the optimizers lesson.

Why it matters

Separating the loss L from the model F is the move that makes the rest of the paper general: every algorithm that follows is the same loop with a different L plugged in.

2 Numerical optimization in function space · original

Everyday picture

Instead of turning the knobs of a machine and watching all its outputs move at once, imagine reaching in and moving each output directly. If the machine predicts four house prices, treat those four predictions themselves as the four knobs. The slope of the loss with respect to each prediction says which way that one prediction should move.

Tiny example

Start from F = (4, 4, 4, 4) for the four houses. With the squared loss ½(y − F)², the slope for each house is F − y = (3, 2, −2, −3). Stepping the other way by a full unit moves the predictions to (1, 2, 6, 7): every price hit exactly. That is suspiciously perfect, and useless for a house of size 2.5, which is not one of the four knobs. §3 fixes this.

In words: “at every input, the next increment moves the prediction against the average slope of the loss there, by a step size ρ.”

With the numbers: for house 4, y = 7 and F₀ = 4, so the slope of ½(7 − F)² at F = 4 is 4 − 7 = −3, and with ρ = 1 the increment is +3: the prediction moves from 4 to 7.

In Python:

y, F, eps = 7.0, 4.0, 1e-6
# ∂L/∂F for L = ½(y − F)², measured by nudging F both ways
g = (0.5 * (y - (F + eps)) ** 2 - 0.5 * (y - (F - eps)) ** 2) / (2 * eps)
round(g, 6)  # → -3.0
rho = 1.0
# f_m = −ρ g_m
round(-rho * g, 6)  # → 3.0

Why it matters

This is the paper's central reframing. In ordinary training the gradient is taken with respect to weights; here it is taken with respect to the predictions. The same idea explains why, in the trees and boosting lesson, the residual y − F is called a negative gradient.

3 Finite data · original

“This permits the replacement of the difficult function minimization problem (9) by least-squares function minimization (11), followed by only a single parameter optimization based on the original criterion (12).”Friedman (2001), §3

Everyday picture

The gradient from §2 exists only at the houses in the training data: four arrows saying “this one up by 3, that one down by 2”. To predict for a house nobody has seen, the paper looks for the small tree whose predictions point most nearly in the same direction as those arrows. The tree is a smooth-enough stand-in for the arrows that also has an answer everywhere else.

Tiny example

The arrows (the negative gradient, called the pseudo-responses ỹ) for the four houses are −3, −2, +2, +3. Among all stumps, the one splitting at size 2.5 and predicting −2.5 on the left and +2.5 on the right fits them best by least squares (leftover squared error 1, against 14 for either other split). The line search then asks how far to move along that stump; for squared loss the answer is ρ = 1.

Reading it: each arrow starts at zero and shows how far, and which way, one house's prediction should move to reduce its squared error fastest: that is the negative gradient, −g. The two flat bars are the stump's answer: one value for every house left of the dashed split, one for every house right of it. The bars cannot match every arrow (they miss houses 1 and 4 by 0.5 each way), but they point the same way everywhere, and unlike the arrows they give an answer for any size at all. Fitting the bars to the arrows is equation (11); deciding how far to move along the bars is equation (12).

The math: Algorithm 1, line by line

Line 3 of Algorithm 1 computes the pseudo-responses:

In words: “for every training example i = 1 … N, work out how the loss would change if its prediction were nudged, and flip the sign: that is the target the next tree is trained on.”

With the numbers: squared loss ½(y − F)² at F₀ = 4: ỹ = y − 4 = (−3, −2, 2, 3).

In Python:

y = [1, 2, 6, 7]
F0 = sum(y) / len(y)
# ỹ_i = −∂L/∂F = y_i − F_0 for squared loss
y_tilde = [y_i - F0 for y_i in y]
y_tilde  # → [-3.0, -2.0, 2.0, 3.0]

Line 4 fits the base learner to them by least squares, whatever the loss:

In words: “choose the tree whose predictions, suitably scaled, come closest to the pseudo-responses in squared distance.”

With the numbers: the three possible stumps split after house 1, 2 or 3. Their leftover squared errors on (−3, −2, 2, 3) are 14, 1 and 14, so the split at 2.5 wins, with leaf values −2.5 and +2.5.

In Python:

y_tilde = [-3.0, -2.0, 2.0, 3.0]
def sse(v):
    mean = sum(v) / len(v)
    return sum((a - mean) ** 2 for a in v)
# leftover squared error after splitting after house t
[sse(y_tilde[:t]) + sse(y_tilde[t:]) for t in (1, 2, 3)]  # → [14.0, 1.0, 14.0]

Line 5 then goes back to the real loss to decide the step size:

In words: “slide along the new tree's direction and stop where the actual loss on the training data is lowest.”

With the numbers: for squared loss the best ρ is Σ ỹ_i h_i / Σ h_i² = (7.5 + 5 + 5 + 7.5) / 25 = 1.0.

In Python:

y_tilde = [-3.0, -2.0, 2.0, 3.0]
h = [-2.5, -2.5, 2.5, 2.5]
# for squared loss the line search has a closed form: Σ ỹ h / Σ h²
sum(a * b for a, b in zip(y_tilde, h)) / sum(b * b for b in h)  # → 1.0

Reading it: start at the top: the model begins as the single best constant. Everything inside the dashed frame repeats once per tree. Only the first and third steps know which loss you chose: the pseudo-responses are that loss's negative gradient, and the step size is chosen on that loss. The middle step is always the same least-squares tree fit, which is why one tree-growing routine serves every loss. The arrow on the left carries the updated model back to the top for the next round; after M rounds the model is the constant plus M shrunken trees.

Why it matters

“Least squares to find the direction, the real loss to find the step” is the paper's engineering insight. It means one fast tree-growing routine can boost any differentiable loss. The lesson's GradientBoosting.negative_gradient is line 3; its call to DecisionTree with the squared criterion is line 4.

4 Applications: additive modeling · original

Everyday picture

Algorithm 1 is a machine with a slot for the loss. This section feeds it four different losses and writes out what comes out: least squares (a sanity check), least absolute deviation (robust), Huber (a compromise) and the logistic likelihood (classification).

import statistics
y = [1, 2, 6, 70]
# least squares starts at the mean, and the typo dominates the residuals
F0_ls = sum(y) / len(y)
[round(y_i - F0_ls, 2) for y_i in y]  # → [-18.75, -17.75, -13.75, 50.25]
# LAD starts at the median; its pseudo-responses are signs
F0 = statistics.median(y)
[(y_i > F0) - (y_i < F0) for y_i in y]  # → [-1, -1, 1, 1]

# the left leaf: houses 1 and 2, current prediction 4
y_leaf, F_prev = [1, 2], 4.0
def leaf_loss(gamma):
    return sum(0.5 * (y_i - (F_prev + gamma)) ** 2 for y_i in y_leaf)
[leaf_loss(g) for g in (-3.0, -2.5, -2.0)]  # → [0.5, 0.25, 0.5]

residuals = [-3.0, -2.0, 2.0, 66.0]
delta = 5.0
# −∂L/∂F: the residual, clipped to ±δ
[max(-delta, min(delta, r)) for r in residuals]  # → [-3.0, -2.0, 2.0, 5.0]

import math
y = [1, 1, 1, -1]
y_bar = sum(y) / len(y)
# F_0 = ½ log((1 + ȳ) / (1 − ȳ)): half the log-odds of the base rate
F0 = 0.5 * math.log((1 + y_bar) / (1 - y_bar))
round(F0, 3)  # → 0.549
# ỹ_i = 2 y_i / (1 + exp(2 y_i F_0))
[round(2 * y_i / (1 + math.exp(2 * y_i * F0)), 3) for y_i in y]  # → [0.5, 0.5, 0.5, -1.5]

import math
F0 = 0.5 * math.log(3)
def newton_step(leaf):
    # Σ ỹ / Σ |ỹ|(2 − |ỹ|)
    return sum(leaf) / sum(abs(t) * (2 - abs(t)) for t in leaf)
round(newton_step([0.5, 0.5]), 3)  # → 0.667
round(newton_step([-1.5]), 3)  # → -2.0
# probabilities from F: p = 1 / (1 + e^(−2F))
[round(1 / (1 + math.exp(-2 * (F0 + g))), 3) for g in (newton_step([0.5, 0.5]), newton_step([-1.5]))]  # → [0.919, 0.052]

# w_i = |ỹ_i| (2 − |ỹ_i|)
[round(abs(t) * (2 - abs(t)), 4) for t in (0.5, -1.5, 0.02)]  # → [0.75, 0.75, 0.0396]

import math
F = [0.0, 0.0, 0.0]
# p_k = exp(F_k) / Σ_l exp(F_l)
p = [math.exp(F_k) / sum(math.exp(F_l) for F_l in F) for F_k in F]
y = [1, 0, 0]
# ỹ_k = y_k − p_k
[round(y_k - p_k, 3) for y_k, p_k in zip(y, p)]  # → [0.667, -0.333, -0.333]

5 Regularization · original

“As illustrated here decreasing the learning rate clearly improves performance, usually dramatically. The reason for this is less clear.”Friedman (2001), §5

Everyday picture

Every tree is fitted to the training data, noise included. Taking each tree's full correction lets the first few trees commit hard to whatever they saw. Taking only a fraction of each correction means many trees must agree before the model moves far, which averages some of the noise away. The price is more trees.

Tiny example

With a learning rate ν = 0.1, house 4 moves from 4 to 4 + 0.1 × 2.5 = 4.25 after the first tree, instead of all the way to 6.5. The remaining gap is left for later trees, which will see a slightly different picture.

In words: “each round adds only a fraction ν of the best step the new tree offers.”

With the numbers: house 4: 4 + 0.1 × 1.0 × 2.5 = 4.25.

In Python:

F_prev, nu, rho, h = 4.0, 0.1, 1.0, 2.5
# F_m = F_{m−1} + ν ρ_m h
round(F_prev + nu * rho * h, 2)  # → 4.25

Try it: boost in your browser

Forty noisy points around a sine curve, three of them pushed up by a labelling mistake (hollow circles). Each round adds one small tree with two levels of questions (up to four leaves). Choose the loss, drag ν and the number of trees M, and watch the fit (solid line) against the true curve (dashed). The chart below the picture tracks how far the fit is from the true curve after every round.

Reading it: in the picture, dots are training points and the solid line is the model after M trees; it is always a staircase, because every tree adds flat steps. The chart below plots the fit's distance from the true curve against the number of trees, so its lowest point is the best M for your settings. With least squares, the three hollow outliers pull a bump into the fit near them; switch to least absolute deviation or Huber and the bump mostly disappears, because those losses cap how loudly one point can pull. Now drag ν: at 1.0 the curve bottoms out within a few trees and then worsens as trees chase noise; at 0.1 it falls slowly, reaches a lower minimum, and stays near it much longer. That is the paper's ν-M trade-off in miniature.

What the paper measured

The paper boosts 11-leaf trees on 5,000 simulated examples and scores each model with an error measure relative to the best constant, its equation (37):

In words: “the model's average distance from the true function, as a fraction of the distance a constant guess (the median) would have.” 0 is perfect; 1 is no better than a constant.

With the numbers: treat the four prices as the true function and the one-round fit (2.75, 2.75, 5.25, 5.25) as F̂: the top is (1.75 + 0.75 + 0.75 + 1.75)/4 = 1.25, the bottom (3 + 2 + 2 + 3)/4 = 2.5, so A = 0.5: the model has closed half the gap a constant leaves.

In Python:

import statistics
F_star = [1, 2, 6, 7]
F_hat = [2.75, 2.75, 5.25, 5.25]
top = sum(abs(a - b) for a, b in zip(F_star, F_hat)) / 4
bottom = sum(abs(a - statistics.median(F_star)) for a in F_star) / 4
top, bottom, top / bottom  # → (1.25, 2.5, 0.5)

Reading it: each line is one of the paper's methods from its Table 1: least squares and LAD scored by A, and the logistic model scored by its deviance (−2 log-likelihood). The horizontal axis is the learning rate ν on a log scale, and the vertical axis is the number of trees at which the error was lowest. Every line climbs steeply to the left: for least squares, halving ν roughly doubles the trees needed (15, 43, 77, 146, 326, 855). The table below shows what that buys; the logistic model's deviance falls in the same way, from 0.60 at ν = 1 to 0.45 at ν = 0.125.

Going from ν = 1 to ν = 0.125 cuts the least-squares error by a third (0.48 to 0.32); going further buys nothing but trees. The paper also notes that the misclassification rate kept improving after the likelihood had started to overfit, since a class decision depends only on the sign of F, while the likelihood also cares about its size. And it records a mystery it could not explain: shrinking each step during training (36) works far better than shrinking the finished model towards the mean (38).

Why it matters

The paper's advice survives unchanged: set M as large as you can afford, then pick ν so the validation error reaches its minimum near the end, and stop there (early stopping). The lesson measures the same trade-off with boosting_curves.

6 Simulation studies · original

Everyday picture

One benchmark dataset flatters whichever method happens to suit it. The paper instead generates 100 random target functions of 10 inputs (§6.1: sums of 20 random bumps, each depending on a few inputs), and asks how each method does across all of them.

What it found

• Loss versus noise (§6.2). With normal noise, least squares was best on 73 of the 100 targets and Huber (M_TreeBoost) on the other 27; on average least squares was 0.2% worse than the best method, Huber 0.9% and LAD 7.4%. With the heavy-tailed “slash” noise the picture flips: least squares was 364.6% worse than the best on average, LAD 4.1% and Huber 1.0%. The paper's verdict: Huber is the method of choice, never far from the best.
• Steps versus smooth curves (§6.3). Boosted trees are staircases; MARS fits smooth curves. On the 100 smooth targets the two were comparable by average absolute error. But for MARS the root-mean-squared error was typically about 30% higher than its average absolute error, while for the boosted trees the two were close: MARS tends to be either very close or far off, and the boosted trees' errors are more even.
• Classification (§6.4). On five-class problems, L_K_TreeBoost was best on 78 of 100 targets, 0.6% above the best on average, against 3.5% for LogitBoost and 15% for AdaBoost.MH. Once LogitBoost and AdaBoost.MH were also given shrinkage (ν = 0.1), all three came within a few percent of each other, which suggests the learning rate matters more than the structural differences.

Why it matters

Two lessons that still hold: the loss should match the noise, and shrinkage helps every boosting method, not just this one.

7 Tree boosting · original

Everyday picture

A house's price may depend on its size and its location separately (each adds its own amount), or on how they combine (a big house matters more downtown). A tree with J leaves can express combinations of at most J − 1 inputs. Stumps (J = 2) can only add up separate effects; bigger trees can capture combinations.

Tiny example

Write any function as a sum of parts that each depend on one input, parts that depend on pairs, parts on triples, and so on:

In words: “a function is its one-input effects, plus its two-input interactions, plus higher ones”; a tree with J leaves reaches at most order min(J − 1, n).

With the numbers: the paper's simulations have n = 10 inputs. Trees with J = 2, 3, 6, 11 and 21 leaves reach interactions of order 1, 2, 5, 10 and 10.

In Python:

n = 10
# highest interaction order a J-leaf tree can express: min(J − 1, n)
[min(J - 1, n) for J in (2, 3, 6, 11, 21)]  # → [1, 2, 5, 10, 10]

Reading it: each bar is one tree size J from the paper's §7 (the numbers behind its Figure 6), over the 100 random targets: its length is how much larger, on average, that tree size's error was than the best tree size's on the same target. The number on the right is how many of the 100 targets that tree size won. Stumps (J = 2) won 8 times but were 23.2% worse on average: superb when the target happens to be additive, poor otherwise. From J = 6 up, the sizes are nearly tied and each wins about 30 targets.

Why it matters

This is why boosting uses small trees: a handful of leaves captures the interactions most real targets have, and the paper finds it unlikely that large trees would ever be necessary or desirable. The lesson boosts depth-2 trees (four leaves) in boosting_curves.

8 Interpretation · original

Everyday picture

A model made of hundreds of trees cannot be read tree by tree. The paper offers two summaries: which inputs matter most (relative importance), and how the prediction moves as one or two inputs change with everything else averaged out (partial dependence).

w_l, w_r, y_l, y_r = 2, 2, -2.5, 2.5
# i² = w_l w_r / (w_l + w_r) · (ȳ_l − ȳ_r)²
w_l * w_r / (w_l + w_r) * (y_l - y_r) ** 2  # → 25.0

import math
# (splitting variable v_t, improvement î_t²) for the three splits
splits = [(1, 4.0), (2, 1.0), (1, 0.5)]
I = {j: math.sqrt(sum(i2 for v_t, i2 in splits if v_t == j)) for j in (1, 2)}
round(I[1], 2)  # → 2.12
# relative importance, top input = 100
round(100 * I[2] / I[1], 1)  # → 47.1

def F_hat(x1, x2):
    return 2 * x1 + x1 * x2
# z_{i,\l}: the other input's value in each data row
x2_rows = [0, 1, 2]
# F̄(z) = (1/N) Σ F̂(z, x2_i)
[sum(F_hat(z, x2) for x2 in x2_rows) / len(x2_rows) for z in (1, 2)]  # → [3.0, 6.0]

8.3 A randomly generated function

Applied to the first of the §6 random targets, the importances show no small dominant subset of inputs, the one-input plots show smooth trends built from fine steps, and the two-input plots show interactions of varying strength, all consistent with how the target was generated.

9 Real data · original

Garnet data (§9.1)

13,317 garnets from rocks around the world; the task is to predict titanium dioxide concentration from 10 other chemical concentrations and a three-valued tectonic setting. The data are skewed with many outliers.

Stumps are clearly worse, so the inputs interact; six leaves are enough; and the robust losses beat least squares on this outlier-heavy data. Two inputs, gallium and zirconium, dominate the importances, and the two-input partial dependence shows a strong interaction between them: the prediction barely depends on either while the other is at its lowest values.

Demographic data (§9.2)

9,409 shopping-mall questionnaires; the task is to predict income from 13 other answers, most categorical with a few values, many missing. Here tree size hardly matters (errors 0.57 to 0.63 across all sizes and losses), so stumps suffice: the relationship is additive, and the partial dependence plots then describe the model completely. Income rises with education and peaks at age 45 to 54.

Why it matters

The same diagnostic still works: compare stumps with larger trees on validation data. If stumps are as good, the problem is additive and every effect can be plotted one input at a time.

10 Data mining · original

“All TreeBoost procedures are invariant under all (strictly) monotone transformations of the individual input variables.”Friedman (2001), §10

Everyday picture

Real data arrive messy: skewed incomes, outliers, irrelevant columns, missing answers. Boosted trees shrug off most of it, because a tree only ever asks “is this input above that threshold?”.

What the paper lists

• No input transformations needed. Using x, log x, e^x or x^a gives the same result, since the order of values is all a split sees; long-tailed inputs and input outliers stop mattering. LAD_TreeBoost is also fully robust to outliers in the output, and M_TreeBoost largely so.
• Built-in feature selection and robustness to irrelevant inputs, and trees' handling of missing values.
• The weaknesses of single trees are softened. Coarse staircases become fine ones, instability is averaged away, and the interaction order is controlled by the tree size.
• Speed. After sorting the inputs, the cost grows linearly with the number of examples, inputs and trees. M_TreeBoost on the garnet data (500 trees of 6 leaves) took 20 seconds on a 933 MHz Pentium III.
• Incremental use. Good-enough models often appear by about 100 trees, and boosting can be resumed later, or continued on new data, from where it stopped.

Why it matters

This list is, almost word for word, why boosted trees became the default on tables: the same reasons the trees and boosting lesson gives for trees beating neural networks there (units don't matter, thresholds are native, irrelevant columns are ignored).

What changed since 2001

The loop is unchanged: pseudo-responses, a least-squares tree, per-leaf steps, shrinkage. The refinements around it:

Glossary

Every term with hover guidance on this page, in one place.