Gradient boosting

What it is and when you reach for it

Gradient boosting builds a strong predictor as a sum of many weak ones, F(x) = F₀ + ν Σ_m h_m(x), where each h_m is a shallow decision tree trained to correct the errors of everything before it. It sits opposite random forests on the ensemble map: forests reduce variance by averaging independent deep trees trained in parallel, while boosting reduces bias by adding dependent shallow trees trained sequentially, each one aimed precisely at what the current ensemble still gets wrong. Against neural networks the division of labor is by data type: on images, audio, and text, representation learning wins, but on medium-sized heterogeneous tabular data, the mix of counts, amounts, categories, and ratios that fills most business databases, gradient-boosted trees have been the strongest general-purpose model for two decades and still routinely beat deep learning baselines. You reach for it whenever the data is a table, the features mean different things, and you want near-best accuracy with minutes of training and little preprocessing: no scaling, no imputation ceremony, no architecture search. Its weak spots are the same as any tree model's, no extrapolation and no transfer of learned representations, plus one of its own: the sequential construction means mistakes compound if you overfit, which is why the learning rate and early stopping matter so much here.

The math

Boosting as gradient descent in function space

Start from the objective: minimize L(F) = Σ_i ℓ(y_i, F(x_i)) over functions F, not parameter vectors. Since the training loss only touches F through its values on the n training points, we can treat F as a vector F = (F(x₁), ..., F(x_n)) and ask what plain gradient descent would do: step in the direction of the negative gradient, whose i-th component is g_i = ∂ℓ(y_i, F(x_i)) / ∂F(x_i). The catch is that this update is defined only at the training points, and a model must predict everywhere. Friedman's resolution, in the 1999 paper that named the algorithm, is to approximate the gradient step with a function from a restricted class: fit a regression tree h_m to the targets −g_i by least squares, so h_m is the projection of the negative gradient onto the space of small trees, then update F_m = F_m−1 + ν h_m. Each boosting round is therefore literally one step of gradient descent in function space, with the tree acting as a smoother that generalizes the per-sample gradient to all of feature space, and ν playing exactly the role of a learning rate. Everything else in the algorithm is a choice of loss or a refinement of the step.

Squared loss: boosting is residual fitting

For ℓ = ½(y − F)² the negative gradient is −g_i = y_i − F(x_i), the plain residual. So least-squares boosting has a completely elementary description: fit a small tree to the residuals, add it, recompute residuals, repeat. A worked round makes it concrete. Take four points with x = (1, 2, 3, 4) and y = (1.0, 1.5, 3.0, 6.0). Initialize with the best constant, F₀ = mean(y) = 2.875, giving residuals r = (−1.875, −1.375, 0.125, 3.125) and squared-error 15.19. Fit a stump: of the three candidate splits, x ≤ 1 leaves error 10.50, x ≤ 2 leaves 4.63, and x ≤ 3 leaves 2.17, so the stump splits at x ≤ 3, predicting the left mean −1.042 and the right mean 3.125. With shrinkage ν = 0.1 the update is F₁ = F₀ + 0.1 · h₁, i.e. 2.771 on the left and 3.188 on the right, and the residual error drops from 15.19 to 12.71. Run two hundred more rounds of exactly this and the additive staircase hugs the data:

F0: constant        F0 + v*h1: one step        after many rounds
      *                    *                          *
  ----------          ______*----                ____*
      *                *   /                       _*
   *                *_____/                     __*
 *__________       *_________          *______*

Log loss: gradients for classification

For binary classification keep the model on the logit scale, p = σ(F), and use the log loss, which for y ∈ {0, 1} simplifies to ℓ(y, F) = log(1 + e^F) − yF. Differentiating with respect to F gives g = σ(F) − y = p − y, so the pseudo-residual is y − p: the tree at each round is fit to the gap between the label and the current predicted probability, exactly parallel to the regression case but on the logit scale. One subtlety earns a flag because it is the classic implementation bug: the leaf value for log-loss boosting is not the mean of the pseudo-residuals in the leaf. The tree's structure is found by least-squares fitting of the pseudo-residuals, but the value assigned to each leaf should be the step that actually minimizes the log loss, which one Newton step gives as Σ(y − p) / Σ p(1 − p) over the leaf's samples. The second derivative h = p(1 − p) appears here for the first time and takes center stage two sections down. Initialization follows the same logic as regression: F₀ is the log-odds of the base rate, log(p̄/(1 − p̄)).

