37 - Boosting: AdaBoost, Gradient Boosting and a note on XGBoost

You are preparing for an exam with a workbook of exercises. On the first pass you do everything; after that you do not redo the exercises you already know, you concentrate on the ones you got wrong. On the next pass, again: you spend your attention on the errors that remain. After ten passes you have a study method that, exercise after exercise, has corrected its own gaps. Boosting does the same thing: it builds a sequence of simple models in which each one concentrates on the errors of the previous ones.

The difference from the Random Forest is one of principle. The forest trains many deep trees in parallel and independently and averages them to reduce variance. Boosting trains small trees (often of depth 1-4) one after the other and adds them up to reduce bias. In this lesson you derive and implement by hand the two fundamental algorithms: AdaBoost and Gradient Boosting.

Words to clear up right away

• Weak learner: a model only just better than chance. The stump is the smallest possible decision tree: a single question, two leaves (see decision trees).
• Example weight: a number that says how much that example counts in the error function. High weight = “getting it wrong costs more”.
• Residual: the difference y − prediction, that is, how much the model still gets wrong on an example.
• Gradient: for a function of many variables, the vector of partial derivatives; it points in the direction in which the function grows fastest. Going in the opposite direction makes it fall fastest: that is the idea of gradient descent, which you find in Gradient descent.

The four words form a chain: the weak learner is the brick, the weight tells each brick which examples to look at, the residual and the gradient tell you what is left to correct. AdaBoost uses weights; Gradient Boosting uses residuals and gradient. If you remember this division, the rest of the lesson is the same idea written in two ways.

The drawing: at each round the model (the dashed line) gets some examples wrong; in the next round those examples are drawn bigger (higher weight), and the new stump is pushed to classify them correctly. The final result is a weighted sum of all the stumps.

AdaBoost: reweighting the examples

AdaBoost (Adaptive Boosting) works on binary classification with labels y ∈ {−1, +1}. Every example has a weight w_i, at first equal for all (1/n). At each round t:

1. Train a weak learner h_t (typically a stump: a tree with a single test) taking the weights into account.
2. Compute the weighted error: ε_t = Σ_i w_i · [h_t(x_i) ≠ y_i].
3. Compute the learner's weight: α_t = ½ · ln( (1 − ε_t) / ε_t ).
4. Update the weights: w_i ← w_i · exp( −α_t · y_i · h_t(x_i) ), then normalize so that they sum to 1.

The four steps are a cycle, not a list: the error of step 2 decides the weight of step 3, and the weight decides how the data changes for the learner of the next round. No learner knows the others; the only channel of communication between rounds is the example weights.

The final prediction is a weighted vote: H(x) = sign( Σ_t α_t · h_t(x) ).

Why exactly these formulas

The factor y_i·h_t(x_i) is +1 if the learner is right and −1 if it is wrong. So the update multiplies the weight of the correct examples by e^(−α) (less than 1, they go down) and that of the errors by e^(+α) (greater than 1, they go up). With α_t = ½ ln((1−ε)/ε) we have e^α = √((1−ε)/ε). Let us check what happens to the total mass of weight, before normalizing:

mass of the correct examples after the update = (1 − ε) · e^(−α) = (1 − ε) · √(ε/(1−ε)) = √(ε(1−ε))
mass of the wrong examples after the update   = ε · e^(+α)      = ε · √((1−ε)/ε)     = √(ε(1−ε))

The two masses are equal: after the update the total weight of the previous errors is exactly 50%. The next learner faces a problem in which the previous learner is worth no more than a coin, so it is forced to find something new. The weight α_t is large when ε is small (reliable learner), zero when ε = 0.5 (as good as guessing), negative if ε > 0.5 (its vote is inverted).

Implementation by hand

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

class AdaBoostDaZero:
    def __init__(self, n_round=50):
        self.n_round = n_round

    def fit(self, X, y):                       # y in {-1, +1}
        n = len(y)
        w = np.full(n, 1.0 / n)
        self.stump, self.alfa = [], []
        self.storia_pesi = []
        for t in range(self.n_round):
            h = DecisionTreeClassifier(max_depth=1, random_state=0).fit(X, y, sample_weight=w)
            pred = h.predict(X)
            eps = w[pred != y].sum()
            eps = np.clip(eps, 1e-10, 1 - 1e-10)
            alfa = 0.5 * np.log((1 - eps) / eps)
            w = w * np.exp(-alfa * y * pred)
            w = w / w.sum()
            self.stump.append(h); self.alfa.append(alfa); self.storia_pesi.append(w.copy())
        return self

    def decision_function(self, X):
        return sum(a * h.predict(X) for a, h in zip(self.alfa, self.stump))

    def predict(self, X):
        return np.where(self.decision_function(X) >= 0, 1, -1)

