From Scratch · C++ · Machine Learning · Decision Trees

A Regression Tree (CART) in C++

June 2026 · 16 min read

A regression tree is the classification tree from Part 1 with two changes: leaves store a number instead of a class, and the impurity is squared error around the mean instead of Gini. Everything else, the recursion, the stopping rules, the sorted sweep, carries over. This is the CART (Classification And Regression Trees) variant. The finished program fits a noisy sine to a test RMSE of about $0.13$.

Squared-Error Impurity

For a set $S$ of targets with mean $\bar y$, the impurity is the sum of squared deviations, and the leaf prediction is that mean:

$$L(S) = \sum_{(x,y)\in S}(y - \bar y)^2, \qquad \bar y = \frac{1}{n}\sum_{(x,y)\in S} y.$$

A pure regression node is one whose targets are all near a single value. The split that most reduces total squared error wins. The key to scoring splits fast is the identity

$$\sum_i (y_i - \bar y)^2 = \sum_i y_i^2 - \frac{\big(\sum_i y_i\big)^2}{n},$$

which means $L$ can be computed from just the running sum and sum of squares. As the threshold sweeps, both update in $O(1)$:

static double sse(double sum, double sumsq, int n) {
    return n ? sumsq - sum * sum / n : 0;   // sum of squared deviations
}

The Node and the Sweep

The node now carries a value (the leaf's mean). The split search keeps a running sum and sum-of-squares for each side and never rescans the data:

#include <bits/stdc++.h>
using namespace std;
using Matrix = vector<vector<double>>;

struct RNode {
    bool leaf = false;
    int feature = -1; double threshold = 0;
    double value = 0;              // mean target, used at leaves
    RNode *L = nullptr, *R = nullptr;
};

struct RegTree {
    int maxDepth, minSamples;
    RNode* root = nullptr;
    RegTree(int md, int ms) : maxDepth(md), minSamples(ms) {}

    RNode* build(vector<int> idx, const Matrix& X, const vector<double>& y, int depth) {
        RNode* nd = new RNode();
        int n = idx.size();
        double S = 0, SS = 0;
        for (int i : idx) { S += y[i]; SS += y[i] * y[i]; }
        nd->value = n ? S / n : 0;              // leaf prediction = mean
        double parentSSE = sse(S, SS, n);

        if (depth >= maxDepth || n < minSamples || parentSSE < 1e-9) {
            nd->leaf = true; return nd;
        }

        int D = X[0].size();
        double bestGain = 1e-9; int bestFeat = -1; double bestThr = 0;
        for (int f = 0; f < D; ++f) {
            vector<int> ord = idx;
            sort(ord.begin(), ord.end(),
                 [&](int a, int b){ return X[a][f] < X[b][f]; });
            double sL = 0, ssL = 0, sR = S, ssR = SS;       // left empty, right full
            for (int t = 0; t + 1 < n; ++t) {
                int i = ord[t];
                sL += y[i]; ssL += y[i] * y[i];
                sR -= y[i]; ssR -= y[i] * y[i];
                if (X[ord[t]][f] == X[ord[t + 1]][f]) continue;
                int nl = t + 1, nr = n - nl;
                double gain = parentSSE - (sse(sL, ssL, nl) + sse(sR, ssR, nr));
                if (gain > bestGain) {
                    bestGain = gain; bestFeat = f;
                    bestThr = (X[ord[t]][f] + X[ord[t + 1]][f]) / 2;
                }
            }
        }

        if (bestFeat < 0) { nd->leaf = true; 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(const Matrix& X, const vector<double>& y) {
        vector<int> idx(X.size());
        iota(idx.begin(), idx.end(), 0);
        root = build(idx, X, y, 0);
    }
    double predict(const vector<double>& x) const {
        RNode* n = root;
        while (!n->leaf) n = (x[n->feature] < n->threshold) ? n->L : n->R;
        return n->value;
    }
};

Any loss, no gradients. Because CART simply tries every candidate split value and keeps the best, the impurity does not need to be differentiable. Squared error is the common choice, but absolute error or quantile loss would slot in the same way. The original CART paper evaluates a split at every value occurring in the data for that feature.

Running It

Fit $y = \sin(2\pi x) + \varepsilon$ with small Gaussian noise:

int main() {
    mt19937 rng(1);
    normal_distribution<double> noise(0, 0.1);
    uniform_real_distribution<double> ux(0, 1);
    Matrix X; vector<double> y; int N = 400;
    for (int i = 0; i < N; ++i) {
        double x = ux(rng);
        X.push_back({x});
        y.push_back(sin(2 * M_PI * x) + noise(rng));
    }
    int ntr = 0.8 * N;
    Matrix Xtr(X.begin(), X.begin() + ntr), Xte(X.begin() + ntr, X.end());
    vector<double> ytr(y.begin(), y.begin() + ntr), yte(y.begin() + ntr, y.end());

    RegTree t(/*maxDepth*/7, /*minSamples*/5);
    t.fit(Xtr, ytr);
    auto rmse = [&](const Matrix& Xs, const vector<double>& ys) {
        double s = 0;
        for (size_t i = 0; i < Xs.size(); ++i) {
            double d = t.predict(Xs[i]) - ys[i]; s += d * d;
        }
        return sqrt(s / Xs.size());
    };
    printf("train RMSE = %.3f   test RMSE = %.3f\n", rmse(Xtr, ytr), rmse(Xte, yte));
}

train RMSE = 0.068 test RMSE = 0.134

The tree approximates the sine as a staircase of constant pieces. With the noise standard deviation at $0.1$, a test RMSE of $0.134$ is close to the noise floor: the tree has recovered the signal and is mostly fitting the irreducible noise. A single tree's staircase is jagged, though, which is exactly the high-variance behavior that the next part smooths away by averaging.