Shrinkage and subsampling

The learning rate ν scales every tree's contribution, and the empirical law, established already in Friedman's papers, is that smaller ν with proportionally more trees almost always generalizes better: ν in 0.03 to 0.1 with hundreds to thousands of trees is the standard operating point, versus ν = 1 with few trees, which overfits the early trees' mistakes. The intuition is regularization by small steps. Each tree overfits its pseudo-residuals a little; shrinking its contribution means later trees, fit on data-recomputed gradients, get to vote on and partially cancel that noise before it hardens into the model. Stochastic gradient boosting (Friedman, 2002) adds the second knob: fit each tree on a random subsample of rows (without replacement, typically 50 to 80 percent). This decorrelates consecutive trees, adds the same kind of variance reduction bagging provides, and speeds up each round; column subsampling per tree or per split, borrowed from random forests, does the same across features and is standard in all the modern libraries.

Second-order boosting: the XGBoost step

First-order boosting uses only the gradient and relies on the least-squares fit and a fixed ν to size the step. The refinement that XGBoost systematized is to expand the loss to second order around the current prediction: writing f for the new tree's step, ℓ(y, F + f) ≈ ℓ(y, F) + g·f + ½ h·f², with per-sample gradient g_i and hessian h_i, and to choose the tree to minimize this quadratic plus explicit regularization λ‖w‖² on the leaf values. The quadratic solves in closed form per leaf: with G = Σ g_i and H = Σ h_i over a leaf's samples, the optimal leaf value is w* = −G / (H + λ), and plugging it back gives the leaf's contribution to the objective, −½ G² / (H + λ). That yields a split criterion that replaces impurity entirely: the gain of a split is ½ [ G_L²/(H_L + λ) + G_R²/(H_R + λ) − G²/(H + λ) ] − γ, where γ prices each additional leaf. The conceptual upgrade is that the tree grower now optimizes the actual boosting objective, curvature included, rather than a least-squares proxy, and every leaf takes a Newton step rather than a gradient step scaled by a global constant. For squared loss h ≡ 1 and Newton boosting collapses back to residual fitting, which is a good sanity check on the whole derivation. The implementations below use this (g, h) formulation directly, because it handles regression and classification with the same twenty lines.

Implementation, twice

The same honesty note as on the decision trees page applies: trees are not what autodiff frameworks are for, and these implementations exist for pedagogy, with the recursion for tree growth in plain Python and the split search vectorized as a histogram over quantile thresholds, now accumulating gradient and hessian sums instead of class counts. But boosting adds one place where the framework genuinely earns its keep: the per-sample gradients and hessians themselves. The JAX version gets g and h by differentiating the loss with jax.grad composed twice and vectorized with vmap, so switching the entire algorithm to a new loss means writing one scalar function; the PyTorch version does the same for g via autograd and uses the closed-form hessian. That is the function-space-gradient-descent story made executable: autodiff supplies the direction, and the tree is the step.

import torch

class Node:
    __slots__ = ("feature", "threshold", "left", "right", "value")
    def __init__(self, value):
        self.feature = None
        self.threshold = None
        self.left = None
        self.right = None
        self.value = value