X, y01 = make_classification(n_samples=1500, n_features=10, n_informative=5, n_redundant=0,
                             class_sep=1.0, random_state=0)
y = 2 * y01 - 1
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

stump = DecisionTreeClassifier(max_depth=1).fit(Xtr, ytr)
ada = AdaBoostDaZero(n_round=100).fit(Xtr, ytr)
lib = AdaBoostClassifier(n_estimators=100, random_state=0).fit(Xtr, ytr)
print("single stump:       test %.3f" % stump.score(Xte, yte))
print("AdaBoost by hand:   test %.3f" % (ada.predict(Xte) == yte).mean())
print("AdaBoost sklearn:   test %.3f" % lib.score(Xte, yte))
print("alpha of the first 5 rounds:", np.round(ada.alfa[:5], 3))

single stump:       test 0.800
AdaBoost by hand:   test 0.849
AdaBoost sklearn:   test 0.849
alpha of the first 5 rounds: [0.708 0.411 0.258 0.284 0.143]

In the listing, fit initializes equal weights (1/n) and for each round: trains a stump passing sample_weight=w (the library gives more importance to the examples with high weight when computing the impurity); computes the weighted error eps by summing the weights of the mistakes; bounds eps with np.clip to avoid log(0) or divisions by zero if the stump were perfect; derives alfa; updates and normalizes the weights; saves stump, alfa and weights. decision_function is the weighted sum of the votes; predict takes its sign. Note y = 2*y01 - 1: AdaBoost requires −1/+1 labels, because the update formula uses the product y·h(x).

A stump alone scores 0.800 on the test set: stumps are “weak” by construction. A hundred stumps weighted together reach 0.849, identical to the library. Look at the alphas of the first rounds (0.708, 0.411, 0.258, ...): the first learners get the strongest vote and the weight goes down as the problems get harder. The library's weights are double ours (measured: estimator_weights_ gives 1.416, 0.822, 0.516, ... for the first rounds, that is, 2 × 0.708, 2 × 0.411, 2 × 0.258): scikit-learn uses the SAMME variant, which in the binary case scales the vote by a constant factor. The factor changes neither the sign of the weighted vote nor the normalized weights of the examples, and that is why the two accuracies coincide.

The dark side: AdaBoost and wrong labels

The weight update has an important side effect: an example that is got wrong round after round sees its weight grow exponentially. If that example is simply mislabeled, the algorithm spends more and more resources “learning” a mistake. Let us measure it: we flip 15% of the training labels and observe where the weight ends up.

rng = np.random.default_rng(0)
sporche = rng.random(len(ytr)) < 0.15                 # 15% of the training labels flipped
ytr_sporco = np.where(sporche, -ytr, ytr)

ada_s = AdaBoostDaZero(n_round=100).fit(Xtr, ytr_sporco)
w_finale = ada_s.storia_pesi[-1]
print("share of examples with a wrong label: %.1f%%" % (100 * sporche.mean()))
print("share of the total weight they carry after 100 rounds: %.1f%%" % (100 * w_finale[sporche].sum()))
print("test, trained on clean labels: %.3f   on dirty labels: %.3f" % (
    (ada.predict(Xte) == yte).mean(), (ada_s.predict(Xte) == yte).mean()))

share of examples with a wrong label: 13.7%
share of the total weight they carry after 100 rounds: 27.7%
test, trained on clean labels: 0.849   on dirty labels: 0.833

The listing: rng.random(len(ytr)) < 0.15 draws a random number for each example and marks it “dirty” if it is below 0.15: so about 15% (here 13.7%, being a random choice). The dirty labels are flipped with np.where. AdaBoost is retrained and the weight vector of the last round is read: the sum of the weights on the dirty examples says how much attention the algorithm devotes to label errors.

