Part VI — Nonlinear Methods · Lecture 11

Decision Trees, Random Forests, and Gradient Boosting

Leave the hyperplane behind. A tree carves the input space into boxes; a forest averages many de-correlated trees; boosting fits each new tree to what the ensemble still gets wrong. Two of these are still the best thing to run on tabular data.

By the end you can

Derive the optimal leaf value for a fixed partition, and see why finding the best tree is intractable
Derive why misclassification error fails as a splitting criterion and Gini and entropy do not
Derive \(\Var\!\left(\bar f\right)=\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\) and read random forests off it
Derive gradient boosting as gradient descent in function space
Derive the second-order leaf value \(w^\star=-G/(H+\lambda)\) and the split-gain formula behind XGBoost
Say precisely which of bias and variance each ensemble method attacks

Assumed

Phase B, Lecture 7: bias–variance decomposition, cross-validation
Phase B, Lecture 3: gradient descent, step sizes, why small steps can beat large ones
Phase A, Lecture 13: variance of a sum, covariance

1A completely different way to be nonlinear

Everything in Part II fitted a function of the form \(f(\vx)=\vw^\top\phi(\vx)\): you choose the basis \(\phi\) up front and the model chooses the weights. Trees abandon that entirely. A tree makes no smoothness assumption, has no coefficients, and is invariant to any monotone transformation of any individual feature — take logs of half your columns and the tree is unchanged.

Intuition

A tree is a game of twenty questions. Each internal node asks one question of the form "is \(x_j \le t\)?", and after a few questions you land in a leaf and answer with whatever the training points in that leaf averaged to.

Geometrically: each question is an axis-aligned cut, and the questions together carve \(\RR^d\) into rectangular boxes. The fitted function is constant on each box. That is the entire model class — piecewise constant on axis-aligned boxes.

Regression tree

A partition of the input space into \(M\) disjoint regions \(R_1,\dots,R_M\) formed by axis-aligned splits, together with a value per region:

\[ f(\vx) \;=\; \sum_{m=1}^{M} c_m \,\mathbb{1}[\vx \in R_m]. \]

For classification the leaf holds a class distribution \(\hat p_{mk}\) instead of a single number, and predicts \(\argmax_k \hat p_{mk}\).

Given the boxes, the leaf values are immediate0 / 4

\[\hat R(f)=\sum_{i=1}^{n}\big(y_i-f(\vx_i)\big)^2 = \sum_{m=1}^{M}\ \sum_{i:\,\vx_i\in R_m}\big(y_i-c_m\big)^2\]

The objective splits over regions. Because the regions are disjoint and cover everything, each training point contributes to exactly one inner sum. So the \(M\) leaf values are \(M\) completely independent one-parameter problems — no coupling at all, unlike the normal equations where every coefficient depends on every other.

\[\frac{\partial}{\partial c_m}\sum_{i\in R_m}(y_i-c_m)^2 = -2\sum_{i\in R_m}(y_i-c_m)=0 \;\Longrightarrow\; \hat c_m = \frac{1}{|R_m|}\sum_{i\in R_m} y_i\]

Each leaf is the mean. Phase A Lecture 13's result again: the constant minimising squared deviation is the mean. The second derivative is \(2|R_m|>0\), so it is a minimum.

\[\text{0–1 loss} \Rightarrow \hat c_m = \text{majority class}, \qquad \text{absolute loss} \Rightarrow \hat c_m = \text{median}\]

The leaf value is decided by the loss, not by the tree. Exactly the correspondence from Lecture 4: squared loss → mean, absolute loss → median, 0–1 loss → mode. Change the loss and only this step changes; the partition machinery is untouched.

\[\text{but choosing } \{R_m\} \text{ optimally is NP-hard}\]

The hard half. The number of distinct binary trees over \(n\) points with \(d\) features is astronomically large, and finding the partition that minimises \(\hat R\) is NP-hard (Hyafil & Rivest, 1976). So every tree algorithm in existence is greedy: choose the single best split now, recurse, never reconsider. Nothing guarantees the result is close to optimal — it is a heuristic, and calling it one is more honest than most textbooks are.

2Choosing a split: impurity, and a criterion that quietly fails

For regression, "best split" means the one that most reduces the sum of squared errors, and the code below shows the \(O(n\log n)\) prefix-sum trick that makes it fast. For classification we need a measure of how mixed a node is. Three candidates, all functions of the class proportions \(p_k\) in a node:

Impurity measures

\[ \text{misclassification: } \ I_{\text{err}} = 1-\max_k p_k, \qquad \text{Gini: } \ I_G = \sum_k p_k(1-p_k), \qquad \text{entropy: } \ I_H = -\sum_k p_k\log_2 p_k. \]

All three are zero for a pure node and maximal at the uniform distribution. The gain of a split is the parent impurity minus the size-weighted average of the children's:

\[ \Delta = I(\text{parent}) - \frac{n_L}{n}I(\text{left}) - \frac{n_R}{n}I(\text{right}). \]

Why misclassification error is the wrong criterion0 / 5

\[\Delta \;=\; I(\bar p) - \Big[\tfrac{n_L}{n}I(p_L)+\tfrac{n_R}{n}I(p_R)\Big], \qquad \bar p=\tfrac{n_L}{n}p_L+\tfrac{n_R}{n}p_R\]

Set up the geometry. First check the identity \(\bar p = \frac{n_L}{n}p_L+\frac{n_R}{n}p_R\): the parent's class-1 proportion is the size-weighted average of the children's, because the counts simply add. So the gain compares the height of \(I\) at a weighted average against the weighted average of the heights.

\[\Delta \ge 0 \ \text{ for all splits} \iff I \text{ is concave}\]

That is exactly Jensen's inequality. \(I(\text{weighted average}) \ge \text{weighted average of } I\) is the definition of concavity. So any concave impurity gives non-negative gain, and the gain is the vertical gap between the curve and the chord joining \((p_L,I(p_L))\) to \((p_R,I(p_R))\). Draw that chord in the figure below — it is the whole story in one picture.

\[I_{\text{err}}(p)=\min(p,1-p) \ \text{ is concave but } \textbf{piecewise linear}\]

