Gradient Boosting — Correcting Errors in Sequence

Descending the Loss, One Tree at a Time

Boosting is the opposite idea to bagging. A forest averages many independent deep trees to cut variance; boosting adds many small, dependent trees in sequence, each fitted to the errors the running ensemble still makes, to cut bias. Friedman's insight is that this is gradient descent in function space: each round fits a tree to the negative gradient of the loss — the ordinary residual under squared loss, $y-\sigma(F)$ under log-loss — and takes a shrunken step in that direction.

The step that is easy to leave out

One detail separates a working implementation from a slow one, and it is easy to miss. The step is not simply the tree's output. Friedman follows the tree fit with a per-leaf line search, choosing the constant in each terminal region that best reduces the loss there. Under squared loss that constant is the leaf mean, so the raw tree can be used directly. Under log-loss it is a Newton step, and using the raw gradient fit instead understates every update — the ensemble still converges, but climbs measurably more slowly.

That difference is worth seeing rather than asserting. Without the per-leaf step the from-scratch booster trailed scikit-learn by 0.016 AUC at 50 trees and was still 0.004 short at 300. With it, the two agree to the fourth decimal at every tree count — and the regression case was already exact, because squared loss does not need the correction. AdaBoost, the historical ancestor, matches exactly too (0.7655 both), which is fitting since it is gradient boosting under exponential loss with stumps.

Learning rate against tree count

The two knobs interact in a way the forest's do not. A low learning rate (0.02) climbs slowly but reaches the best and most stable AUC; $\nu=1.0$ shoots up and then declines as it overfits. Bias falls with every tree added, but past the sweet spot variance creeps back — which is why production practice is a small learning rate with early stopping, and why the tree count must be tuned here when a forest is simply robust to having too many. R's gbm makes the same point with its cross-validation curve, which bottoms out at 356 trees and turns upward after.

Are the probabilities honest?

Every comparison to this point is AUC, which cares only about ranking. A model can rank perfectly and still report probabilities that are systematically wrong — and the folklore says boosting is the offender, because stage-wise fitting of the log-odds keeps pushing confident cases outward until they pile up near 0 and 1. Binning the test set by predicted probability and plotting the observed default rate in each bin, that expectation inverts twice.

Gradient boosting is the best calibrated of the three — expected calibration error 0.012, never off by more than 0.021 in a decile. And the logistic regression is comfortably the worst at 0.059, missing by up to 0.11 — the one model in the table that estimates a probability by maximum likelihood. The classic boosting pathology comes from overfitting, and this fit is regularised into behaving: depth 3, learning rate 0.05, 21,000 rows. The logistic fails in the opposite direction — underfit, unable to bend to the nonlinearity in PAY_1, so its probabilities are wrong at both ends however well it ranks in the middle.

Which sharpens the general point past the folklore: miscalibration is a property of a fit, not of a model class. AUC and calibration measure different things, and the ranking order here is not the calibration order — the forest wins on AUC by 0.002 and loses on ECE. Judge a model on the job it is doing: a sort key for a review queue needs the first column, a probability feeding an expected-loss calculation needs the diagonal, and the logistic would fail that job here despite being the only model designed for it.

Where it lands

Against everything built so far, on the same split: boosting and the forest are neck and neck, and both clear the single tree and the parametric baseline comfortably. The two ensembles arrive at similar flexibility from opposite directions — the forest by averaging away variance, boosting by chipping away at bias — and boosting's fitted curve shows it, bending as finely as the forest but less damped, with sharper wiggles in the sparse tails where its overfitting tendency shows through.

Notebooks

Downloads

Boosting Module — Source Code

"""
gboost.py -- Boosting from scratch: AdaBoost (Freund & Schapire 1997) and Gradient Boosting (Friedman 2001).

Backs the notebooks in  Z-ML Trees-Gradient Boosting.

Boosting is the opposite idea to bagging. A random forest averages many INDEPENDENT deep trees to cut
VARIANCE; boosting adds many SMALL, DEPENDENT trees in sequence, each fit to the errors the running
ensemble still makes, to cut BIAS. The from-scratch core here is the boosting LOOP itself -- the trees
used as base learners are shallow regression trees (built from scratch already in the CART notebook, so
we reuse scikit-learn's for speed; the novelty is the additive, gradient-driven loop).

GRADIENT BOOSTING (Friedman) views the ensemble F(x) as a function optimised by gradient descent in
FUNCTION space. Given a loss L(y, F), each round fits a tree to the negative gradient ("pseudo-residuals")
g_i = -dL/dF, then steps  F <- F + nu * gamma, where nu is the learning rate (shrinkage). For squared loss
the pseudo-residual is just the ordinary residual y - F; for the binary log-loss it is y - sigmoid(F).

The step is not simply the tree's output. Friedman follows the tree fit with a PER-LEAF LINE SEARCH,
choosing the constant in each terminal region that best reduces the loss there. For squared loss that
constant IS the leaf mean, so the tree output can be used directly; for log-loss it is a Newton step

    gamma_L = sum_{i in L} (y_i - p_i) / sum_{i in L} p_i (1 - p_i),

and using the raw gradient-fit instead understates every update -- the ensemble still converges, but
climbs noticeably more slowly. It is the difference between trailing scikit-learn by 0.016 AUC at 50
trees and matching it to the fourth decimal.

ADABOOST (Freund & Schapire) is the historical ancestor: it re-weights the training points after each
weak learner, up-weighting the ones still misclassified, and combines the weak learners by a weighted
vote. It is exactly gradient boosting with the EXPONENTIAL loss and stumps -- the special case that
started the field.
"""

import numpy as np
from sklearn.tree import DecisionTreeRegressor, DecisionTreeClassifier


def _sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -30, 30)))