The 13.7% of examples (those actually flipped) end up carrying 27.7% of the total weight: double their share. Accuracy on the test set falls from 0.849 to 0.833. The effect here is moderate; it grows with the share of noise and with the number of rounds. AdaBoost is sensitive to label noise: if your data is dirty, prefer a boosting with more robust losses or with strong regularization, or clean the data first.

Gradient Boosting: gradient descent in function space

Gradient Boosting generalizes the idea to any differentiable loss function. The model is a sum of functions F_M(x) = F₀ + η · Σ_m h_m(x). At each step a new tree h_m is added that moves us in the direction that makes the loss fall fastest: it is a gradient descent in which the “parameters” are the values predicted for each example.

With the squared loss L = ½ (y − F)² the gradient with respect to the prediction F is

∂L/∂F = −(y − F)     →     −∂L/∂F = y − F    (the residual!)

So, for the squared loss, “following the negative gradient” simply means training the new tree on the residuals of the current model. For any other loss, the tree is trained on the negative gradient (called the pseudo-residual). The algorithm:

1. Start from a constant: F₀ = mean of y (for the squared loss it is the optimal constant).
2. For m = 1...M: compute the residuals r_i = y_i − F_{m−1}(x_i); train a small tree h_m on (x_i, r_i); update F_m = F_{m−1} + η · h_m.

The initial constant is the simplest possible model; each following tree shifts the predictions a little toward the residuals, scaled by η. The final model is the constant plus the sum of the trees, so to predict you have to go through all of them: a small M means a lean model, a large M an accurate but slower model at risk of overfitting, as the experiment that follows shows.

The learning rate η (between 0 and 1, often 0.05-0.1) is the shrinkage: each tree contributes only a fraction of what it suggests. Small steps and many trees generalize better than few big steps, at the price of more trees to train. Let us implement it for regression and compare with scikit-learn:

from sklearn.datasets import make_friedman1
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error

Xf, yf = make_friedman1(n_samples=1000, n_features=10, noise=1.0, random_state=0)
Xa, Xb, ya, yb = train_test_split(Xf, yf, test_size=0.3, random_state=0)

class GBRegDaZero:
    def __init__(self, n_alberi=200, eta=0.1, prof=2):
        self.n_alberi, self.eta, self.prof = n_alberi, eta, prof

    def fit(self, X, y):
        self.f0 = y.mean()
        F = np.full(len(y), self.f0)
        self.alberi = []
        for m in range(self.n_alberi):
            residuo = y - F                                     # negative gradient of the squared loss
            h = DecisionTreeRegressor(max_depth=self.prof, random_state=0).fit(X, residuo)
            F = F + self.eta * h.predict(X)
            self.alberi.append(h)
        return self

    def predict(self, X, fino_a=None):
        F = np.full(len(X), self.f0)
        for h in self.alberi[:fino_a]:
            F = F + self.eta * h.predict(X)
        return F

gb = GBRegDaZero(n_alberi=200, eta=0.1, prof=2).fit(Xa, ya)
lib = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1, max_depth=2, random_state=0).fit(Xa, ya)

print("test MSE by hand: %.4f" % mean_squared_error(yb, gb.predict(Xb)))
print("test MSE sklearn: %.4f" % mean_squared_error(yb, lib.predict(Xb)))
d = np.abs(gb.predict(Xb) - lib.predict(Xb))
print("difference between predictions: mean %.3f, max %.3f (scale of y: std %.2f)" % (d.mean(), d.max(), yb.std()))

test MSE by hand: 2.3872
test MSE sklearn: 2.3882
difference between predictions: mean 0.005, max 0.223 (scale of y: std 5.16)

The listing: make_friedman1 generates a classic regression problem with 10 features (only 5 really useful) and a nonlinear relationship. GBRegDaZero.fit starts from f0, the mean of y (the best constant), then for each tree: computes the residuo of the current model, trains a depth-2 DecisionTreeRegressor on it, and updates F by adding eta times the tree's prediction. predict redoes the sum on new data; the fino_a parameter lets you use only the first M trees, handy for tracing the error during construction.

Our loop and GradientBoostingRegressor reach the same error (2.387 against 2.388). The predictions are not identical to the decimal (mean difference 0.005 against a standard deviation of y of 5.16, maximum 0.223; a plausible cause is ties between equivalent splits, not isolated), but they are the same idea: for the squared loss the library does nothing other than this loop on the residuals. Now we notice something that did not happen in the Random Forest.

