Ensemble Theory: Bagging, Random Forests & Boosting

Ensemble learning combines multiple base models (weak learners) to produce a unified predictor whose predictive performance and robustness exceed that of any individual constituent. The two foundational paradigms in ensemble theory are Bagging (Variance Reduction) and Boosting (Bias Reduction).

Variance AsymptoteStatistical Physics

rho * sigma^2

Variance of an ensemble of B estimators converges to pairwise correlation rho * sigma^2 as B approaches infinity.

1. Mathematical Proof: Variance Reduction in Bagging#

Consider an ensemble of $B$ base estimators $f_1(x), f_2(x), \dots, f_B(x)$ , each with variance $\sigma^2$ and pairwise correlation coefficient $\rho = \text{Corr}(f_i, f_j)$ for $i \ne j$ .

The ensemble prediction is the simple arithmetic mean:

\bar{f}(x) = \frac{1}{B} \sum_{b=1}^B f_b(x)

Architecture Diagram

Rendering diagram...

Step 1: Expand Ensemble Variance

By the definition of the variance of a linear combination of random variables:

\text{Var}(\bar{f}) = \text{Var}\left( \frac{1}{B} \sum_{b=1}^B f_b \right) = \frac{1}{B^2} \left[ \sum_{b=1}^B \text{Var}(f_b) + \sum_{i=1}^B \sum_{j \ne i}^B \text{Cov}(f_i, f_j) \right]

Step 2: Substitute Variance and Covariance

Since each estimator has variance $\sigma^2$ and covariance $\text{Cov}(f_i, f_j) = \rho \sigma^2$ :

\text{Var}(\bar{f}) = \frac{1}{B^2} \left[ B \sigma^2 + B(B - 1) \rho \sigma^2 \right]

Step 3: Simplify the Expression

Divide both terms by $B^2$ :

\text{Var}(\bar{f}) = \frac{\sigma^2}{B} + \frac{B(B - 1) \rho \sigma^2}{B^2} = \frac{\sigma^2}{B} + \left( 1 - \frac{1}{B} \right) \rho \sigma^2

Rearranging the terms:

\text{Var}(\bar{f}) = \rho \sigma^2 + \frac{1 - \rho}{B} \sigma^2

Asymptotic Implications & Random Forests#

Taking the limit as the number of estimators $B \to \infty$ :

\lim_{B \to \infty} \text{Var}(\bar{f}) = \rho \sigma^2

Why Random Forests Outperform Simple Bagging

Simple bagging of decision trees yields high pairwise correlation

\rho

because strong features dominate the top splits of every tree. Random Forests eliminate this bottleneck by introducing Feature Subspace Sampling: at each split, only a random subset of

m \approx \sqrt{d}

features is evaluated. This forces trees to explore alternative predictive pathways, driving pairwise correlation

\rho \to 0

and drastically reducing overall ensemble variance

\rho \sigma^2

2. Derivation of Out-of-Bag (OOB) Error#

In bootstrap sampling, $N$ instances are drawn uniformly at random with replacement from a training set of size $N$ .

The probability that a specific training sample $x_i$ is not selected in a single draw is $1 - \frac{1}{N}$ .

The probability that $x_i$ is completely excluded across all $N$ independent draws is:

P(x_i \notin \text{Bootstrap}) = \left( 1 - \frac{1}{N} \right)^{\!N}

Taking the calculus limit as dataset size $N \to \infty$ :

\lim_{N \to \infty} \left( 1 - \frac{1}{N} \right)^{\!N} = \frac{1}{e} \approx 0.367879 \dots \approx 36.8\%

Thus, approximately $36.8\%$ of training samples are withheld from each individual tree. These Out-Of-Bag (OOB) samples serve as built-in cross-validation test sets, evaluating generalization error without dedicated validation splits.

3. First-Principles Derivation of AdaBoost#

AdaBoost sequentially trains base learners $h_t(x) \in \{-1, +1\}$ by reweighting training samples to focus on instances misclassified by earlier rounds.

Architecture Diagram

Rendering diagram...

Deriving the Estimator Weight $\alpha_t$ #