class GradientBoosting:
    """Newton gradient boosting with short trees (binary clf or regression).

    Split search is a histogram of (grad, hess) sums over per-feature
    quantile thresholds, scored with the XGBoost gain; leaf values are
    Newton steps -G / (H + lam). Autograd supplies g; h is closed-form.
    Pedagogical: the tree recursion is Python, not a fused C++ grower.
    """

    def __init__(self, n_trees=300, lr=0.1, max_depth=3, subsample=0.8,
                 lam=1.0, min_leaf=5, n_bins=32, loss="logloss"):
        self.n_trees, self.lr, self.max_depth = n_trees, lr, max_depth
        self.subsample, self.lam, self.min_leaf = subsample, lam, min_leaf
        self.n_bins, self.loss = n_bins, loss

    def _loss(self, F, y):
        if self.loss == "mse":
            return 0.5 * (y - F) ** 2
        return torch.nn.functional.softplus(F) - y * F  # log-loss on logits

    def _grads(self, F, y):
        F = F.detach().requires_grad_(True)
        g = torch.autograd.grad(self._loss(F, y).sum(), F)[0]
        if self.loss == "mse":
            h = torch.ones_like(g)
        else:
            p = torch.sigmoid(F.detach())
            h = p * (1 - p)
        return g.detach(), h

    def _best_split(self, X, g, h):
        left = (X.unsqueeze(-1) <= self.thresholds.unsqueeze(0)).float()
        GL = torch.einsum("ndt,n->dt", left, g)      # (d, T) grad sums
        HL = torch.einsum("ndt,n->dt", left, h)
        nL = left.sum(0)
        G, H = g.sum(), h.sum()
        gain = (GL ** 2 / (HL + self.lam)
                + (G - GL) ** 2 / (H - HL + self.lam)
                - G ** 2 / (H + self.lam))
        valid = (nL >= self.min_leaf) & (len(g) - nL >= self.min_leaf)
        gain = torch.where(valid, gain, torch.tensor(float("-inf")))
        j, t = divmod(int(gain.argmax()), gain.shape[1])
        if gain[j, t] <= 0:
            return None
        return j, float(self.thresholds[j, t])

    def _grow(self, X, g, h, depth):
        node = Node(value=float(-g.sum() / (h.sum() + self.lam)))  # Newton step
        if depth >= self.max_depth or len(g) < 2 * self.min_leaf:
            return node
        split = self._best_split(X, g, h)
        if split is None:
            return node
        node.feature, node.threshold = split
        m = X[:, node.feature] <= node.threshold
        node.left = self._grow(X[m], g[m], h[m], depth + 1)
        node.right = self._grow(X[~m], g[~m], h[~m], depth + 1)
        return node

    def _tree_predict(self, node, X):
        out = torch.empty(len(X))
        stack = [(node, torch.arange(len(X)))]
        while stack:
            nd, idx = stack.pop()
            if nd.feature is None:
                out[idx] = nd.value
                continue
            m = X[idx, nd.feature] <= nd.threshold
            stack.append((nd.left, idx[m]))
            stack.append((nd.right, idx[~m]))
        return out

    def fit(self, X, y):
        n = len(y)
        y = y.float()
        qs = torch.linspace(0, 1, self.n_bins + 1)[1:-1]
        self.thresholds = torch.quantile(X, qs, dim=0).T   # (d, T), computed once
        if self.loss == "mse":
            self.f0 = float(y.mean())
        else:
            p = y.mean().clamp(1e-6, 1 - 1e-6)
            self.f0 = float(torch.log(p / (1 - p)))        # log-odds base rate
        F = torch.full((n,), self.f0)
        self.trees = []
        for _ in range(self.n_trees):
            g, h = self._grads(F, y)
            idx = torch.randperm(n)[: int(self.subsample * n)]
            tree = self._grow(X[idx], g[idx], h[idx], depth=0)
            F = F + self.lr * self._tree_predict(tree, X)  # update ALL rows
            self.trees.append(tree)
        return self

    def decision_function(self, X):
        F = torch.full((len(X),), self.f0)
        for tree in self.trees:
            F = F + self.lr * self._tree_predict(tree, X)
        return F

    def predict(self, X):
        F = self.decision_function(X)
        return (F > 0).long() if self.loss == "logloss" else F

import jax
import jax.numpy as jnp
from functools import partial

def logloss(f, y):
    """Per-sample log loss on the logit scale: softplus(f) - y*f."""
    return jnp.logaddexp(0.0, f) - y * f

def mse(f, y):
    return 0.5 * (y - f) ** 2

# The one place autodiff genuinely earns its keep in a tree ensemble:
# per-sample g and h come from differentiating the scalar loss twice
# and vmapping. Swap the loss function and the whole algorithm follows.
def make_grads(loss):
    g_fn = jax.vmap(jax.grad(loss))
    h_fn = jax.vmap(jax.grad(jax.grad(loss)))
    return lambda F, y: (g_fn(F, y), h_fn(F, y))