Boosting can overfit, and has to be stopped

Adding trees to a forest is always harmless. In a boosted model, every added tree reduces the training error, but after a certain point it starts to fit the noise. Let us look at the test error as the trees grow, with a high learning rate:

gb_alto = GBRegDaZero(n_alberi=600, eta=0.5, prof=4).fit(Xa, ya)
for M in (1, 5, 20, 50, 100, 300, 600):
    print("trees %3d   train MSE %.3f   test MSE %.3f" % (
        M, mean_squared_error(ya, gb_alto.predict(Xa, M)), mean_squared_error(yb, gb_alto.predict(Xb, M))))

trees   1   train MSE 11.938   test MSE 13.652
trees   5   train MSE 2.336   test MSE 4.544
trees  20   train MSE 0.587   test MSE 3.803
trees  50   train MSE 0.091   test MSE 3.645
trees 100   train MSE 0.006   test MSE 3.747
trees 300   train MSE 0.000   test MSE 3.768
trees 600   train MSE 0.000   test MSE 3.768

The listing trains a model with a high learning rate (0.5) and deeper trees (4) to make the overfitting visible, and evaluates the error at a growing number of trees by using predict(X, fino_a=M).

The training error falls to zero, while the test error touches its minimum (3.645) around 50 trees and then climbs back (3.747 at 100, 3.768 beyond): it is overfitting in its clearest form, the same you saw with trees and with SVMs. The climb here is small because the function is smooth, but the signal is clear: after the minimum every extra tree worsens generalization. It is controlled with three levers: a lower learning rate (with more trees), shallow trees, and above all early stopping, that is, stopping when the error on a validation set stops improving (in scikit-learn: n_iter_no_change and validation_fraction). The number of trees is a hyperparameter to be chosen in validation, not to be set “high enough”.

From Gradient Boosting to XGBoost: a sketch

XGBoost (Chen and Guestrin, 2016) is an implementation of gradient boosting with a regularized objective and a lot of engineering; the second-order expansion it uses is, as the paper itself says, already present in the literature (Friedman et al.). Here you see the idea; we do not implement it.

• Second-order expansion: it uses not only the gradient g_i of the loss but also the second derivative h_i. For a tree with leaves j, the optimal weight of a leaf is w*_j = − G_j / (H_j + λ), with G_j and H_j the sums of g_i and h_i of the examples that end up in the leaf, and λ an L2 regularization term.
• Regularization in the objective (the gain formula is the one in the original paper): it penalizes the number of leaves (γ per leaf) and the size of the weights (λ). The gain of a split is ½ [ G_L²/(H_L+λ) + G_R²/(H_R+λ) − (G_L+G_R)²/(H_L+H_R+λ) ] − γ; if it is negative, the split is not made.
• Engineering: search for splits on histograms, native handling of missing values, parallelism, subsampling of rows and columns.

The first two items change the mathematics of the tree: the leaf no longer holds the mean of the residuals but a value that also takes the curvature into account, and the useless splits are discarded by the γ penalty. The third item does not change the algorithm, but it is what makes it fast on millions of rows. For a comparison of principle, the Gradient Boosting we have built is enough.

LightGBM and CatBoost follow the same idea with different choices (the LightGBM documentation explains that it grows the tree leaf by leaf, choosing the leaf with the largest loss reduction). If you do not want to install anything else, scikit-learn includes HistGradientBoostingClassifier, of the same histogram-based family:

from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier

Xc, yc = make_classification(n_samples=3000, n_features=20, n_informative=6, flip_y=0.03, random_state=0)
Xc_tr, Xc_te, yc_tr, yc_te = train_test_split(Xc, yc, test_size=0.3, stratify=yc, random_state=0)

hgb = HistGradientBoostingClassifier(early_stopping=True, validation_fraction=0.2,
                                     n_iter_no_change=15, random_state=0).fit(Xc_tr, yc_tr)
rf = RandomForestClassifier(n_estimators=300, random_state=0).fit(Xc_tr, yc_tr)
print("trees used by boosting with early stopping: %d (maximum allowed 100)" % hgb.n_iter_)
print("test boosting: %.3f    test random forest: %.3f" % (hgb.score(Xc_te, yc_te), rf.score(Xc_te, yc_te)))

