Decision Trees & Ensembles

A decision tree recursively partitions the feature space. Each internal node tests one feature against a threshold (for continuous features) or membership (for categorical ones). Each leaf node holds a prediction — a class label or a real value. Training is greedy: at each node, we search over all features and all thresholds for the split that maximizes some impurity reduction.

Formal definition: Given dataset $D = {(x_i, y_i)}_{i=1}^n$, a node splits $D$ into $D_L$ and $D_R$. The quality of the split is: $$\Delta I = I(D) - \frac{|D_L|}{|D|} I(D_L) - \frac{|D_R|}{|D|} I(D_R)$$ where $I$ is an impurity measure. We pick the split that maximizes $\Delta I$.

Analogy: Think of sorting a messy pile of colored balls into boxes. At each step, you draw a line — "balls bigger than 5cm go left, smaller go right" — trying to make each box as pure (single-color) as possible. A decision tree automates this process, trying every possible line in every possible dimension.

The algorithm terminates when a stopping criterion is met: maximum depth, minimum samples per node, or zero impurity. The resulting tree is a piecewise-constant function: it carves the feature space into rectangular regions (axis-aligned hyperrectangles) and assigns one prediction per region.

def build_tree(X, y, depth=0, max_depth=5):
    # Base cases
    if depth == max_depth or len(set(y)) == 1:
        return {'leaf': True, 'value': majority_class(y)}
    
    best_gain, best_feat, best_thresh = 0, None, None
    for feat in range(X.shape[1]):
        for thresh in sorted(set(X[:, feat])):
            left = y[X[:, feat] <= thresh]
            right = y[X[:, feat] > thresh]
            gain = info_gain(y, left, right)  # defined below
            if gain > best_gain:
                best_gain, best_feat, best_thresh = gain, feat, thresh
    
    if best_feat is None:  # no informative split
        return {'leaf': True, 'value': majority_class(y)}
    
    mask = X[:, best_feat] <= best_thresh
    return {
        'leaf': False, 'feat': best_feat, 'thresh': best_thresh,
        'left':  build_tree(X[mask],  y[mask],  depth+1, max_depth),
        'right': build_tree(X[~mask], y[~mask], depth+1, max_depth)
    }

The impurity measure used in the classic ID3 and C4.5 algorithms is Shannon entropy, borrowed from information theory. For a discrete distribution $p = (p_1, \ldots, p_k)$ over $k$ classes: $$H(p) = -\sum_{c=1}^k p_c \log_2 p_c$$ Intuition: entropy measures surprise. A fair coin has maximum entropy (1 bit — you learn one full bit of information when it lands). A two-headed coin has zero entropy — there is no surprise.

Information Gain of a split is the entropy reduction: $$\text{IG}(D, \text{split}) = H(D) - \sum_{\text{child}} \frac{|D_\text{child}|}{|D|} H(D_\text{child})$$ We greedily maximize this at each node.

Gini impurity is an alternative that is faster to compute and produces very similar trees: $$\text{Gini}(p) = 1 - \sum_{c=1}^k p_c^2$$ Both reach 0 for a pure node and are maximized at uniform distributions. sklearn uses Gini by default for classification.

For regression, the impurity is variance (or equivalently, MSE from the node mean): $$I(D) = \frac{1}{|D|}\sum_{i \in D}(y_i - \bar y)^2$$ and we split to minimize weighted child variance — identical logic, different criterion.

import numpy as np

def entropy(y):
    """Shannon entropy of a label array."""
    classes, counts = np.unique(y, return_counts=True)
    p = counts / len(y)
    # convention: 0 * log2(0) = 0
    return -np.sum(p * np.log2(p + 1e-12))

def info_gain(parent, left, right):
    n = len(parent)
    return (entropy(parent)
            - (len(left)/n)  * entropy(left)
            - (len(right)/n) * entropy(right))

# Example:
y = np.array([1,1,1,0,0,0,0,0])  # 3 positives, 5 negatives
print(entropy(y))  # ~0.954 bits — close to maximum

