4.6.2 Random Forests

A random forest is a bagged ensemble of decision trees with an extra source of randomness: each split sees only a random subset of the features. This makes the trees' errors less similar, so averaging them becomes more useful.

In practice, random forests are among the most dependable models for tabular data. They need little preparation, work with features on any scale, capture interactions between features, and perform reasonably well with their default settings.

That makes a random forest a good first strong model after a simple baseline, especially when contest time is short. Later, it is also a useful point of comparison for gradient boosting, which can be more accurate but needs more care.

A fully grown decision tree can be unstable: small changes in its training data can change its splits and predictions. Random forests reduce this instability while retaining the tree's ability to model non-linear relationships and feature interactions.

Bootstrap samples and bagging

Bagging trains each tree on a different bootstrap sample. This gives the trees different views of the data. Rows not chosen for a particular tree are out of bag for that tree and become useful later for evaluation.

Before adding the feature randomness that defines a random forest, the following comparison shows what bootstrap sampling and averaging alone do. It compares one unrestricted tree with 100 bagged trees on the breast cancer training split from Decision Trees, using 5-fold cross-validation.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import BaggingClassifier
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.tree import DecisionTreeClassifier

cancer = load_breast_cancer()
X_train, X_val, y_train, y_val = train_test_split(
    cancer.data, cancer.target, test_size=0.3,
    stratify=cancer.target, random_state=0,
)
tree = DecisionTreeClassifier(random_state=0)
bagging = BaggingClassifier(DecisionTreeClassifier(),
                            n_estimators=100, random_state=0)
for name, model in [("tree", tree), ("bagging", bagging)]:
    scores = cross_val_score(model, X_train, y_train, cv=5)
    print(name, scores.round(3), round(scores.mean(), 3))

tree [0.95  0.888 0.912 0.899 0.975] 0.925
bagging [0.962 0.95  0.962 0.924 0.975] 0.955

Bagging raises the mean accuracy from 0.925 to 0.955, and the worst fold from 0.888 to 0.924.

Random feature subsets

Bagged trees remain strongly correlated, because the strongest features win the upper splits in most bootstrap samples. When the bagging model above is fitted to the whole training split, worst perimeter is the root feature of 50 of its 100 trees, and only five different features appear at a root.

The number of features chosen per tree is the hyperparameter max_features, introduced in Decision Trees, and $m = d$, where d is the number of features, gives bagging. By default the scikit-learn classifier uses $m = \sqrt{d}$ rounded down, 5 of the 30 breast cancer features. Its 100 trees then have 14 different root features, and worst perimeter is the root of only 19.

function random_forest(X, y, n_trees, m):
    n = number of samples in X
    trees = []
    repeat n_trees times:
        rows = n random row indices, drawn with replacement
        tree = full CART tree grown on X[rows], y[rows],
               trying only m random features at each split
        trees.append(tree)
    return trees

function forest_predict(trees, x):
    if regression:
        return mean of tree.predict(x) over all trees
    else:  # classification
        return the class with the highest mean probability

Measuring the correlation

The error-correlation formula from the ensemble introduction can be tested on the validation split. For a sample with label $y \in \{0, 1\}$, tree $b$ predicts a probability $p_b$ of class 1 (benign) and has error $e_b = p_b - y$. The code estimates $\sigma^2$ and $\rho$ from these errors for three values of max_features and compares the formula with the forest's actual squared error.

import numpy as np
from sklearn.ensemble import RandomForestClassifier

B = 100
for max_features in [None, "sqrt", 1]:
    forest = RandomForestClassifier(n_estimators=B,
                                    max_features=max_features,
                                    random_state=0)
    forest.fit(X_train, y_train)
    # Row b holds the errors of tree b on the validation samples
    E = np.array([t.predict_proba(X_val)[:, 1] - y_val
                  for t in forest.estimators_])
    sigma2 = (E ** 2).mean()
    # Mean correlation over all pairs of trees (upper triangle)
    rho = np.corrcoef(E)[np.triu_indices(B, k=1)].mean()
    formula = rho * sigma2 + (1 - rho) * sigma2 / B
    p = forest.predict_proba(X_val)[:, 1]
    print(max_features, round(sigma2, 4), round(rho, 3),
          round(formula, 4), round(((p - y_val) ** 2).mean(), 4),
          round(forest.score(X_val, y_val), 3))