@partial(jax.jit, static_argnames="min_leaf")
def best_split(X, g, h, thresholds, lam=1.0, min_leaf=5):
    """XGBoost gain over every (feature, threshold) pair in one shot."""
    left = (X[:, :, None] <= thresholds[None]).astype(g.dtype)  # (n, d, T)
    GL = jnp.einsum("ndt,n->dt", left, g)
    HL = jnp.einsum("ndt,n->dt", left, h)
    nL = left.sum(0)
    G, H = g.sum(), h.sum()
    gain = GL**2 / (HL + lam) + (G - GL)**2 / (H - HL + lam) - G**2 / (H + lam)
    valid = (nL >= min_leaf) & (X.shape[0] - nL >= min_leaf)
    gain = jnp.where(valid, gain, -jnp.inf)
    j, t = jnp.unravel_index(jnp.argmax(gain), gain.shape)
    return j, thresholds[j, t], gain[j, t]

def grow(X, g, h, thresholds, depth, max_depth=3, lam=1.0, min_leaf=5):
    value = -g.sum() / (h.sum() + lam)          # Newton leaf value
    if depth >= max_depth or X.shape[0] < 2 * min_leaf:
        return {"value": float(value)}
    j, thr, gain = best_split(X, g, h, thresholds, lam=lam, min_leaf=min_leaf)
    if not bool(gain > 0):
        return {"value": float(value)}
    m = X[:, j] <= thr
    kw = dict(max_depth=max_depth, lam=lam, min_leaf=min_leaf)
    return {"feature": int(j), "threshold": float(thr),
            "left":  grow(X[m], g[m], h[m], thresholds, depth + 1, **kw),
            "right": grow(X[~m], g[~m], h[~m], thresholds, depth + 1, **kw)}

def tree_predict(tree, X):
    # Branch-free: evaluate every leaf for every row, select with where.
    if "feature" not in tree:
        return jnp.full(X.shape[0], tree["value"])
    m = X[:, tree["feature"]] <= tree["threshold"]
    return jnp.where(m, tree_predict(tree["left"], X),
                        tree_predict(tree["right"], X))

def fit(X, y, key, loss=logloss, n_trees=300, lr=0.1, max_depth=3,
        subsample=0.8, lam=1.0, min_leaf=5, n_bins=32):
    n = X.shape[0]
    y = y.astype(jnp.float32)
    qs = jnp.linspace(0, 1, n_bins + 1)[1:-1]
    thresholds = jnp.quantile(X, qs, axis=0).T
    if loss is mse:
        f0 = float(y.mean())
    else:
        p = jnp.clip(y.mean(), 1e-6, 1 - 1e-6)
        f0 = float(jnp.log(p / (1 - p)))
    grads = make_grads(loss)
    F = jnp.full(n, f0)
    trees = []
    for _ in range(n_trees):
        key, sub = jax.random.split(key)
        g, h = grads(F, y)
        idx = jax.random.choice(sub, n, (int(subsample * n),), replace=False)
        tree = grow(X[idx], g[idx], h[idx], thresholds, depth=0,
                    max_depth=max_depth, lam=lam, min_leaf=min_leaf)
        F = F + lr * tree_predict(tree, X)     # update all rows, not the subsample
        trees.append(tree)
    return f0, trees

def decision_function(f0, trees, X, lr=0.1):
    F = jnp.full(X.shape[0], f0)
    for t in trees:
        F = F + lr * tree_predict(t, X)
    return F

Notice what changed from the classification tree on the decision trees page and what did not: the histogram machinery is identical, but class counts became (G, H) sums, Gini became the Newton gain, and majority-vote leaves became −G/(H + λ) steps. Boosting is a different objective flowing through the same tree grower.

Using it on a real shape of problem

A realistic shape: 6,000 rows, 25 features, binary label with a nonlinear decision rule and noise, roughly a small risk-model dataset.

import torch