left  = np.array([1,1,1])  # pure positive
right = np.array([0,0,0,0,0])  # pure negative
print(info_gain(y, left, right))  # ~0.954 — perfect split

A fully-grown decision tree is a universal approximator — it can perfectly memorize any training set. Given $n$ training points, a tree of depth $\log_2 n$ can create one leaf per example, achieving zero training error. But this tree has memorized the data, not learned from it. On a new test point, it extrapolates wildly from a single training example.

This is the overfitting problem in its purest form. The tree has high variance: tiny changes to the training set produce dramatically different trees.

How to regularize a decision tree:

1. Max depth — the simplest control. A shallow tree is a low-variance, high-bias model (a stump of depth 1 is just a threshold on one feature).

2. Min samples per leaf — do not split a node if it has fewer than $k$ examples. This prevents splits driven by noise in small subsets.

3. Min impurity decrease — only split if $\Delta I \geq \epsilon$. Ignores marginally useful splits.

4. Post-pruning (reduced-error pruning) — grow the full tree, then walk back up and eliminate subtrees that do not improve validation-set accuracy. This is the C4.5 approach and tends to produce better trees than pre-pruning.

The bias-variance decomposition makes this precise. For any estimator: $$\mathbb{E}[(\hat y - y)^2] = \text{Bias}^2 + \text{Variance} + \text{Noise}$$ A shallow tree: high bias (underfits), low variance. A deep tree: low bias, high variance. There is always a sweet spot — but finding it for a single tree is hard. Ensemble methods sidestep the problem entirely.

from sklearn.tree import DecisionTreeClassifier
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=500, n_features=20,
                            n_informative=5, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

for max_depth in [1, 3, 5, 10, None]:
    clf = DecisionTreeClassifier(max_depth=max_depth, random_state=42)
    clf.fit(X_train, y_train)
    print(f"depth={max_depth:>4}  "
          f"train={clf.score(X_train,y_train):.3f}  "
          f"test={clf.score(X_test,y_test):.3f}")
# depth=   1  train=0.691  test=0.673   <- underfitting
# depth=   3  train=0.854  test=0.813
# depth=   5  train=0.929  test=0.833
# depth=  10  train=0.993  test=0.813
# depth=None  train=1.000  test=0.780   <- overfitting

Leo Breiman's key insight (1996): if a model has high variance, average several independent copies. Averaging $B$ independent estimators each with variance $\sigma^2$ gives combined variance $\sigma^2/B$. But you only have one training set — how do you get many independent models?

Bootstrap sampling: draw $n$ examples from your training set with replacement. Each bootstrap sample is size $n$ but contains ~63.2% unique examples (the rest are duplicates; about 36.8% of examples are left out — these are the out-of-bag (OOB) samples). Train one tree on each bootstrap sample. Predict by majority vote (classification) or averaging (regression).

This is Bagging (Bootstrap Aggregating). Why does it help?

• Each tree sees slightly different data, so the trees make different errors.
• When errors are uncorrelated, averaging cancels them.
• The bias of each tree is preserved (bagging does not reduce bias — this is crucial).
• Variance is reduced by a factor of up to $B$ if the trees were independent (they are not — they share the same training distribution — so in practice variance reduces by less).

OOB error: because each tree was trained without seeing ~37% of the data, you can evaluate it on those OOB examples. The OOB error is a free, unbiased estimate of generalization error — no separate validation set needed. Empirically it closely tracks $k$-fold CV.

import numpy as np
from sklearn.tree import DecisionTreeClassifier

def bagging_predict(X_train, y_train, X_test, B=100, max_depth=None):
    n = len(X_train)
    predictions = np.zeros((B, len(X_test)))
    rng = np.random.default_rng(42)
    
    for b in range(B):
        # Bootstrap sample (with replacement)
        idx = rng.integers(0, n, size=n)
        X_b, y_b = X_train[idx], y_train[idx]
        
        tree = DecisionTreeClassifier(max_depth=max_depth)
        tree.fit(X_b, y_b)
        predictions[b] = tree.predict(X_test)
    
    # Majority vote across B trees
    return np.apply_along_axis(
        lambda col: np.bincount(col.astype(int)).argmax(),
        axis=0, arr=predictions
    )