The formula nearly equals the forest's squared error in every row. From bagging to "sqrt", the correlation falls from 0.527 to 0.460 while each tree becomes only slightly worse, and the forest's accuracy rises from 0.924 to 0.953. With one feature per split, the correlation falls further, but the trees become so much worse that the forest loses accuracy. max_features thus trades the accuracy of single trees against their correlation. The scikit-learn User Guide recommends "sqrt" for classification and all features (1.0) for regression as defaults, to be checked by cross-validation.

The number of trees

The code records the 5-fold cross-validated accuracy of forests of increasing size, their fitting time and the time to predict the 171 validation samples.

import time

for B in [1, 10, 50, 100, 200, 500, 1000]:
    forest = RandomForestClassifier(n_estimators=B, random_state=0)
    cv_acc = cross_val_score(forest, X_train, y_train, cv=5).mean()
    start = time.perf_counter()
    forest.fit(X_train, y_train)
    fit_s = time.perf_counter() - start
    start = time.perf_counter()
    forest.predict(X_val)
    pred_ms = 1000 * (time.perf_counter() - start)
    print(B, round(cv_acc, 3), f"{fit_s:.3f} s", f"{pred_ms:.1f} ms")

The accuracy rises steeply up to about 50 trees and then stays between 0.960 and 0.965, a difference of two training samples. It does not decrease as trees are added: Breiman (2001) proved that the error of a forest converges to a limit as the number of trees grows, so more trees do not cause overfitting. The cost, however, grows in proportion to the number of trees, so n_estimators should be just large enough for the score to settle.

The trees are independent, so n_jobs=-1 fits them and predicts with them in parallel on all processor cores, without changing the result. On a synthetic dataset of 20,000 samples and 40 features, n_jobs=-1 makes fitting 200 trees more than five times faster; on the small breast cancer data, the overhead of parallelism cancels most of the gain.

Out-of-bag evaluation

Each OOB prediction uses only trees trained without the sample, so the OOB score estimates performance on new data without a validation split or refitting.

forest = RandomForestClassifier(oob_score=True, random_state=0)
forest.fit(X_train, y_train)
print(round(forest.oob_score_, 3))                # 0.97
print(forest.oob_decision_function_[:2].round(3))
# [[0.033 0.967]
#  [0.03  0.97 ]]

in_bag = np.zeros((100, len(X_train)), dtype=bool)
for b, idx in enumerate(forest.estimators_samples_):
    in_bag[b, idx] = True
oob_share = (~in_bag).mean(axis=0)  # share of trees, per sample
print(round(oob_share.mean(), 3), oob_share.min(), oob_share.max())
# 0.369 0.21 0.54

cv = cross_val_score(RandomForestClassifier(random_state=0),
                     X_train, y_train, cv=5)
print(round(cv.mean(), 3))                        # 0.96
print(round(forest.score(X_val, y_val), 3))       # 0.953

oob_decision_function_ holds the OOB class probabilities and estimators_samples_ the indices drawn for each tree. A sample is out-of-bag for 36.9% of the trees on average, as the bootstrap calculation predicts. The OOB accuracy of 0.970 is close to the cross-validated accuracy of 0.960 and the validation accuracy of 0.953, but needs one fit instead of five.

The OOB score is accuracy for a classifier and $R^2$ for a regressor; oob_score also accepts a function metric(y_true, y_pred) for the task's metric. OOB estimates require bootstrap=True, and with few trees they are noisy, because each prediction averages only about a third of the trees.

Feature importances

The impurity-based importance (MDI) of a tree is defined in Decision Trees, and a forest's feature_importances_ is the mean over its trees. Random feature subsets spread importance across correlated features: the single tree in that module gives worst perimeter 0.8 and worst radius and worst area 0, whereas the forest gives them 0.171, 0.128 and 0.066. MDI is still computed from training data and still favours features with many distinct values.

from sklearn.inspection import permutation_importance