torch.manual_seed(0)
n, d = 6000, 25
X = torch.randn(n, d)
logit = X[:, 0] * X[:, 1] + X[:, 2] ** 2 - 1 + 0.5 * torch.randn(n)
y = (logit > 0).long()                      # interaction + curvature: tree food
X_tr, y_tr, X_te, y_te = X[:4500], y[:4500], X[4500:], y[4500:]

model = GradientBoosting(n_trees=300, lr=0.1, max_depth=3).fit(X_tr, y_tr)
acc = (model.predict(X_te) == y_te).float().mean()
print(f"test accuracy {acc:.3f}")

# staged test loss: the curve to actually look at
F = torch.full((len(y_te),), model.f0)
for i, tree in enumerate(model.trees):
    F = F + model.lr * model._tree_predict(tree, X_te)
    if (i + 1) % 50 == 0:
        loss = (torch.nn.functional.softplus(F) - y_te * F).mean()
        print(f"trees {i+1:3d}  test logloss {loss:.4f}")

import jax
import jax.numpy as jnp

key = jax.random.PRNGKey(0)
n, d = 6000, 25
kx, kn, kf = jax.random.split(key, 3)
X = jax.random.normal(kx, (n, d))
logit = X[:, 0] * X[:, 1] + X[:, 2] ** 2 - 1 + 0.5 * jax.random.normal(kn, (n,))
y = (logit > 0).astype(jnp.int32)           # interaction + curvature: tree food
X_tr, y_tr, X_te, y_te = X[:4500], y[:4500], X[4500:], y[4500:]

f0, trees = fit(X_tr, y_tr, kf, n_trees=300, lr=0.1, max_depth=3)
F = decision_function(f0, trees, X_te)
acc = ((F > 0).astype(jnp.int32) == y_te).mean()
print(f"test accuracy {acc:.3f}")

# staged test loss: the curve to actually look at
F = jnp.full(len(y_te), f0)
for i, tree in enumerate(trees):
    F = F + 0.1 * tree_predict(tree, X_te)
    if (i + 1) % 50 == 0:
        loss = jnp.mean(logloss(F, y_te.astype(jnp.float32)))
        print(f"trees {i+1:3d}  test logloss {loss:.4f}")

Expect, with machine- and seed-dependent numbers: test log loss falling steeply over the first 50 to 100 trees, then flattening, and eventually creeping upward if you keep adding trees, which is the overfitting signature that early stopping exists to catch. Test accuracy on data like this should land in the high 0.80s to low 0.90s, clearly above a depth-limited single tree and far above logistic regression, which cannot represent the x₀x₁ interaction at all. Halve the learning rate and the curve shifts right, needing roughly twice the trees to reach a similar or slightly better floor: shrinkage and ensemble size trade off almost exactly.

Applications

Gradient boosting's home turf is tabular machine learning, and the dominance is not folklore: for a decade the winning solutions of most tabular Kaggle competitions have been boosted-tree ensembles, and systematic comparisons (Grinsztajn, Oyallon, and Varoquaux's 2022 NeurIPS benchmark study is the standard reference) keep finding that tuned tree ensembles beat tuned deep models on typical medium-sized tables. Ranking is the second pillar: LambdaMART, the method behind production web search rankers for years (it came out of Microsoft research and powered Bing's ranker, and its descendants live in LightGBM's ranking objectives), is gradient boosting with a twist worth knowing at concept level. Ranking metrics like NDCG are non-differentiable, so LambdaMART defines the gradients directly: for each pair of documents ordered incorrectly, a lambda proportional to the NDCG change from swapping them, and boosting proceeds on those synthetic gradients, a vivid demonstration that the algorithm only ever needs a gradient signal, not an explicit loss. The third pillar is risk: fraud detection at payment processors, credit default scoring, insurance pricing, and churn models are overwhelmingly boosted trees, because the data is tabular, the classes are imbalanced (handled cleanly through the loss), and monotonic constraints let modelers encode business rules like "risk must not decrease as delinquencies increase". My credit risk model project is exactly this shape of problem end to end. Beyond those, boosted trees show up as the tabular head on top of learned embeddings, as rerankers in recommendation pipelines, and in particle physics classifiers, where BDTs were standard at the LHC before deep learning arrived.