Bagging helps, but the trees are still correlated: they all use the same features. If one feature is very dominant (say, income in a loan dataset), nearly every tree will split on it near the root, producing trees that look alike. Averaging correlated trees reduces variance less than averaging uncorrelated ones.

Random Forests (Breiman 2001) fix this by adding feature randomization: at each split, instead of considering all $p$ features, consider a random subset of size $m$ (typically $m \approx \sqrt{p}$ for classification, $m \approx p/3$ for regression). This forces different trees to rely on different features, decorrelating them.

Why this works mathematically: if each tree has variance $\sigma^2$ and pairwise correlation $\rho$, the ensemble variance is: $$\text{Var}(\text{RF}) = \rho \sigma^2 + \frac{1-\rho}{B}\sigma^2$$ As $B \to \infty$, the second term vanishes, but the first term $\rho\sigma^2$ is an irreducible floor. Reducing $\rho$ (via feature subsampling) reduces this floor.

Feature importance: a natural byproduct of Random Forests is a measure of each feature's importance. The Mean Decrease in Impurity (MDI) sums, across all trees and all nodes, the impurity decrease attributable to each feature, weighted by the number of samples reaching that node. Features that are used at high-level splits with large impurity reductions score high. This gives a fast, heuristic feature ranking (though it can be biased toward high-cardinality features — permutation importance is a more reliable alternative).

from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_breast_cancer
import numpy as np

data = load_breast_cancer()
X, y = data.data, data.target

rf = RandomForestClassifier(
    n_estimators=200,   # B trees
    max_features='sqrt',  # m = sqrt(p) at each split
    max_depth=None,       # grow full trees
    oob_score=True,       # compute OOB error for free
    n_jobs=-1,
    random_state=42
)
rf.fit(X, y)
print(f"OOB accuracy: {rf.oob_score_:.4f}")

# Feature importances (MDI)
importances = rf.feature_importances_
top5 = np.argsort(importances)[::-1][:5]
for i in top5:
    print(f"  {data.feature_names[i]:30s}  {importances[i]:.4f}")

Bagging trains trees in parallel, each on a random subsample. Boosting takes the opposite approach: train trees sequentially, with each tree focusing on the examples the previous trees got wrong.

AdaBoost (Freund & Schapire, 1997) is the classic formulation. Maintain a weight $w_i$ for each training example, initially $w_i = 1/n$. At each round $t$:

1. Train a weak learner $h_t$ (typically a stump) on the weighted dataset.
2. Compute its weighted error $\epsilon_t = \sum_i w_i \cdot \mathbf{1}[h_t(x_i) \neq y_i]$.
3. Assign it a vote weight $\alpha_t = \frac{1}{2}\ln\frac{1-\epsilon_t}{\epsilon_t}$ (larger for more accurate classifiers).
4. Increase weights on misclassified examples, decrease on correct ones.

Final prediction: $H(x) = \text{sign}\left(\sum_t \alpha_t h_t(x)\right)$.

The remarkable theorem (Schapire): as long as each weak learner does better than random ($\epsilon_t < 0.5$), the training error of AdaBoost converges to zero exponentially fast. More strikingly, even after training error is zero, test error often continues to improve as more trees are added — boosting increases the classification margin even after zero training error, which is why it often avoids overfitting despite zero training error.

Bias-variance view: Unlike bagging (which reduces variance), boosting primarily reduces bias — it takes a high-bias model (stumps) and iteratively corrects its systematic errors, building up a low-bias ensemble.

# AdaBoost core loop (simplified, binary ±1 labels)
import numpy as np