AdaBoost minimizes the exponential loss over additive ensemble $F_t(x) = F_{t-1}(x) + \alpha_t h_t(x)$ :

E = \sum_{i=1}^N \exp\Big( -y_i \big( F_{t-1}(x_i) + \alpha_t h_t(x_i) \big) \Big) = \sum_{i=1}^N w_i^{(t)} \exp\Big( -\alpha_t y_i h_t(x_i) \Big)

Where sample weights $w_i^{(t)} = \exp(-y_i F_{t-1}(x_i))$ .

Step 1: Split Loss into Correct and Incorrect Predictions

Because $y_i, h_t(x_i) \in \{-1, +1\}$ , the product $y_i h_t(x_i)$ is $+1$ for correct predictions and $-1$ for incorrect predictions:

E = \sum_{y_i = h_t(x_i)} w_i^{(t)} e^{-\alpha_t} + \sum_{y_i \ne h_t(x_i)} w_i^{(t)} e^{\alpha_t}

Step 2: Express in Terms of Weighted Error $\epsilon_t$

Let $W = \sum_{i=1}^N w_i^{(t)}$ denote total sample weight. The weighted training error is:

\epsilon_t = \frac{\sum_{y_i \ne h_t(x_i)} w_i^{(t)}}{W} \implies \sum_{y_i \ne h_t(x_i)} w_i^{(t)} = W \epsilon_t

The weight of correctly classified samples is:

\sum_{y_i = h_t(x_i)} w_i^{(t)} = W (1 - \epsilon_t)

Substituting these into the loss expression:

E = W (1 - \epsilon_t) e^{-\alpha_t} + W \epsilon_t e^{\alpha_t}

Step 3: Differentiate with Respect to $\alpha_t$ and Solve for Zero

To minimize exponential loss with respect to $\alpha_t$ :

\frac{dE}{d\alpha_t} = -W (1 - \epsilon_t) e^{-\alpha_t} + W \epsilon_t e^{\alpha_t} = 0

Divide by $W$ and multiply by $e^{\alpha_t}$ :

-(1 - \epsilon_t) + \epsilon_t e^{2\alpha_t} = 0 \implies \epsilon_t e^{2\alpha_t} = 1 - \epsilon_t

e^{2\alpha_t} = \frac{1 - \epsilon_t}{\epsilon_t}

Taking the natural logarithm and dividing by 2:

\alpha_t = \frac{1}{2} \ln\left( \frac{1 - \epsilon_t}{\epsilon_t} \right)

When $\epsilon_t < 0.5$ , $\alpha_t > 0$ . As error $\epsilon_t \to 0$ , estimator weight $\alpha_t \to \infty$ .

4. XGBoost: Second-Order Taylor Expansion#

Modern Gradient Boosted Machines like XGBoost optimize arbitrary loss functions using second-order Taylor approximations.

At boosting iteration $t$ , the regularized objective is:

\mathcal{L}^{(t)} = \sum_{i=1}^N l\Big(y_i, \hat{y}_i^{(t-1)} + f_t(x_i)\Big) + \Omega(f_t)

Where tree regularization $\Omega(f_t) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2$ .

Expanding $l(y_i, \hat{y}_i^{(t-1)} + f_t(x_i))$ via a second-order Taylor series around $\hat{y}_i^{(t-1)}$ :

\mathcal{L}^{(t)} \approx \sum_{i=1}^N \left[ l(y_i, \hat{y}_i^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) \right] + \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2

Where first and second order gradients are defined as:

g_i = \left[ \frac{\partial l(y_i, \hat{y})}{\partial \hat{y}} \right]_{\hat{y} = \hat{y}_i^{(t-1)}}, \qquad h_i = \left[ \frac{\partial^2 l(y_i, \hat{y})}{\partial \hat{y}^2} \right]_{\hat{y} = \hat{y}_i^{(t-1)}}

Grouping samples mapped to leaf node $j$ ( $I_j = \{i : q(x_i) = j\}$ ):

\mathcal{L}^{(t)} = \sum_{j=1}^T \left[ \left( \sum_{i \in I_j} g_i \right) w_j + \frac{1}{2} \left( \sum_{i \in I_j} h_i + \lambda \right) w_j^2 \right] + \gamma T