trees used by boosting with early stopping: 72 (maximum allowed 100)
test boosting: 0.917    test random forest: 0.908

The listing: HistGradientBoostingClassifier is a histogram-based boosting (the features are discretized, by default into at most 255 bins, to speed up the search for thresholds). With early_stopping=True it sets aside 20% of the training set (validation_fraction) and stops when the validation score does not improve for 15 consecutive iterations (n_iter_no_change). n_iter_ reports how many trees were actually used.

The boosting stopped by itself at 72 trees and scores 0.917 against 0.908 for the 300-tree forest: a difference of 0.009 on 900 samples, that is, about 8 examples. There is no universal winner: on tabular data the two families are often close and the difference depends on the dataset and on the care put into the hyperparameter search.

Trade-off: boosting versus Random Forest

• Random Forest: simpler to configure, robust to default values, parallelizable per tree, tolerates label noise better, free OOB estimate.
• Boosting: often more accurate after tuning; sequential training (not parallelizable per tree); more hyperparameters (learning rate, depth, number of trees, regularization); risk of overfitting if it is not stopped.

The two families trade simplicity for potential accuracy: the forest works well almost without tuning, boosting asks you to choose learning rate, depth and number of trees with care, and in exchange it can get higher. The comparison above (0.917 against 0.908 on a single dataset) decides nothing: you need to try both with validation.

Real-world traps

• Number of trees chosen on the test set: use early stopping with a dedicated validation set.
• Dirty labels with AdaBoost: the weights explode on the wrong examples (see above).
• Learning rate and number of trees must be tuned together: halving η requires more trees to obtain the same result.
• Time order ignored: in time series the validation for early stopping must come after the training in time, otherwise it is leakage.
• Comparing models with untuned configurations: a tuned Random Forest against a boosting with default values is not a fair comparison (and vice versa).

The five traps have a common root: boosting adapts as far as you let it. The number of trees, the label noise, the order of time and the tuning are all forms of the same error, that is, letting the model win on the training set without an independent check.

Learning rate and number of trees: two knobs that compensate each other

In Gradient Boosting the overall effect of what you have learned is η · M, the learning rate times the number of trees: a small η needs more trees to reach the same fit. Why then not always use a very small η? Because it costs more trees (more time) and, above all, because the literature (J. Friedman, 2001) and practice show that small steps generalize better: each tree makes a modest correction, so the model is not thrown off balance by single noisy trees. We verify it with staged_predict, a scikit-learn method that returns the predictions after 1, 2, ..., M trees, to obtain the whole test-error curve in a single training. Three learning rates, up to 800 trees of depth 2.

print("learning rate   trees at minimum   minimum MSE   MSE at 50 trees   at 200   at 800")
for eta in (0.5, 0.1, 0.02):
    g = GradientBoostingRegressor(n_estimators=800, learning_rate=eta, max_depth=2, random_state=0).fit(Xa, ya)
    mse = np.array([mean_squared_error(yb, p) for p in g.staged_predict(Xb)])
    print("%-14g  %14d          %10.3f        %8.3f          %6.3f       %6.3f" % (eta, mse.argmin() + 1, mse.min(), mse[49], mse[199], mse[799]))

learning rate   trees at minimum   minimum MSE   MSE at 50 trees   at 200   at 800
0.5                         81               2.346           2.598           2.401        2.510
0.1                        742               2.001           4.202           2.388        2.004
0.02                       800               2.445          14.327           5.288        2.445

For each learning rate the listing trains 800 trees once and, with staged_predict, computes the squared error on the test set after every tree; then it reports the number of trees at which the error is minimum and the minimum, plus three points of the curve. With η = 0.5 the minimum (2.346) is reached already with 81 trees and then rises: the model learns fast and starts to overfit. With η = 0.1, 742 trees are needed and the minimum is better, 2.001. With η = 0.02 the curve is still going down at 800 trees (2.445): many more trees would be needed, and the minimum has not been reached.

Two readings. The first: the combination 0.1 with almost 750 trees beats 0.5 with 81 by almost 15% less error: the small step with many trees is better. The second, practical: the right number of trees depends on the learning rate, so it makes no sense to fix M in advance; you fix η (typically 0.05-0.1), use a high number of trees and leave it to early stopping to stop.