def adaboost(X, y, T=50):
    n = len(y)
    w = np.ones(n) / n          # uniform weights
    alphas, stumps = [], []
    
    for t in range(T):
        stump = fit_stump(X, y, w)  # weighted decision stump
        pred = stump.predict(X)     # predictions in {-1, +1}
        
        err = np.sum(w * (pred != y))  # weighted error
        alpha = 0.5 * np.log((1 - err) / (err + 1e-10))
        
        # Update weights: wrong examples get higher weight
        w = w * np.exp(-alpha * y * pred)
        w /= w.sum()               # renormalize
        
        alphas.append(alpha)
        stumps.append(stump)
    
    def predict(X_new):
        votes = sum(a * s.predict(X_new) for a, s in zip(alphas, stumps))
        return np.sign(votes)
    
    return predict

Gradient Boosting (Friedman 2001) reframes boosting as gradient descent in the space of functions. This unifying view explains AdaBoost as a special case and allows you to use any differentiable loss function.

The key idea: suppose you want to minimize a loss $L(y, F(x))$ over functions $F$. Gradient descent would move in the direction of the negative gradient. But $F$ is not a parameter vector — it is a function over the input space. The "gradient" at each training point is the negative partial derivative: $$r_i = -\left[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}\right]{F=F{t-1}}$$ These $r_i$ are called pseudo-residuals. They tell us, for each example, in which direction and by how much to adjust our current prediction.

Algorithm:

1. Initialize $F_0(x) = \arg\min_\gamma \sum_i L(y_i, \gamma)$ (e.g., the mean for MSE).
2. For $t = 1, \ldots, T$: a. Compute pseudo-residuals $r_i$ from the current model $F_{t-1}$. b. Fit a tree $h_t$ to the pseudo-residuals (regression tree). c. Compute the optimal step size $\gamma_t$ and update: $F_t = F_{t-1} + \nu \gamma_t h_t$, where $\nu$ is the learning rate.

For MSE loss ($L = (y-F)^2/2$), the pseudo-residuals are just the ordinary residuals $r_i = y_i - F_{t-1}(x_i)$. Each tree fits the residuals of the previous ensemble — this is the clearest explanation of why boosting "corrects mistakes".

For log-loss (classification), pseudo-residuals involve probabilities and link functions. This recovers LogitBoost.

Learning rate / shrinkage: scaling each tree's contribution by $\nu < 1$ slows convergence but dramatically improves generalization. More trees + smaller $\nu$ almost always beats fewer trees + larger $\nu$.

import numpy as np
from sklearn.tree import DecisionTreeRegressor

def gradient_boosting_mse(X_train, y_train, X_test,
                          n_trees=100, max_depth=3, lr=0.1):
    # Initialize: predict the mean
    F = np.full(len(X_train), y_train.mean())
    F_test = np.full(len(X_test), y_train.mean())
    
    for _ in range(n_trees):
        # Pseudo-residuals for MSE: r_i = y_i - F(x_i)
        residuals = y_train - F
        
        # Fit a tree to residuals
        tree = DecisionTreeRegressor(max_depth=max_depth)
        tree.fit(X_train, residuals)
        
        # Update predictions (learning rate = shrinkage)
        F      += lr * tree.predict(X_train)
        F_test += lr * tree.predict(X_test)
    
    return F_test

# sklearn's GradientBoostingRegressor does the same thing,
# but with efficient histogram binning and regularization.

XGBoost (Chen & Guestrin 2016) takes gradient boosting and makes it rigorous, fast, and regularized. It is the algorithm that won the most Kaggle competitions between 2014 and 2019.

The objective function at each step $t$ is a second-order Taylor approximation of the loss plus a regularization term: $$\mathcal{L}^{(t)} \approx \sum_{i=1}^n \left[g_i f_t(x_i) + \frac{1}{2} h_i f_t(x_i)^2\right] + \Omega(f_t)$$ where $g_i = \partial_{\hat y} L(y_i, \hat y_{t-1})$ and $h_i = \partial^2_{\hat y} L(y_i, \hat y_{t-1})$ are the first and second derivatives of the loss at the current prediction. The regularization term is: $$\Omega(f) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2$$ where $T$ is the number of leaves, $w_j$ are leaf weights, $\gamma$ penalizes tree complexity, and $\lambda$ is L2 regularization on leaf weights.