class GradientBoosting:
    """Friedman (2001) gradient boosting. loss='ls' (squared, regression) or 'log' (binary log-loss).
    Base learner: a depth-limited regression tree fit to the negative gradient."""
    def __init__(self, n_estimators=200, learning_rate=0.1, max_depth=3,
                 subsample=1.0, loss="ls", random_state=0):
        self.n_estimators = n_estimators; self.learning_rate = learning_rate
        self.max_depth = max_depth; self.subsample = subsample; self.loss = loss
        self.random_state = random_state

    def fit(self, X, y):
        X = np.asarray(X, float); y = np.asarray(y, float); n = len(y)
        rng = np.random.default_rng(self.random_state); self.trees = []
        if self.loss == "log":
            p0 = np.clip(y.mean(), 1e-3, 1 - 1e-3); self.F0 = np.log(p0 / (1 - p0))
        else:
            self.F0 = y.mean()
        F = np.full(n, self.F0)
        self.train_loss_ = []
        for m in range(self.n_estimators):
            if self.loss == "log":
                grad = y - _sigmoid(F)                                  # negative gradient of log-loss
            else:
                grad = y - F                                            # negative gradient of squared loss
            if self.subsample < 1.0:
                idx = rng.choice(n, int(self.subsample * n), replace=False)
            else:
                idx = np.arange(n)
            tree = DecisionTreeRegressor(max_depth=self.max_depth, random_state=self.random_state)
            tree.fit(X[idx], grad[idx])
            if self.loss == "log":
                # Friedman's per-leaf line search: a Newton step in each terminal region,
                # gamma_L = sum(y - p) / sum(p(1-p)). For squared loss the leaf means are
                # already optimal, so the raw tree output is used there.
                p = _sigmoid(F); lf = tree.apply(X); lfi = lf[idx]
                num = np.bincount(lfi, weights=grad[idx], minlength=lf.max() + 1)
                den = np.bincount(lfi, weights=(p * (1 - p))[idx], minlength=lf.max() + 1)
                gamma = num / np.maximum(den, 1e-12)
                step = gamma[lf]
            else:
                gamma = None; step = tree.predict(X)
            F += self.learning_rate * step
            self.trees.append((tree, gamma))
            if self.loss == "log":
                p = _sigmoid(F); self.train_loss_.append(-np.mean(y * np.log(p + 1e-9) + (1 - y) * np.log(1 - p + 1e-9)))
            else:
                self.train_loss_.append(np.mean((y - F) ** 2))
        return self

    def staged_decision_function(self, X):
        """running ensemble output after each tree -- for learning curves."""
        X = np.asarray(X, float); F = np.full(X.shape[0], self.F0)
        for tree, gamma in self.trees:
            F += self.learning_rate * (gamma[tree.apply(X)] if gamma is not None else tree.predict(X))
            yield F

    def decision_function(self, X):
        X = np.asarray(X, float); F = np.full(X.shape[0], self.F0)
        for tree, gamma in self.trees:
            F += self.learning_rate * (gamma[tree.apply(X)] if gamma is not None else tree.predict(X))
        return F

    def predict_proba(self, X):
        p = _sigmoid(self.decision_function(X)); return np.column_stack([1 - p, p])

    def predict(self, X):
        if self.loss == "log":
            return (self.decision_function(X) > 0).astype(int)
        return self.decision_function(X)


class AdaBoost:
    """Discrete AdaBoost (Freund & Schapire 1997) with decision stumps, for binary y in {0,1}."""
    def __init__(self, n_estimators=200, random_state=0):
        self.n_estimators = n_estimators; self.random_state = random_state

    def fit(self, X, y):
        X = np.asarray(X, float); yy = 2 * np.asarray(y, int) - 1; n = len(yy)      # {0,1} -> {-1,+1}
        w = np.full(n, 1.0 / n); self.stumps = []; self.alphas = []
        for m in range(self.n_estimators):
            st = DecisionTreeClassifier(max_depth=1, random_state=self.random_state)
            st.fit(X, yy, sample_weight=w)
            pred = st.predict(X)
            err = np.clip(np.sum(w * (pred != yy)) / np.sum(w), 1e-10, 1 - 1e-10)
            alpha = 0.5 * np.log((1 - err) / err)
            w *= np.exp(-alpha * yy * pred); w /= w.sum()
            self.stumps.append(st); self.alphas.append(alpha)
        return self

    def staged_decision_function(self, X):
        X = np.asarray(X, float); F = np.zeros(X.shape[0])
        for a, st in zip(self.alphas, self.stumps):
            F = F + a * st.predict(X); yield F

    def decision_function(self, X):
        X = np.asarray(X, float)
        return np.sum([a * st.predict(X) for a, st in zip(self.alphas, self.stumps)], axis=0)

    def predict(self, X):
        return (self.decision_function(X) > 0).astype(int)

    def predict_proba(self, X):
        p = _sigmoid(2 * self.decision_function(X)); return np.column_stack([1 - p, p])   # calibrated-ish score

References

• Friedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of Statistics 29(5), 1189–1232. — function-space gradient descent, and the per-leaf line search
• Freund, Y. & Schapire, R. E. (1997). A decision-theoretic generalization of on-line learning and an application to boosting. JCSS 55(1), 119–139. — AdaBoost
• Friedman, J., Hastie, T. & Tibshirani, R. (2000). Additive logistic regression: a statistical view of boosting. Annals of Statistics 28(2), 337–407. — AdaBoost as stagewise fitting under exponential loss
• Friedman, J. H. (2002). Stochastic gradient boosting. Computational Statistics & Data Analysis 38(4), 367–378. — the subsampling option
• Ridgeway, G. (2007). Generalized Boosted Models: A Guide to the gbm Package. — the R implementation and its cross-validated tree selection