Let $G_j = \sum_{i \in I_j} g_i$ and $H_j = \sum_{i \in I_j} h_i$ . Setting the derivative with respect to leaf weight $w_j$ to zero:

G_j + (H_j + \lambda) w_j = 0 \implies w_j^* = -\frac{G_j}{H_j + \lambda}

Substituting optimal weights $w_j^*$ back gives the optimal split quality metric:

\tilde{\mathcal{L}}^{(t)} = -\frac{1}{2} \sum_{j=1}^T \frac{G_j^2}{H_j + \lambda} + \gamma T

5. From-Scratch Implementation: Zero Dependencies#

The following implementations construct an AdaBoost Classifier with Decision Stumps from scratch in Python and Rust using only built-in language primitives.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import math
from typing import List, Tuple

class DecisionStump:
    def __init__(self):
        self.polarity: int = 1
        self.feature_idx: int = 0
        self.threshold: float = 0.0
        self.alpha: float = 0.0

    def predict(self, X: List[List[float]]) -> List[int]:
        num_samples = len(X)
        predictions = [1] * num_samples
        for i in range(num_samples):
            val = X[i][self.feature_idx]
            if (self.polarity == 1 and val < self.threshold) or (self.polarity == -1 and val > self.threshold):
                predictions[i] = -1
        return predictions

class AdaBoostFromScratch:
    def __init__(self, n_estimators: int = 10):
        self.n_estimators = n_estimators
        self.estimators: List[DecisionStump] = []

    def fit(self, X: List[List[float]], y: List[int]) -> None:
        num_samples = len(X)
        num_features = len(X[0])
        # Initialize uniform sample weights
        w = [1.0 / num_samples] * num_samples

        for _ in range(self.n_estimators):
            stump = DecisionStump()
            min_error = float("inf")

            # Search across all features and unique thresholds
            for feat in range(num_features):
                unique_vals = set(row[feat] for row in X)
                for thresh in unique_vals:
                    for polarity in [1, -1]:
                        # Generate predictions for this candidate stump
                        error = 0.0
                        for i in range(num_samples):
                            val = X[i][feat]
                            pred = -1 if ((polarity == 1 and val < thresh) or (polarity == -1 and val > thresh)) else 1
                            if pred != y[i]:
                                error += w[i]

                        if error < min_error:
                            min_error = error
                            stump.polarity = polarity
                            stump.threshold = thresh
                            stump.feature_idx = feat

            # Clip error to avoid division by zero or negative logarithms
            eps = max(1e-10, min(min_error, 1.0 - 1e-10))
            stump.alpha = 0.5 * math.log((1.0 - eps) / eps)

            # Update sample weights: w_i * exp(-alpha * y_i * h(x_i))
            preds = stump.predict(X)
            for i in range(num_samples):
                w[i] *= math.exp(-stump.alpha * y[i] * preds[i])

            # Normalize sample weights
            w_sum = sum(w)
            for i in range(num_samples):
                w[i] /= w_sum

            self.estimators.append(stump)

    def predict(self, X: List[List[float]]) -> List[int]:
        # Form additive sum: sum alpha_t * h_t(x)
        num_samples = len(X)
        final_scores = [0.0] * num_samples
        for stump in self.estimators:
            preds = stump.predict(X)
            for i in range(num_samples):
                final_scores[i] += stump.alpha * preds[i]
        return [1 if s >= 0.0 else -1 for s in final_scores]

if __name__ == "__main__":
    X_train = [[1.0, 2.0], [1.5, 1.0], [2.0, 2.5], [6.0, 7.0], [7.0, 8.0], [8.0, 6.5]]
    y_train = [-1, -1, -1, 1, 1, 1]

    clf = AdaBoostFromScratch(n_estimators=5)
    clf.fit(X_train, y_train)
    preds = clf.predict([[1.2, 1.8], [7.5, 7.0]])
    print(f"Predictions: {preds}")

6. Active Knowledge Checks#

Concept check

In Random Forests, why does feature subspace sampling (evaluating only sqrt(d) random features per node) reduce ensemble test error compared to standard bagging?