Here is the defect. Concave, so the gain is never negative — but only weakly concave. On either side of \(p=\tfrac12\) it is a straight line, and a chord between two points on a straight line lies exactly on it. Gain zero.

\[p_L,p_R \ \text{ on the same side of } \tfrac12 \ \Longrightarrow\ \Delta_{\text{err}}=0\]

The consequence. Any split whose two children still have the same majority class scores exactly zero under misclassification error, no matter how much it purifies them. The greedy search sees a flat landscape and stops. Concretely: a parent with 400/400 split into (300,100) and (100,300) and split into (200,0) and (200,400) both score \(\Delta_{\text{err}}=0.25\) — indistinguishable, although the second produces a perfectly pure leaf.

\[I_G(p)=2p(1-p) \ \text{ and } \ I_H(p) \ \text{ are } \textbf{strictly} \text{ concave}\]

The fix. \(I_G''(p)=-4<0\) everywhere and \(I_H''(p)=-\frac{1}{\ln 2}\big(\frac1p+\frac1{1-p}\big)<0\) everywhere. Strict concavity means the chord is strictly below the curve whenever \(p_L\neq p_R\), so every non-trivial split gets strictly positive credit and the greedy search always has a gradient to follow. This is why sklearn offers gini and entropy and not error — and why you still use misclassification error for pruning, where you are comparing whole trees rather than searching for a split.

Impurity as a concave function — and the gain as the gap to the chordInteractive

Gini 2p(1−p)entropy (÷2)misclassification min(p,1−p)chord = weighted child impurity

class-1 proportion in the LEFT child0.20

class-1 proportion in the RIGHT child0.44

fraction of points going left0.50

show chord for: 0 Gini · 1 entropy · 2 misclassification2

The gain of a split is the vertical distance from the curve down to the chord, measured at the parent's proportion p̄. Set both children to the same side of ½ and switch the chord to misclassification: the chord lies flat on the curve and the gain is exactly zero, while Gini and entropy still register a real improvement. That is the whole argument for strict concavity.

3Growing, and then not growing too far

The recursion is: find the best split by scanning every feature and every candidate threshold, apply it, recurse on both children, stop at a depth limit or a minimum leaf size. A fully grown tree with \(n\) distinct points has \(n\) leaves and interpolates the training data exactly — training error zero, test error terrible. Depth is the capacity knob, and Lecture 7's U-shaped test curve appears in its cleanest form.

Practical engineering note

Pre-pruning vs post-pruning. Stopping early on a depth or minimum-leaf-size rule (pre-pruning) is fast but short-sighted: a split that looks useless can enable an excellent split beneath it — the XOR pattern is the standard example, where neither first split reduces impurity at all. Cost-complexity pruning (post-pruning) grows the tree fully and then prunes back to minimise \(\hat R(T)+\alpha|T|\) over subtrees, choosing \(\alpha\) by cross-validation. The penalty \(\alpha|T|\) is the same regularisation idea as Lecture 6, with the number of leaves playing the role of the norm.

In practice, when a tree is going into an ensemble, people mostly skip pruning: the ensemble handles the variance. Pruning matters when you are shipping the single tree because someone wants to read it.

Grow a tree and watch it overfitInteractive

class 0class 1learned boundary (axis-aligned boxes)

maximum depth2

label noise0.10

minimum points per leaf1

rotate the true boundary (degrees)0

Every region is a rectangle — that is not a drawing convenience, it is the model class. Push the depth up and the shaded regions fragment into slivers that chase individual noisy points while the test error climbs. Then rotate the true boundary off the axes and watch the tree approximate a diagonal line with a staircase: the one thing an axis-aligned model cannot do gracefully.

Common trap

Trees are unstable, and this is their defining flaw. Change one training point and the top split can change; if it does, the entire tree below it is rebuilt and the fitted function can be unrecognisable. In Lecture 7's terms: very low bias (a deep tree can fit anything), very high variance.

Everything in the rest of this lecture is a response to that one sentence. Bagging averages the variance away. Boosting takes the opposite route — deliberately use high-bias, low-variance stumps and remove the bias by addition. Two opposite fixes for two opposite starting points, which is why you must know which problem you have before choosing.

4Bagging: averaging away the variance

If \(B\) independent estimators each have variance \(\sigma^2\), their average has variance \(\sigma^2/B\). We do not have \(B\) independent training sets, so we fake them: draw \(B\) bootstrap resamples of the training data (sample \(n\) rows with replacement), fit a tree to each, average the predictions. That is bagging — bootstrap aggregating.

But the trees are not independent — they see overlapping data. The next derivation is the single most important formula in ensemble learning, because it tells you exactly what that dependence costs.

The variance of an average of correlated estimators0 / 6

\[\Var\Big(\frac1B\sum_{b=1}^{B} f_b\Big) = \frac{1}{B^2}\Big[\sum_{b}\Var(f_b) + \sum_{b\ne b'}\Cov(f_b,f_{b'})\Big]\]

Expand the variance of a sum. Phase A Lecture 13: \(\Var(\sum_b f_b)=\sum_b\Var(f_b)+\sum_{b\ne b'}\Cov(f_b,f_{b'})\), and the \(1/B\) outside squares to \(1/B^2\). The cross terms are the whole point — dropping them is the mistake that makes people expect \(\sigma^2/B\).