Against the real libraries

Four libraries cover essentially all production use. XGBoost is the system that made second-order boosting standard: the regularized Newton objective derived above, sparsity-aware split finding with a learned default direction for missing values, approximate and histogram tree methods, GPU training, and monotonic and interaction constraints. LightGBM attacked speed and won: histogram-based splits (the trick this page's implementations borrow), leaf-wise rather than level-wise growth, which spends the leaf budget where the gain is instead of filling depth uniformly, GOSS (gradient-based one-side sampling, which keeps all large-gradient rows and subsamples the small-gradient ones, since well-fit rows carry little signal about where to split next), and EFB, which bundles mutually exclusive sparse features into single columns. CatBoost contributed two ideas about leakage: ordered target statistics for categorical features, encoding each row's categories using only rows that come before it in a random permutation so a category's encoding never contains its own label, and ordered boosting, which fights the subtler self-leakage where residuals are computed by a model that has already seen the very labels the residuals are fit to. Finally, scikit-learn's HistGradientBoostingClassifier and its regressor twin, explicitly inspired by LightGBM, bring histogram boosting with native missing-value and categorical support into the standard sklearn API, and are the right default when you are already in that ecosystem.

What all of them add over this page's code is the same category of thing: cache-aware C++ histogram kernels, parallel and GPU-resident tree growing, sparse inputs, early stopping, monotonic constraints, and battle-tested handling of edge cases, typically two to four orders of magnitude faster than the Python recursion here. The from-scratch version is enough for understanding, for verifying a paper's claim about a loss function, or for small experiments where a custom (g, h) pair is easier to inject than to configure. Verification is concrete and worth doing in two steps. First, the exact check: on a small dataset run XGBoost with tree_method="exact", max_depth=3, learning_rate=1.0, n_estimators=1, reg_lambda=1.0, subsample=1.0, and base_score=0.5, then dump the single tree with get_booster().get_dump(); your first tree's splits and leaf values should match to within threshold-binning differences, and if you feed both the same pre-binned features they should match nearly exactly, since both compute −G/(H + λ) leaves. Second, the statistical check: with matched depth, learning rate, and tree count, test log loss on held-out data should track XGBoost's within a small margin on data like the example above. When gradient boosting is still the right call over a neural network: heterogeneous features, fewer than a few million rows, latency budgets in microseconds per row on CPU, and teams that need retraining to be a five-minute deterministic job. Neural networks win when there is scale plus structure boosting cannot see: raw text, images, sequences, or embeddings shared across tasks.

Traps and misconceptions

"Boosting means bagging with extra steps." They are opposite mechanisms. Bagging trains independent deep trees in parallel and averages away variance; boosting trains dependent shallow trees sequentially to remove bias, and each tree is meaningless outside the sum. That is why forest trees are deep and boosted trees are shallow, and why a forest is robust to adding more trees while boosting can overfit with them.

"Leaf values are the mean residual." Only for squared loss. For log loss and every other loss with curvature, the least-squares leaf value is wrong; the correct step is the Newton value Σg / Σh (with signs per your convention), and implementations that skip this converge slowly or miscalibrate. This is the single most common bug in from-scratch boosters.

"Tune the number of trees like any hyperparameter." The tree count is a path, not a point: every prefix of the ensemble is a model, so evaluate the staged predictions on a validation set and stop early. Grid-searching n_estimators while refitting from scratch wastes exactly the structure that makes boosting cheap to tune.

"A lower learning rate is just slower." It is a regularizer, not merely a speed knob. Small ν changes which function you converge to, generally a better-generalizing one, because each tree's overfit contribution gets shrunk and later corrected. The standard recipe is to fix ν around 0.05, make the tree budget generous, and let early stopping choose the count.

"The predicted probabilities are trustworthy." Boosted classifiers optimized for accuracy-like objectives are often miscalibrated, and early-stopped log-loss models can be too. If the probabilities feed a decision threshold or an expected-loss computation, as in credit and fraud, check a calibration curve and, if needed, fit isotonic or Platt scaling on held-out data.