Tree depth and subsampling

The trees in boosting are “weak” on purpose: too deep, and each one fits too much and the corrective mechanism loses its meaning. Depth also controls the order of the interactions between features that the model can represent: trees of depth 1 (stumps) are additive models, in which each feature contributes on its own; with depth 3, interactions among three features are captured. Subsampling (subsample) trains each tree on a random fraction of the examples, introducing the beneficial noise of bagging (Friedman 2002, stochastic gradient boosting). Let us measure the effect of both, with η = 0.1.

print("depth        trees at minimum   minimum MSE   MSE after 300 trees")
for depth in (1, 2, 3, 5):
    g = GradientBoostingRegressor(n_estimators=300, learning_rate=0.1, max_depth=depth, random_state=0).fit(Xa, ya)
    mse = np.array([mean_squared_error(yb, p) for p in g.staged_predict(Xb)])
    print("%-10d  %15d   %11.3f   %17.3f" % (depth, mse.argmin() + 1, mse.min(), mse[-1]))

print("\nsubsample   mean MSE over 5 seeds")
for subsample in (1.0, 0.5, 0.25):
    errors = [mean_squared_error(yb, GradientBoostingRegressor(n_estimators=300, learning_rate=0.1, max_depth=3, subsample=subsample, random_state=s).fit(Xa, ya).predict(Xb)) for s in range(5)]
    print("%-9g   %.3f +- %.3f" % (subsample, np.mean(errors), np.std(errors)))

depth        trees at minimum   minimum MSE   MSE after 300 trees
1                       300         2.892               2.892
2                       297         2.241               2.243
3                       266         1.957               1.974
5                       161         2.525               2.534

subsample   mean MSE over 5 seeds
1           1.974 +- 0.005
0.5         1.933 +- 0.087
0.25        1.926 +- 0.102

The first loop varies the depth with 300 trees and η = 0.1. Stumps (depth 1) reach an error of 2.89 and have not yet stopped improving at 300 trees; depth 2 reaches 2.24, depth 3 the minimum, 1.96, with 266 trees; depth 5 is worse, 2.53, and its minimum comes earlier (161 trees) because each tree learns a lot. The data function (Friedman 1) contains interaction terms between pairs of features and a product, so additive stumps are not enough: this is the reason why a bit of depth helps. Too much, however, brings the overfitting back. The practical rule: depth 3-6 for tabular data, to be chosen with validation.

The second loop varies the subsampling with depth 3. With 1.0 (all the examples) the result is almost deterministic (1.974 ± 0.005: the only randomness is in the ties); with 0.5 and 0.25 the mean error is a little lower (1.933 and 1.926) but the deviation across seeds goes up (0.087 and 0.102): the difference is not significant with five seeds. Subsampling is, here, an option that does no harm and that speeds things up (you train on half of the data), not a demonstrated improvement. A trick from the literature does not always bring a measurable benefit on a given problem.

Gradient Boosting for classification: building it by hand

For regression with squared loss the pseudo-residual was the residual y − F. For binary classification, the model sums the scores F(x) on a logit scale and the probability is p = σ(F). The loss is the log-loss of logistic regression, and its gradient with respect to F is, as you derived there, p − y: so the pseudo-residual (negative gradient) is r = y − p. The construction step is the same: a regression tree is trained on r and F is updated. There is a subtlety (Friedman 2001): instead of using the mean of the residuals as the value of each leaf, a Newton step is taken on the loss and the leaf value becomes γ = Σ r / Σ p(1 − p), the ratio between the gradient and the curvature of the log-loss (p(1 − p) is the second derivative). We implement it all by hand and compare with GradientBoostingClassifier.

from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import log_loss, accuracy_score

X_cl, y_cl = make_classification(n_samples=1500, n_features=10, n_informative=4, flip_y=0.05, class_sep=0.9, random_state=0)
X_cl_tr, X_cl_te, y_cl_tr, y_cl_te = train_test_split(X_cl, y_cl, test_size=0.3, stratify=y_cl, random_state=0)
sigmoid = lambda z: 1 / (1 + np.exp(-z))

