← All Posts
From Scratch · C++ · Machine Learning · Random Forest

A Random Forest in C++

A single deep tree overfits: it can drive training accuracy to $100\%$ while generalizing poorly. A random forest fixes that by averaging many decorrelated trees. Two ingredients on top of Part 1's tree do it: each tree trains on a bootstrap sample, and each split considers only a random subset of features. The payoff in this part: an overfit single tree at $70\%$ test accuracy becomes a forest at $79\%$, with a free $85\%$ out-of-bag estimate.

Step 1: Feature Sub-Sampling in the Tree

The only change to the tree is that, at every node, it shuffles the features and considers just the first $\texttt{maxFeatures}$ of them. The standard choice for classification is

$$\texttt{mtry} = \big\lceil \sqrt{D}\, \big\rceil.$$

Each tree also carries its own random generator so the forest is reproducible.

#include <bits/stdc++.h>
using namespace std;
using Matrix = vector<vector<double>>;
struct Node { bool leaf=false; int feature=-1; double threshold=0; int prediction=-1; Node *L=nullptr,*R=nullptr; };

struct Tree {
    int maxDepth, minSamples, K, maxFeatures;
    mt19937 rng;
    Node* root = nullptr;
    Tree(int maxDepth, int minSamples, int K, int maxFeatures, unsigned seed)
        : maxDepth(maxDepth), minSamples(minSamples), K(K),
          maxFeatures(maxFeatures), rng(seed) {}

    static double gini(const vector<int>& counts, int n) {
        if (n == 0) return 0.0;
        double impurity = 1.0;
        for (int c : counts) {
            double p = (double)c / n;
            impurity -= p * p;
        }
        return impurity;
    }

    Node* build(vector<int> idx, const Matrix& X, const vector<int>& y, int depth) {
        Node* nd = new Node();
        int n = idx.size();
        vector<int> cnt(K, 0);
        for (int i : idx) cnt[y[i]]++;
        double parentGini = gini(cnt, n);

        if (depth >= maxDepth || n < minSamples || parentGini == 0.0) {
            nd->leaf = true;
            nd->prediction = max_element(cnt.begin(), cnt.end()) - cnt.begin();
            return nd;
        }

        // consider only a random subset of mtry features at this node
        int D = X[0].size();
        vector<int> feats(D);
        iota(feats.begin(), feats.end(), 0);
        shuffle(feats.begin(), feats.end(), rng);
        feats.resize(min(maxFeatures, D));

        double bestGain = 0.0;
        int bestFeat = -1;
        double bestThr = 0.0;
        for (int f : feats) {
            vector<int> ord = idx;
            sort(ord.begin(), ord.end(),
                 [&](int a, int b) { return X[a][f] < X[b][f]; });

            vector<int> left(K, 0), right = cnt;
            for (int t = 0; t + 1 < n; ++t) {
                int i = ord[t];
                left[y[i]]++;
                right[y[i]]--;
                if (X[ord[t]][f] == X[ord[t + 1]][f]) continue;

                int nl = t + 1, nr = n - nl;
                double weighted = (double)nl / n * gini(left,  nl)
                                + (double)nr / n * gini(right, nr);
                double gain = parentGini - weighted;
                if (gain > bestGain) {
                    bestGain = gain;
                    bestFeat = f;
                    bestThr = (X[ord[t]][f] + X[ord[t + 1]][f]) / 2;
                }
            }
        }

        if (bestFeat < 0) {
            nd->leaf = true;
            nd->prediction = max_element(cnt.begin(), cnt.end()) - cnt.begin();
            return nd;
        }

        nd->feature = bestFeat;
        nd->threshold = bestThr;
        vector<int> li, ri;
        for (int i : idx)
            (X[i][bestFeat] < bestThr ? li : ri).push_back(i);
        nd->L = build(li, X, y, depth + 1);
        nd->R = build(ri, X, y, depth + 1);
        return nd;
    }

    void fit(vector<int> idx, const Matrix& X, const vector<int>& y) {
        root = build(idx, X, y, 0);
    }

    int predict(const vector<double>& x) const {
        Node* n = root;
        while (!n->leaf)
            n = (x[n->feature] < n->threshold) ? n->L : n->R;
        return n->prediction;
    }
};