\[\Var(f_b)=\sigma^2 \ \ \forall b, \qquad \Cov(f_b,f_{b'})=\rho\sigma^2 \ \ \forall b\ne b'\]

The exchangeability assumption. The trees are built by an identical random procedure on the same training set, so they are exchangeable: same marginal variance \(\sigma^2\), same pairwise correlation \(\rho\). Reasonable, and the code below checks it holds well enough to predict the answer to a few percent.

\[= \frac{1}{B^2}\Big[B\sigma^2 + B(B-1)\rho\sigma^2\Big]\]

Count the terms. There are \(B\) diagonal terms and \(B(B-1)\) off-diagonal ordered pairs. The cross terms outnumber the diagonal ones by a factor of \(B-1\) — which is why they dominate for large \(B\).

\[\boxed{\ \Var\Big(\frac1B\sum_b f_b\Big) \;=\; \rho\sigma^2 \;+\; \frac{1-\rho}{B}\,\sigma^2\ }\]

Simplify. \(\frac{1}{B}\sigma^2+\frac{B-1}{B}\rho\sigma^2 = \rho\sigma^2+\frac{1-\rho}{B}\sigma^2\). Read the two terms: the second vanishes as \(B\to\infty\), the first does not.

\[\lim_{B\to\infty}\Var = \rho\sigma^2 \;>\;0 \quad\text{whenever } \rho>0\]

The floor. No number of trees ever gets you below \(\rho\sigma^2\). Adding trees past the point where \(\frac{1-\rho}{B}\sigma^2 \ll \rho\sigma^2\) — that is, past \(B\approx\frac{1-\rho}{\rho}\) — is wasted computation. This is why \"just use more trees\" stops helping, and it also tells you the only two ways forward: lower \(\rho\), or lower \(\sigma^2\).

\[\text{random forest: lower } \rho \text{ by choosing splits from a random subset of } m<d \text{ features}\]

And that is a random forest. Bagging alone leaves \(\rho\) high because every tree finds the same dominant feature and splits on it first. Restricting each node to a random subset of \(m\) features forces different trees down different paths, cutting \(\rho\). But it also makes each individual tree worse, raising \(\sigma^2\) — so \(m\) is a hyperparameter with an interior optimum, and the code below measures both effects separately instead of taking the usual \(m=\sqrt d\) on faith.

The variance floor: why more trees stop helpingInteractive

ρ = 0.02ρ = 0.10ρ = 0.30selected ρ

correlation between trees ρ0.10

number of trees B50

Vertical axis is variance as a multiple of one tree's variance σ². The dashed line is the floor ρσ². Drag ρ toward zero and the curve keeps falling toward 1/B; push it up and the curve flattens almost immediately. The readout says how many trees actually buy you anything — past that point you are paying compute for nothing, and the only remaining lever is de-correlation.

Where this shows up in ML

Out-of-bag error, free of charge. A bootstrap sample of size \(n\) omits each particular row with probability \((1-1/n)^n \to e^{-1}\approx 0.368\). So every training row is out-of-bag for about 37% of the trees, and you can predict it using exactly those trees. Averaging the error over all rows gives an almost-unbiased estimate of test error with no held-out set and no cross-validation refits. Lecture 7 spent a whole section on how to buy an honest error estimate; the bootstrap hands one over as a by-product.

The caveat: OOB uses an ensemble of about \(0.37B\) trees rather than \(B\), so it is very slightly pessimistic. With \(B\) in the hundreds the difference is negligible.

Check yourself

You bag 500 fully-grown decision trees and the test error barely improves over a single tree. Which explanation is most likely, and what should you try?

Not quite. The formula says variance is \(\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\). At \(B=500\) the second term is already tiny for any realistic \(\rho\), so you are sitting on the floor. Ten times more trees changes essentially nothing — and the figure above lets you confirm that by eye.

Correct. This is the diagnosis the formula forces. High \(\rho\) usually means one feature dominates every split, so every tree starts the same way. Random feature selection at each node is the standard cure and is precisely what turns bagging into a random forest.

Not quite. Fully-grown trees are already at maximum depth — they interpolate the training data. The problem is variance, not bias, and deeper is not available anyway.

Not quite. Possible, and worth checking — but implausible for fully-grown trees, which have almost no bias. Diagnose before switching: measure \(\rho\) as the code block does. If \(\rho\) is high, de-correlate; if the single tree's training error is also high, then and only then is bias the story.

5Boosting: fitting what the ensemble still gets wrong

Bagging builds its trees in parallel and independently. Boosting builds them sequentially, each one trained on the mistakes of everything built so far. The models are deliberately weak — depth-1 or depth-3 stumps with enormous bias — and the ensemble removes that bias by accumulation. Same ingredient, opposite recipe, opposite failure mode.

Gradient boosting is gradient descent in function space0 / 6

\[\min_{F}\ \hat R(F) = \sum_{i=1}^{n} L\big(y_i, F(\vx_i)\big)\]

Change the variable we optimise over. Until now we minimised over a parameter vector \(\vw\). Here the variable is the function \(F\). But \(\hat R\) only ever touches \(F\) through the \(n\) numbers \(F(\vx_1),\dots,F(\vx_n)\), so for the purpose of differentiating we can treat \(\vF=(F(\vx_1),\dots,F(\vx_n))\in\RR^n\) as an ordinary vector. That reduction is what makes the whole idea work.

\[g_i \;=\; \frac{\partial L(y_i,F(\vx_i))}{\partial F(\vx_i)}, \qquad \vF \leftarrow \vF - \nu\,\vg\]

Take one gradient-descent step in \(\RR^n\), with step size \(\nu\). Exactly Lecture 3, with the \(n\) function values as the parameters. For squared loss \(L=\tfrac12(y-F)^2\) we get \(g_i=F(\vx_i)-y_i\), so \(-g_i=y_i-F(\vx_i)\) is the residual. The famous \"fit the residual\" rule is not a heuristic — it is the negative gradient of squared loss, and nothing more.

\[\text{problem: } -\vg \ \text{is defined only at the } n \text{ training points}\]

The obstacle. A step in \(\RR^n\) tells us how to change \(F\) at the training inputs and says nothing about anywhere else. But we need a function we can evaluate on new data. Descent in \(\RR^n\) alone would just memorise the training set.

\[h_m \;=\; \argmin_{h\in\mathcal{H}} \sum_{i=1}^{n}\big(-g_i - h(\vx_i)\big)^2\]

Project the gradient onto the model class. Fit a small tree to the negative gradients by least squares — find the member of \(\mathcal H\) that best approximates the step we wanted to take. The tree generalises the direction to the whole input space. Note this is a least-squares fit regardless of what \(L\) is: the loss enters only through the targets \(-g_i\), which is exactly why one implementation handles every differentiable loss.

\[F_m(\vx) \;=\; F_{m-1}(\vx) \;+\; \nu\, h_m(\vx), \qquad F_0 = \argmin_c \textstyle\sum_i L(y_i,c)\]

Take the step. \(\nu\in(0,1]\) is the learning rate, called shrinkage here. It is Lecture 3's step size, and it plays the same role: small steps in a direction that is only approximately right end up ahead of large ones. Empirically \(\nu\approx0.1\) with many rounds beats \(\nu=1\) with few, and the code block measures the gap.

\[F_M = F_0 + \nu\sum_{m=1}^{M} h_m \quad\text{— an additive model, fitted greedily, stage by stage}\]

What we have built. \(M\) and \(\nu\) trade off against each other (roughly \(\nu M \approx\) constant), so tune one and fix the other. Unlike bagging, more rounds eventually overfit: the ensemble keeps chasing residuals that are pure noise, so \(M\) must be chosen by a validation set — early stopping is not optional here, it is the main regulariser.

Gradient boosting on one dimension: watch the residual get eatenInteractive

training pointstrue functionboosted fit F_Mresiduals y − F_M

boosting rounds M0

learning rate ν0.30

depth of each tree1

noise level0.25

At M = 0 the model is a horizontal line at ȳ and the residuals are the data itself. Each round fits one small tree to the residuals and adds ν times it. Watch the residuals (lower band) collapse toward zero. Set ν = 1 and it lurches; set ν = 0.05 and it crawls but ends up smoother. Push M to 150 with ν = 1 and the fit starts tracking individual noise points — that is overfitting, and no amount of extra rounds undoes it.

6The second-order version: where XGBoost comes from

Gradient boosting takes a first-order step. If we keep the second-order term as well, we get a Newton step in function space — and out of that falls a closed form for the leaf values and a principled split criterion. This derivation is the core of XGBoost, LightGBM and CatBoost, and it is short.

Newton boosting: the leaf value and the gain formula0 / 6

\[\hat R(F_{m-1}+h) \approx \sum_{i=1}^{n}\Big[L(y_i,F_{m-1}(\vx_i)) + g_i h(\vx_i) + \tfrac12 h_i\, h(\vx_i)^2\Big] + \Omega(h)\]

Second-order Taylor expansion in the function value, Phase A Lecture 6 applied per point, with \(g_i=\partial_F L\) and \(h_i=\partial^2_F L\) evaluated at the current \(F_{m-1}(\vx_i)\). \(\Omega\) is a regulariser on the new tree. The constant first term does not involve \(h\), so drop it.

\[\Omega(h)=\gamma T + \tfrac{\lambda}{2}\sum_{j=1}^{T} w_j^2, \qquad h(\vx)=\sum_{j=1}^{T} w_j\,\mathbb{1}[\vx\in R_j]\]

Say what the tree is and how we penalise it. \(T\) leaves, value \(w_j\) in leaf \(j\); \(\gamma\) charges per leaf (complexity control, as in cost-complexity pruning) and \(\lambda\) is a ridge penalty on the leaf values (Lecture 6, now applied to leaves instead of coefficients).

\[\text{obj} = \sum_{j=1}^{T}\Big[G_j w_j + \tfrac12 (H_j+\lambda) w_j^2\Big] + \gamma T, \qquad G_j=\!\!\sum_{i\in R_j}\!g_i,\ \ H_j=\!\!\sum_{i\in R_j}\!h_i\]

Regroup by leaf. Every training point lies in exactly one leaf, so \(\sum_i\) becomes \(\sum_j\sum_{i\in R_j}\), and inside a leaf \(h(\vx_i)=w_j\) is constant and comes out of the sum. The \(T\) leaves are now decoupled — the same decoupling as in the very first derivation of this lecture, and the reason everything closes in one line.

\[\frac{\partial}{\partial w_j}\Big[G_j w_j+\tfrac12(H_j+\lambda)w_j^2\Big]=G_j+(H_j+\lambda)w_j=0 \;\Longrightarrow\; \boxed{\,w_j^\star=-\frac{G_j}{H_j+\lambda}\,}\]

One-dimensional quadratic, minimised exactly. Second derivative \(H_j+\lambda>0\) for a convex loss, so it is a minimum. Sanity-check with squared loss: \(g_i=F-y_i\), \(h_i=1\), so \(G_j=-\sum_{i\in R_j}(y_i-F_i)\) and \(H_j=|R_j|\), giving \(w_j^\star=\frac{\sum(y_i-F_i)}{|R_j|+\lambda}\) — the mean residual, shrunk. At \(\lambda=0\) it is exactly the mean residual, so this reduces to ordinary gradient boosting, as it must. The code block verifies this to machine precision.

\[\text{obj}^\star = -\frac12\sum_{j=1}^{T}\frac{G_j^2}{H_j+\lambda} + \gamma T\]

Substitute \(w^\star\) back. \(G w^\star+\tfrac12(H+\lambda)w^{\star2} = -\frac{G^2}{H+\lambda}+\frac12\frac{G^2}{H+\lambda} = -\frac12\frac{G^2}{H+\lambda}\). This is a score for a tree structure with the leaf values already optimised out — smaller is better.

\[\text{Gain} = \frac12\left[\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{(G_L{+}G_R)^2}{H_L{+}H_R+\lambda}\right]-\gamma\]

The split criterion, derived rather than chosen. Difference of the structure scores before and after splitting one leaf. Two things to notice. First, \(-\gamma\): a split must improve the score by at least \(\gamma\) or it is rejected — pruning built into the search. Second, at \(\lambda=\gamma=0\) with squared loss this reduces to \(\frac{(\sum_L r)^2}{n_L}+\frac{(\sum_R r)^2}{n_R}-\frac{(\sum r)^2}{n}\), which is precisely the prefix-sum criterion coded in best_split below. Everything is consistent.

Practical engineering note

What to actually do on a tabular problem in 2026.

Start with gradient-boosted trees (XGBoost / LightGBM / CatBoost). On tabular data with mixed types and a few thousand to a few million rows, they still beat neural networks in most careful benchmarks. This is not nostalgia — it is the one large domain where deep learning did not take over.
Tune \(\nu\) and \(M\) together, with early stopping. Fix \(\nu=0.05\)–\(0.1\), set \(M\) large, and stop when validation error stops improving. Then tune depth (3–8), min_child_weight (that is \(H_j\), the leaf's total Hessian) and subsampling.
Random forests when you want robustness with no tuning. Defaults work, it is trivially parallel, and OOB gives a free error estimate. It will rarely be the best model and rarely be a bad one.
Distrust impurity-based feature importances. They are biased toward high-cardinality and continuous features, which get more candidate split points and so more chances to look good on noise. Use permutation importance on held-out data, or SHAP values, and remember both still describe the model, not the world.
Trees cannot extrapolate. Every prediction is an average of training labels, so a tree's output is bounded by \([\min y_i, \max y_i]\) forever. On a time series with a trend, a tree will predict a flat line beyond the training range. This is a hard structural limit, not a tuning problem.

7Build all of it

CART, bagging, the variance formula, and boosting — from scratchPython · NumPy

import numpy as np, time
np.set_printoptions(precision=3, suppress=True)
_t0 = time.time()

# ==================== a CART regression tree, from scratch ====================
class Node:
    __slots__ = ("j", "t", "lo", "hi", "val")
    def __init__(self): self.j = -1; self.t = 0.0; self.lo = self.hi = None; self.val = 0.0

def best_split(X, y, min_leaf):
    """O(n log n) per feature: sort once, then prefix sums.
       Minimising SSE == maximising  (sum_L)^2/n_L + (sum_R)^2/n_R  (see the gain
       formula derivation with lambda = gamma = 0)."""
    n, d = X.shape
    if n < 2*min_leaf: return None
    best_score, bj, bt = y.sum()**2/n, -1, 0.0        # score of NOT splitting
    for j in range(d):
        o = np.argsort(X[:, j], kind="stable")
        xs, ys = X[o, j], y[o]
        cs = np.cumsum(ys); tot = cs[-1]
        k  = np.arange(1, n)                           # size of the left child
        ok = (xs[1:] > xs[:-1]) & (k >= min_leaf) & (n - k >= min_leaf)
        if not ok.any(): continue
        score = np.where(ok, cs[:-1]**2/k + (tot - cs[:-1])**2/(n - k), -np.inf)
        i = int(score.argmax())
        if score[i] > best_score:
            best_score, bj, bt = score[i], j, 0.5*(xs[i] + xs[i+1])
    return None if bj < 0 else (bj, bt)

def grow(X, y, depth, min_leaf=1, mtry=None, rg=None):
    nd = Node(); nd.val = float(y.mean())              # leaf value = mean (derived above)
    if depth == 0 or len(y) < 2*min_leaf or np.all(y == y[0]): return nd
    cols = np.arange(X.shape[1]) if mtry is None else rg.choice(X.shape[1], mtry, replace=False)
    s = best_split(X[:, cols], y, min_leaf)
    if s is None: return nd
    j = int(cols[s[0]]); t = s[1]; m = X[:, j] <= t
    if m.sum() == 0 or (~m).sum() == 0: return nd
    nd.j, nd.t = j, t
    nd.lo = grow(X[m],  y[m],  depth-1, min_leaf, mtry, rg)
    nd.hi = grow(X[~m], y[~m], depth-1, min_leaf, mtry, rg)
    return nd

def pred(nd, X):
    out = np.empty(len(X)); stack = [(nd, np.arange(len(X)))]
    while stack:
        node, ii = stack.pop()
        if node.j < 0: out[ii] = node.val; continue
        m = X[ii, node.j] <= node.t
        stack.append((node.lo, ii[m])); stack.append((node.hi, ii[~m]))
    return out

def nleaves(nd): return 1 if nd.j < 0 else nleaves(nd.lo) + nleaves(nd.hi)

def f_true(X): return np.sin(3*X[:,0]) + 0.6*X[:,1]**2 - 0.5*X[:,0]*X[:,1]
rg  = np.random.default_rng(0)
Xtr = rg.uniform(-2, 2, (400, 2));  ytr = f_true(Xtr) + rg.normal(0, 0.3, 400)
Xte = rg.uniform(-2, 2, (1500, 2)); yte = f_true(Xte) + rg.normal(0, 0.3, 1500)

print("--- 1. ONE TREE: depth is the capacity knob ---")
print("  depth   leaves   train MSE   test MSE")
for dep in (1, 2, 3, 5, 8, 12, 20):
    T = grow(Xtr, ytr, dep)
    print(f"  {dep:5d}   {nleaves(T):6d}   {np.mean((pred(T,Xtr)-ytr)**2):9.4f}  "
          f"{np.mean((pred(T,Xte)-yte)**2):9.4f}")
print("  Irreducible noise = 0.09. Training MSE goes to 0; test MSE bottoms out")
print("  around depth 8 and then rises. Lecture 7's U-curve, with depth as capacity.")

# ============ 2. the impurity criteria, on the classic 400/400 example ========
def gini(p): return 2*p*(1-p)
def ent(p):
    p = np.clip(p, 1e-12, 1-1e-12); return -(p*np.log2(p) + (1-p)*np.log2(1-p))
def mis(p):  return np.minimum(p, 1-p)

print()
print("--- 2. WHY MISCLASSIFICATION ERROR CANNOT CHOOSE A SPLIT ---")
print("  parent: 400 of class 1, 400 of class 0   (impurity at p = 0.5)")
for name, (a, b, c, e) in (("split A", (300, 100, 100, 300)),
                           ("split B", (200,   0, 200, 400))):
    nl, nr = a+b, c+e
    row = f"  {name}  left {a}/{nl}, right {c}/{nr}   "
    for f, lab in ((mis, "misclass"), (gini, "Gini"), (ent, "entropy")):
        g = f(0.5) - (nl/800)*f(a/nl) - (nr/800)*f(c/nr)
        row += f"{lab} {g:.4f}   "
    print(row)
print("  Identical misclassification gain, although split B produces a PURE leaf.")
print("  Gini and entropy both prefer B. Strict concavity is doing that work.")

# ==== 3. bagging: verify  Var(average) = rho*sigma^2 + (1-rho)*sigma^2/B ======
# The correlation rho comes from the SHARED TRAINING SET, not from the bootstrap.
# Two trees inside one forest are conditionally independent given the data, so
# you cannot see rho at all without several independent training sets.
print()
print("--- 3. THE CORRELATED-AVERAGE VARIANCE FORMULA, VERIFIED ---")
D = 10
def g_true(X): return 2.0*np.sign(X[:,0]) + np.sin(3*X[:,1]) + 0.8*X[:,2]**2 + 0.7*X[:,3]
Xq  = np.random.default_rng(77).uniform(-2, 2, (250, D))
yq  = g_true(Xq) + np.random.default_rng(78).normal(0, 0.4, 250)
R, B = 10, 10

def forest_study(mtry):
    Fb = np.zeros((R, B, len(Xq))); rgd = np.random.default_rng(0)
    for r in range(R):                                   # R independent training sets
        Xa = rgd.uniform(-2, 2, (200, D)); ya = g_true(Xa) + rgd.normal(0, 0.4, 200)
        for b in range(B):                               # B bootstrap trees on each
            idx = rgd.integers(0, 200, 200)
            Fb[r, b] = pred(grow(Xa[idx], ya[idx], 6, 3, mtry, rgd), Xq)
    s2 = float(Fb.reshape(R*B, -1).var(0, ddof=1).mean())
    def paircorr(i, j):                                  # rho, from DISTINCT tree pairs
        a = Fb[:, i, :] - Fb[:, i, :].mean(0)
        b = Fb[:, j, :] - Fb[:, j, :].mean(0)
        return ((a*b).mean(0)/np.sqrt((a*a).mean(0)*(b*b).mean(0) + 1e-30)).mean()
    rho = float(np.mean([paircorr(i, j) for i in range(4) for j in range(4) if i < j]))
    return s2, rho, Fb

study = {m: forest_study(m) for m in (None, 5, 2)}
s2, rho, Fb = study[None]
print(f"  sigma^2 (variance of ONE tree's prediction)            = {s2:.4f}")
print(f"  rho     (correlation of two trees in the same forest)  = {rho:.4f}")
print("     B    formula rho*s2 + (1-rho)*s2/B      measured variance")
for Bs in (1, 2, 4, 8, 10):
    print(f"  {Bs:4d}              {rho*s2 + (1-rho)*s2/Bs:8.4f}                {float(Fb[:, :Bs, :].mean(1).var(0, ddof=1).mean()):8.4f}")
print("  TWO estimated numbers predict the whole curve. That is the derivation working.")
print(f"  Floor as B -> infinity: rho*sigma^2 = {rho*s2:.4f}. At B=10 the variance is")
print(f"  {rho*s2 + (1-rho)*s2/10:.4f}, still {100*(rho*s2+(1-rho)*s2/10)/(rho*s2)-100:.0f}% ABOVE the floor, so here")
print(f"  more trees do still help. They stop helping past B = (1-rho)/(0.1*rho) = "
      f"{int(np.ceil((1-rho)/(0.1*rho)))}, which is where the 1/B term drops to a tenth of the floor.")

print()
print("--- 4. FEATURE SUBSAMPLING LOWERS rho -- AND RAISES sigma^2 ---")
print("  mtry   sigma^2     rho    floor rho*s2   bias^2   forest MSE (B=10)")
for m in (None, 5, 2):
    s2m, rhom, Fm = study[m]
    bias2 = float(((Fm.reshape(R*B, -1).mean(0) - g_true(Xq))**2).mean())
    print(f"  {str(m):4s}   {s2m:7.3f}  {rhom:6.3f}   {rhom*s2m:10.4f}   {bias2:6.3f}"
          f"   {float(np.mean((Fm.mean(1)-yq)**2)):10.3f}")
print("  rho falls exactly as predicted -- that is the random-forest mechanism, measured.")
print("  But on THIS problem only 4 of 10 features matter, so restricting the choice")
print("  often excludes all of them: sigma^2 and bias^2 both rise, and mtry = d wins.")
print("  m is a hyperparameter with an interior optimum, not a free improvement.")

# ============ 5. bagging attacks variance; boosting attacks bias =============
print()
print("--- 5. BAGGING REDUCES VARIANCE; BOOSTING REDUCES BIAS ---")
stump = grow(Xtr, ytr, 1, 5)
print(f"  one depth-1 stump                      test MSE {np.mean((pred(stump,Xte)-yte)**2):.4f}")
rgs = np.random.default_rng(1); Pb = np.zeros(len(Xte))
for b in range(200):
    idx = rgs.integers(0, 400, 400); Pb += pred(grow(Xtr[idx], ytr[idx], 1, 5), Xte)
print(f"  200 BAGGED stumps                      test MSE {np.mean((Pb/200-yte)**2):.4f}")
for dep in (1, 2):
    F = np.full(400, ytr.mean()); Fv = np.full(len(Xte), ytr.mean())
    for m in range(200):
        T = grow(Xtr, ytr - F, dep, 5)                  # -gradient of squared loss = residual
        F += 0.3*pred(T, Xtr); Fv += 0.3*pred(T, Xte)
    print(f"  200 BOOSTED depth-{dep} trees (nu=0.3)     test MSE {np.mean((Fv-yte)**2):.4f}")
print("  Averaging models that share a bias cannot remove it: bagging stumps barely helps.")
print("  Boosting attacks the bias directly. Note depth-1 boosting is an ADDITIVE model")
print("  and cannot represent the x0*x1 interaction -- depth 2 buys exactly that.")

# ==================== 6. the learning rate / rounds trade-off ===============
print()
print("--- 6. LEARNING RATE vs NUMBER OF ROUNDS (test MSE) ---")
print("   nu      M=20      M=60      M=200")
for nu in (1.0, 0.3, 0.1):
    F = np.full(400, ytr.mean()); Fv = np.full(len(Xte), ytr.mean()); snap = {}
    for m in range(1, 201):
        T = grow(Xtr, ytr - F, 3, 5); F += nu*pred(T, Xtr); Fv += nu*pred(T, Xte)
        if m in (20, 60, 200): snap[m] = np.mean((Fv - yte)**2)
    print(f"  {nu:4.2f}   {snap[20]:8.4f}  {snap[60]:8.4f}  {snap[200]:8.4f}")
print("  nu=1 lands somewhere mediocre in 20 rounds and never improves: each step")
print("  overshoots. Small nu is slower per round and better in the end -- Lecture 3's")
print("  step-size lesson, in function space.")

# ==================== 7. the Newton leaf value, verified ====================
print()
print("--- 7. THE SECOND-ORDER LEAF VALUE  w* = -G/(H + lambda) ---")
res = ytr - ytr.mean()
T   = grow(Xtr, res, 3, 5)
leaf_id = np.zeros(400, int); nxt = [0]
def tag(nd, ii):
    if nd.j < 0:
        leaf_id[ii] = nxt[0]; nxt[0] += 1; return
    m = Xtr[ii, nd.j] <= nd.t; tag(nd.lo, ii[m]); tag(nd.hi, ii[~m])
tag(T, np.arange(400))
K = leaf_id.max() + 1
for lam in (0.0, 5.0, 50.0):
    w_star, shrunk = [], []
    for k in range(K):
        mk = leaf_id == k
        G = float((-res[mk]).sum())      # dL/dF for L = 0.5(y-F)^2 at F=0 is -res
        H = float(mk.sum())              # d2L/dF2 = 1 per point
        w_star.append(-G/(H + lam))
        shrunk.append(res[mk].mean() * mk.sum()/(mk.sum() + lam))
    w_star, shrunk = np.array(w_star), np.array(shrunk)
    print(f"  lambda={lam:5.1f}   max |w* - (shrunken mean residual)| = "
          f"{np.abs(w_star - shrunk).max():.2e}    mean |w*| = {np.abs(w_star).mean():.4f}")
print("  At lambda=0 the Newton value IS the mean residual, so second-order boosting")
print("  reduces to first-order boosting exactly, as the derivation requires.")
print("  Raising lambda shrinks SMALL leaves hardest (H_j is the leaf's point count) --")
print("  which is what you want, since small leaves are the least trustworthy.")
print()
print(f"total run time {time.time() - _t0:.1f} s")

8Exercises

★

A node contains 60 points of class A and 40 of class B. A candidate split sends (40 A, 10 B) left and (20 A, 30 B) right.

Compute the Gini gain, the entropy gain, and the misclassification gain.
A second candidate sends (50 A, 20 B) left and (10 A, 20 B) right. Which split does each criterion prefer?
Verify the identity \(\bar p = \frac{n_L}{n}p_L + \frac{n_R}{n}p_R\) for both splits.

Split 1. Parent \(p=0.6\) (class A). Left \(n_L=50,\ p_L=0.8\); right \(n_R=50,\ p_R=0.4\).

Gini: parent \(2(0.6)(0.4)=0.48\); left \(2(0.8)(0.2)=0.32\); right \(2(0.4)(0.6)=0.48\). Gain \(=0.48-0.5(0.32)-0.5(0.48)=0.48-0.40=\mathbf{0.08}\).

Entropy: parent \(H(0.6)=0.9710\); left \(H(0.8)=0.7219\); right \(H(0.4)=0.9710\). Gain \(=0.9710-0.5(0.7219)-0.5(0.9710)=\mathbf{0.1245}\).

Misclassification: parent \(0.4\); left \(0.2\); right \(0.4\). Gain \(=0.4-0.5(0.2)-0.5(0.4)=\mathbf{0.10}\).

Split 2. Left \(n_L=70,\ p_L=50/70=0.7143\); right \(n_R=30,\ p_R=1/3\).

Gini: left \(2(0.7143)(0.2857)=0.4082\); right \(2(0.3333)(0.6667)=0.4444\). Gain \(=0.48-0.7(0.4082)-0.3(0.4444)=0.48-0.4190=\mathbf{0.0610}\).

Entropy: \(H(0.7143)=0.8631\), \(H(1/3)=0.9183\). Gain \(=0.9710-0.7(0.8631)-0.3(0.9183)=0.9710-0.8797=\mathbf{0.0913}\).

Misclassification: left \(0.2857\), right \(1/3\). Gain \(=0.4-0.7(0.2857)-0.3(0.3333)=0.4-0.3=\mathbf{0.10}\).

2. Gini and entropy both prefer split 1 (0.0800 > 0.0610; 0.1245 > 0.0913). Misclassification scores them identically at 0.10 — both children keep A as the majority in every case, so we are on the linear piece of \(\min(p,1-p)\) and the chord lies on the curve. Third example of the same failure.

3. Split 1: \(0.5(0.8)+0.5(0.4)=0.6\) ✓. Split 2: \(0.7(0.7143)+0.3(0.3333)=0.5+0.1=0.6\) ✓. The identity is just "the counts add".

★★

Prove that a bootstrap sample of size \(n\) drawn with replacement from \(n\) items omits any particular item with probability \((1-1/n)^n\), and find the limit. Then:

What fraction of the training set does one tree actually see, for large \(n\)?
How many of \(B=500\) trees is a given row out-of-bag for, on average?
Explain why OOB error is slightly pessimistic compared with the full forest's test error.
At \(n=10\), how far is \((1-1/n)^n\) from the limit? Does OOB still work at small \(n\)?

Proof. Each of the \(n\) independent draws picks item \(i\) with probability \(1/n\), so misses it with probability \(1-1/n\). The draws are independent, so all \(n\) miss it with probability \((1-1/n)^n\). Taking logs: \(n\log(1-1/n) = n\big(-\tfrac1n-\tfrac{1}{2n^2}-O(n^{-3})\big) = -1-\tfrac{1}{2n}+O(n^{-2}) \to -1\), so the probability \(\to e^{-1}=0.3679\). ∎

1. About \(1-0.368=63.2\%\) of the rows appear at least once (some several times, since the sample still has \(n\) elements).

2. \(500\times0.368\approx 184\) trees. That is a substantial ensemble, which is why the OOB estimate is usable rather than merely suggestive.

3. Because the OOB prediction for a row averages only those ~184 trees, not all 500. By the variance formula, an ensemble of \(B'=184\) has variance \(\rho\sigma^2+\frac{1-\rho}{184}\sigma^2\) versus \(\rho\sigma^2+\frac{1-\rho}{500}\sigma^2\) — strictly larger, so slightly higher expected error. With \(\rho\) of order 0.1 and \(B\) in the hundreds both terms are dominated by \(\rho\sigma^2\) and the gap is negligible; at \(B=20\) it is not.

4. \((1-1/10)^{10}=0.9^{10}=0.3487\) versus \(0.3679\) — about 5% low, matching the \(-\frac{1}{2n}\) correction term. So each row is OOB for ~35% of the trees. The estimator still works; the real problem at small \(n\) is that averaging errors over few rows is noisy, which is a sample-size issue, not a bootstrap issue.

★★

Suppose all \(B\) trees in a forest are identical (\(\rho=1\)).

What does the variance formula give? Verify it directly without the formula.
Now suppose \(\rho<0\) is somehow achievable. What is the smallest possible variance, and what constraint on \(\rho\) makes the formula sensible?
A colleague proposes de-correlating trees by giving each one a completely different random 10% subset of the training rows (no overlap between trees). Using the formula, say what happens to \(\rho\), to \(\sigma^2\), and whether this is a good idea.

1. \(\rho\sigma^2+\frac{1-\rho}{B}\sigma^2 = \sigma^2+0=\sigma^2\). Directly: if every \(f_b\) is the same random variable \(f\), then \(\frac1B\sum_b f_b = f\), whose variance is \(\sigma^2\). Averaging identical things does nothing — which is exactly the failure mode in the quiz.

2. Variance \(\ge 0\) forces \(\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\ge0\), i.e. \(\rho\ge-\frac{1}{B-1}\). At that boundary the variance is exactly zero. This is not an artefact: \(B\) exchangeable variables cannot all be strongly negatively correlated, because their correlation matrix must stay positive semi-definite, and \(-\frac{1}{B-1}\) is precisely where its smallest eigenvalue hits zero. In practice \(\rho>0\) always for bagged trees — they share data.

3. \(\rho\) falls a great deal (disjoint data means the only shared information is the underlying distribution). But each tree now trains on \(n/10\) rows, so its variance \(\sigma^2\) rises sharply and its bias rises too — a tree on 10% of the data is much worse. Total error is \(\text{bias}^2 + \rho\sigma^2 + \frac{1-\rho}{B}\sigma^2\), and the bias term is untouched by averaging. So: usually a bad idea, and it is the same trade-off as \(m\) in a random forest, pushed to an extreme. The code block measures exactly this trade-off for \(m\) and finds the corner solution \(m=d\) on that problem — which is why you tune it rather than trusting a default.

★★★

Derive gradient boosting for logistic loss \(L(y,F)=\log\!\big(1+e^{-yF}\big)\) with \(y\in\{-1,+1\}\) and \(F\) the log-odds.

Compute \(g_i=\partial L/\partial F\) and \(h_i=\partial^2 L/\partial F^2\).
Write the Newton leaf value \(w_j^\star\) and interpret \(H_j\).
What is \(F_0\)?
Show that a point the ensemble already classifies confidently and correctly contributes almost nothing to either \(G_j\) or \(H_j\), and say what that means computationally.
Implement it and check that the training logistic loss decreases monotonically for small \(\nu\).

1. Write \(u=-yF\), so \(L=\log(1+e^{u})\) and \(\frac{\partial L}{\partial u}=\frac{e^u}{1+e^u}=\sigma(u)\). By the chain rule \(\frac{\partial u}{\partial F}=-y\), so

\[ g = -y\,\sigma(-yF) = \frac{-y}{1+e^{yF}}. \]

For the second derivative, using \(\sigma'=\sigma(1-\sigma)\) and another factor of \(-y\), with \(y^2=1\):

\[ h = y^2\,\sigma(-yF)\big(1-\sigma(-yF)\big) = \sigma(yF)\big(1-\sigma(yF)\big) = p(1-p), \quad p=\sigma(F). \]

Note \(h\) does not depend on \(y\) — only on how confident the current model is.

2. \(w_j^\star = -\dfrac{\sum_{i\in R_j} g_i}{\sum_{i\in R_j} h_i + \lambda} = \dfrac{\sum_{i\in R_j} y_i\,\sigma(-y_iF_i)}{\sum_{i\in R_j} p_i(1-p_i)+\lambda}\).

The numerator is a sum of signed "how wrong are we" terms; the denominator is the total curvature. Since \(p(1-p)\le\tfrac14\), \(H_j\) is an effective sample size that counts uncertain points fully and confident points barely at all. That is exactly what XGBoost's min_child_weight constrains: not the number of points in a leaf, but the amount of usable information.

3. \(F_0=\argmin_c\sum_i\log(1+e^{-y_ic})\). Differentiating and setting to zero gives \(\sum_i y_i\sigma(-y_ic)=0\), whose solution is \(F_0=\log\frac{n_+}{n_-}\) — the log-odds of the base rate, the same intercept-only fit as in Lecture 5.

4. Confident and correct means \(y_iF_i \gg 0\). Then \(|g_i|=\sigma(-y_iF_i)\approx e^{-y_iF_i}\to0\) and \(h_i=p(1-p)\to0\), both exponentially fast. So such a point contributes essentially nothing to any leaf's \(G\) or \(H\) and cannot influence the next tree. Computationally, that is the licence to ignore it: it justifies gradient-based one-sided sampling (LightGBM's GOSS), which keeps the large-gradient points and subsamples the rest. It is also the precise sense in which boosting "focuses on hard examples" — a consequence of the derivative, not a design choice.

5. Reuse grow and pred: at each round compute g = -y/(1+np.exp(y*F)) and h = p*(1-p), fit a tree to -g/h weighted by h (or simply to -g for the first-order version), recompute the leaf values as \(-G_j/(H_j+\lambda)\), and update \(F \mathrel{+}= \nu\,h_m\). Track np.log1p(np.exp(-y*F)).mean(): it decreases monotonically for \(\nu\) small. For \(\nu=1\) it need not — the Taylor expansion is only local, exactly as in Lecture 3.

Where this goes next

Trees get nonlinearity by partitioning the input space. The final lecture of Phase B gets it the other way round: keep a single linear boundary, but move to a space where the data is linearly separable — and do it without ever computing the coordinates in that space. That is the kernel trick, and it needs the inner products of Phase A Lecture 8 and the constrained optimisation of Lecture 12.

Previous10 · Naive Bayes and the Generative–Discriminative Divide Next12 · Support Vector Machines: Margins, Duality, and Kernels