names = list(cancer.feature_names)
mdi = forest.feature_importances_
per_tree = [t.feature_importances_ for t in forest.estimators_]
print(np.allclose(mdi, np.mean(per_tree, axis=0)))  # True
for name in ["worst perimeter", "worst radius", "worst area"]:
    print(name, round(mdi[names.index(name)], 3))
# worst perimeter 0.171
# worst radius 0.128
# worst area 0.066

result = permutation_importance(forest, X_val, y_val,
                                n_repeats=10, random_state=0)
for j in result.importances_mean.argsort()[::-1][:5]:
    print(names[j], round(result.importances_mean[j], 3),
          round(result.importances_std[j], 3))

worst radius 0.018 0.005
worst concave points 0.017 0.01
worst area 0.016 0.006
worst texture 0.013 0.005
mean texture 0.009 0.004

Shuffling worst radius lowers the validation accuracy by 0.018 on average, about 3 of the 171 samples, with a standard deviation of 0.005 over the 10 shuffles; the scoring argument selects another metric. Because the drops are measured on validation data, a feature that only helped to fit noise in the training set receives an importance near 0.

All five drops are small, although the forest's validation accuracy is 0.953. Many breast cancer features are strongly correlated, and when one is shuffled the forest obtains similar information from the others, so correlated features can all appear unimportant; the scikit-learn example on multicollinear features in the resources shows this on the same dataset.

In scikit-learn

RandomForestClassifier and RandomForestRegressor in sklearn.ensemble accept the hyperparameters of a single tree, such as criterion, max_depth and min_samples_leaf, with the same defaults, which grow full trees; the API documentation warns that such trees can be very large on some datasets. The table lists max_features, whose default differs from that of a single tree, and the hyperparameters specific to forests.

predict_proba returns the mean of the trees' class probabilities and predict the most probable class; in Breiman's original method, as the User Guide notes, each tree instead votes for one class.

proba = np.mean([t.predict_proba(X_val) for t in forest.estimators_],
                axis=0)
print(np.allclose(forest.predict_proba(X_val), proba))  # True

RandomForestRegressor averages its trees' predictions, each a mean of training labels, so like a regression tree it never predicts outside the range of its training labels: a forest fitted to $y = 2x$ for $x$ between 0 and 10 predicts 19.88 at both $x = 20$ and $x = 100$. Forests need no feature scaling, and in scikit-learn 1.9.1 they accept NaN values in X.

Extremely randomised trees

ExtraTreesClassifier and ExtraTreesRegressor add further randomness. Each split considers a random subset of features, but for each candidate feature it draws a threshold at random instead of searching, and it uses the best of these random splits. By default each tree sees the whole training set (bootstrap=False). According to the User Guide, this usually reduces variance a little more than a random forest, at the cost of slightly more bias.

from sklearn.ensemble import ExtraTreesClassifier

for model in [RandomForestClassifier(random_state=0),
              ExtraTreesClassifier(random_state=0)]:
    scores = cross_val_score(model, X_train, y_train, cv=5)
    print(type(model).__name__, round(scores.mean(), 3))
# RandomForestClassifier 0.96
# ExtraTreesClassifier 0.975

Random forests in olympiad tasks

A random forest is a strong first model for tabular data and for fixed-length feature vectors, such as embeddings from a pretrained network. It needs no scaling, captures interactions and performs reasonably with its defaults, so it gives a reference score soon after the baseline.

• Compare it with other models. In Approaching a Typical Problem, a default forest reached a validation balanced accuracy of 0.944, against 0.991 for logistic regression. Gradient boosting, covered in XGBoost, is the main alternative among tree ensembles.
• Tune little. The User Guide names n_estimators and max_features as the main parameters to adjust. The OOB score evaluates each setting with one fit; Hyperparameter Tuning covers systematic searches.
• Respect the time limit. Measure fitting and prediction time early, as advised in Reproducibility and time limits. Fewer trees, a smaller max_depth or max_samples, a larger min_samples_leaf and n_jobs=-1 all reduce it.
• Fix the seed. With random_state set, the graded notebook rebuilds the same forest, as shown in Scikit-learn Basics.