Note that fit now takes an explicit index list, so the forest can hand each tree its own bootstrap sample.

Step 2: Bootstrap, Vote, and Out-of-Bag

The forest trains nTrees trees, each on $N$ points drawn with replacement from the training set. It records which samples landed in each tree's bag, predicts by majority vote, and uses the left-out samples to estimate test error.

struct Forest {
    int nTrees, K;
    vector<Tree> trees;
    vector<vector<char>> inbag;       // inbag[t][i] = 1 if sample i trained tree t

    Forest(int nt, int k) : nTrees(nt), K(k) {}

    void fit(const Matrix& X, const vector<int>& y) {
        int N = X.size(), D = X[0].size();
        int mtry = max(1, (int)round(sqrt((double)D)));
        mt19937 rng(7);
        uniform_int_distribution<int> pick(0, N - 1);
        for (int t = 0; t < nTrees; ++t) {
            vector<int> idx(N); vector<char> bag(N, 0);
            for (int i = 0; i < N; ++i) { int j = pick(rng); idx[i] = j; bag[j] = 1; }  // bootstrap
            Tree tr(/*depth*/20, /*minSamples*/2, K, mtry, 1000 + t);
            tr.fit(idx, X, y);
            trees.push_back(move(tr));
            inbag.push_back(bag);
        }
    }
    int predict(const vector<double>& x) const {
        vector<int> v(K, 0);
        for (auto& t : trees) v[t.predict(x)]++;            // majority vote
        return max_element(v.begin(), v.end()) - v.begin();
    }
    double oobAccuracy(const Matrix& X, const vector<int>& y) {
        int ok = 0, counted = 0;
        for (size_t i = 0; i < X.size(); ++i) {
            vector<int> v(K, 0); bool any = false;
            for (int t = 0; t < nTrees; ++t)
                if (!inbag[t][i]) { v[trees[t].predict(X[i])]++; any = true; }  // only OOB trees
            if (any) {
                counted++;
                if ((int)(max_element(v.begin(), v.end()) - v.begin()) == y[i]) ok++;
            }
        }
        return 100.0 * ok / counted;
    }
};
Out-of-bag is the free validation set. Because each bootstrap leaves out roughly $37\%$ of the data ($e^{-1}$), every training point is out-of-bag for about a third of the trees. Predicting it with only those trees gives an honest test estimate without ever holding data back.

Running It

To make the difference visible, use ten features where only the first four carry signal and the rest are noise, plus $8\%$ label flips. A single deep tree latches onto the noise features; the forest does not.

int main() {
    mt19937 rng(11);
    normal_distribution<double> g(0, 1);
    uniform_real_distribution<double> u(0, 1);
    int N = 600, D = 10, K = 2;
    Matrix X; vector<int> y;
    for (int i = 0; i < N; ++i) {
        vector<double> r(D);
        for (int f = 0; f < D; ++f) r[f] = g(rng);          // 10 features, 6 are noise
        int lab = (0.9*r[0] + 0.8*r[1] - 1.0*r[2] + 0.5*r[3] > 0) ? 1 : 0;
        if (u(rng) < 0.08) lab ^= 1;                        // label noise
        X.push_back(r); y.push_back(lab);
    }
    /* ... shuffle, 75/25 split into Xtr/ytr, Xte/yte; trainRows = 0..|Xtr|-1 ... */

    Tree single(20, 2, K, /*mtry*/D, 5);                    // all features = a plain tree
    single.fit(trainRows, Xtr, ytr);
    Forest forest(/*nTrees*/120, K);
    forest.fit(Xtr, ytr);
    // ... report accuracies ...
}

Compiled and run, the contrast is the whole point of bagging:

single deep tree: train=100.0% test=70.0% random forest: test=78.7% OOB=84.7% (mtry=3 of 10)

The single tree memorizes the training set ($100\%$) but only reaches $70\%$ on new data. The forest, by averaging $120$ decorrelated trees, lifts test accuracy to $78.7\%$, and its out-of-bag estimate of $84.7\%$ comes for free. Summing each feature's total impurity drop across the forest would, as a bonus, rank the informative features above the noise ones.