class ManualGBClassifier:
    def __init__(self, n_trees=100, eta=0.1, depth=2):
        self.n_trees, self.eta, self.depth = n_trees, eta, depth

    def fit(self, X, y):
        p0 = y.mean()
        self.f0 = np.log(p0 / (1 - p0))                                    # start: the logit of the base rate
        F = np.full(len(y), self.f0)
        self.trees = []
        for _ in range(self.n_trees):
            p = sigmoid(F)
            residual = y - p                                                # pseudo-residual = negative gradient of the log-loss
            tree = DecisionTreeRegressor(max_depth=self.depth, random_state=0).fit(X, residual)
            leaf = tree.apply(X)                                            # which leaf each example falls into
            gamma = {l: residual[leaf == l].sum() / (p[leaf == l] * (1 - p[leaf == l])).sum() for l in np.unique(leaf)}   # Newton step per leaf
            F = F + self.eta * np.array([gamma[l] for l in leaf])
            self.trees.append((tree, gamma))
        return self

    def decision_function(self, X):
        F = np.full(len(X), self.f0)
        for tree, gamma in self.trees:
            F = F + self.eta * np.array([gamma.get(l, 0.0) for l in tree.apply(X)])
        return F

manual = ManualGBClassifier().fit(X_cl_tr, y_cl_tr)
p_manual = sigmoid(manual.decision_function(X_cl_te))
library = GradientBoostingClassifier(n_estimators=100, learning_rate=0.1, max_depth=2, random_state=0).fit(X_cl_tr, y_cl_tr)
p_library = library.predict_proba(X_cl_te)[:, 1]

print("                by hand   scikit-learn")
print("test log-loss    %.4f    %.4f" % (log_loss(y_cl_te, p_manual), log_loss(y_cl_te, p_library)))
print("test accuracy    %.4f    %.4f" % (accuracy_score(y_cl_te, p_manual > 0.5), accuracy_score(y_cl_te, p_library > 0.5)))
print("correlation between the predicted probabilities: %.4f;  maximum difference: %.3f" % (np.corrcoef(p_manual, p_library)[0, 1], np.abs(p_manual - p_library).max()))

                by hand   scikit-learn
test log-loss    0.3471    0.3476
test accuracy    0.8689    0.8689
correlation between the predicted probabilities: 0.9997;  maximum difference: 0.150

The class starts from the logit of the base rate (the best constant), then for each of the 100 rounds computes the current probabilities, the pseudo-residual y − p, trains a depth-2 regression tree on the residuals and, for each leaf, computes the Newton value γ = Σ r / Σ p(1 − p); tree.apply(X) returns the leaf index of each example, and it is what allows the sklearn means to be replaced by the Newton values. The prediction sums the corrections scaled by η and passes through the sigmoid.

The comparison: test log-loss 0.3471 by hand and 0.3476 with scikit-learn, identical accuracy (0.8689), correlation of the probabilities 0.9997. The two implementations do not coincide to the decimal (maximum difference 0.15 on a single example, due to details of tie-breaking between splits and of initialization), but they are the same algorithm: the mean difference is negligible and the ordering of the examples is practically identical. It is the demonstration that boosting for classification is not magic: it is gradient descent in function space, with pseudo-residual y − p, the same “output minus target” quantity as in logistic regression.

Three exercises to try

Exercise 1. Repeat the learning-rate table with η = 0.02 and 3,000 trees: is the minimum now reached, and is it below 2.001? Exercise 2. In the by-hand boosting for classification, replace the Newton step with the mean of the residuals in the leaves: how much worse does the log-loss get at 100 rounds? Exercise 3. Add five noise columns to the Friedman regression and measure whether stumps or depth-3 trees suffer more.

To practice on checking meter readings there is a Python exercise skeleton with three parts to complete: one round of AdaBoost with the check that the weighted error returns to 0.5, the constant-based Gradient Boosting compared with the closed form, and the weight multipliers that explain why a wrongly transcribed reading ruins the ensemble. As long as the points to complete are empty, it prints which exercises are missing.

What you can do now

You can carry out one round of AdaBoost by hand and explain why the weighted error returns to 0.5, you can build a Gradient Boosting for regression and for classification on the pseudo-residuals, and you can recognize when to stop it with validation. You can choose learning rate, depth and subsampling by looking at the error curve, and you can read XGBoost and LightGBM as variants of the same scheme. From the next lesson, clustering, we remove the labels.

Sources