Why second-order? Using the Hessian $h_i$ gives a much better approximation of the true loss surface, enabling larger, more aggressive steps without instability — similar to how Newton's method converges faster than gradient descent.

Key engineering innovations in XGBoost:

• Column subsampling (like Random Forests) at the tree or node level.
• Approximate greedy algorithm: instead of exact split finding (expensive), build a histogram of feature values and consider only quantile-based split candidates.
• Sparsity-aware split finding: handles missing values natively — at each split, both directions are tried and the better default direction is learned.
• Cache-aware access: column blocks are pre-sorted and stored for efficient in-cache access.
• Out-of-core computation: block compression and sharding for datasets that exceed RAM.

These innovations give XGBoost a 10x–100x speedup over naive gradient boosting, enabling it to run on millions of examples in minutes.

import xgboost as xgb
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import numpy as np

X, y = make_regression(n_samples=10000, n_features=50,
                       n_informative=20, noise=10, random_state=42)
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)

model = xgb.XGBRegressor(
    n_estimators=500,
    learning_rate=0.05,    # shrinkage (nu)
    max_depth=6,
    subsample=0.8,         # row subsampling (bagging-like)
    colsample_bytree=0.8,  # feature subsampling per tree
    reg_lambda=1.0,        # L2 on leaf weights
    gamma=0.1,             # minimum gain to split a node
    early_stopping_rounds=20,
    eval_metric='rmse',
    random_state=42
)

model.fit(X_train, y_train,
          eval_set=[(X_val, y_val)],
          verbose=False)

pred = model.predict(X_val)
print(f"RMSE: {np.sqrt(mean_squared_error(y_val, pred)):.3f}")
print(f"Best iteration: {model.best_iteration}")

The choice between bagging/random forests and boosting is fundamentally about the bias-variance tradeoff.

Bagging/Random Forests reduce variance. They work best when you have a high-variance base model (deep trees). If the base model already has low variance, averaging copies of it helps little. Key insight: increasing the number of trees in a random forest never increases bias — the generalization error converges monotonically. Random forests are also robust: outliers in the training set affect only the trees that happen to include them in their bootstrap sample.

Boosting primarily reduces bias. It takes a simple, biased model (stumps) and builds up complexity iteratively. Boosting is more sensitive to noisy labels and outliers, because it deliberately upweights hard examples — which may include mislabeled points. It requires careful tuning of the learning rate and number of trees.

Practical heuristic (from empirical literature):

• If your dataset is large and noisy, prefer Random Forest (more robust, less tuning).
• If your dataset is clean and you want maximum accuracy on tabular data, prefer XGBoost or LightGBM.
• Random Forests are embarrassingly parallelizable; boosting is inherently sequential (though data parallelism is possible within each tree).
• LightGBM uses a leaf-wise (best-first) tree growth strategy instead of the level-wise (breadth-first) strategy used by XGBoost/sklearn, which often produces better results on large datasets.

Modern context (2020s): For purely tabular data, tree ensembles (XGBoost, LightGBM, CatBoost) still routinely outperform deep learning. For images, audio, text: neural networks dominate. For tabular data with high-cardinality categorical features: CatBoost's native category handling is often superior.

# Visualizing bias-variance tradeoff empirically
import numpy as np
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.tree import DecisionTreeRegressor

def bias_variance_experiment(model_factory, X_train_sets, X_test, y_test):
    """
    Train model on many bootstrap training sets.
    Compute bias^2, variance, and noise empirically.
    """
    preds = np.array([
        model_factory().fit(X_tr, y_tr).predict(X_test)
        for X_tr, y_tr in X_train_sets
    ])
    # E[f(x)] across training sets
    mean_pred = preds.mean(axis=0)
    bias2  = np.mean((mean_pred - y_test)**2)
    variance = np.mean(preds.var(axis=0))
    print(f"Bias^2={bias2:.3f}  Variance={variance:.3f}  "
          f"Sum={bias2+variance:.3f}")