From Scratch · C++ · ML · Gradient Boosted Machines · Part 1

A Regression GBM in C++

July 2026 · 18 min read

This post implements the core Gradient Boosted Machine for regression. The base learner is the squared-error CART regression tree from Regression Tree (CART). The new code is the boosting wrapper: initialize with the target mean, compute residuals, fit a shallow tree to those residuals, then add a shrunken tree prediction to the ensemble.

Data Types

The implementation keeps the same data layout used by the decision-tree posts:

#include <bits/stdc++.h>
using namespace std;

using Matrix = vector<vector<double>>;

The existing RegTree exposes two methods:

struct RegTree {
    RegTree(int maxDepth, int minSamples);
    void fit(const Matrix& X, const vector<double>& y);
    double predict(const vector<double>& x) const;
};

Internally it uses the sorted squared-error split sweep from the CART post. GBM does not need a new tree algorithm; it needs a new target for each tree.

The GBM Struct

The ensemble stores an intercept F0, a learning rate, and the fitted trees in the exact order they were trained:

struct RegressionGBM {
    int rounds;
    int treeDepth;
    int minSamples;
    double learningRate;
    double F0 = 0.0;
    vector<RegTree> trees;

    RegressionGBM(int rounds_, int depth_, int minSamples_, double lr_)
        : rounds(rounds_), treeDepth(depth_), minSamples(minSamples_),
          learningRate(lr_) {}

    void fit(const Matrix& X, const vector<double>& y) {
        int n = (int)X.size();
        F0 = accumulate(y.begin(), y.end(), 0.0) / n;
        vector<double> pred(n, F0);

        trees.clear();
        trees.reserve(rounds);

        for (int m = 0; m < rounds; ++m) {
            vector<double> residual(n);
            for (int i = 0; i < n; ++i)
                residual[i] = y[i] - pred[i];

            RegTree tree(treeDepth, minSamples);
            tree.fit(X, residual);

            for (int i = 0; i < n; ++i)
                pred[i] += learningRate * tree.predict(X[i]);

            trees.push_back(move(tree));
        }
    }

    double predict(const vector<double>& x) const {
        double out = F0;
        for (const auto& tree : trees)
            out += learningRate * tree.predict(x);
        return out;
    }
};

The residual is the gradient target. For squared loss, the negative gradient is exactly y[i] - pred[i]. That is why this first GBM looks like repeated residual fitting.

Evaluation Helper

Regression quality is reported as root mean squared error:

double rmse(const RegressionGBM& model,
            const Matrix& X,
            const vector<double>& y) {
    double s = 0.0;
    for (size_t i = 0; i < X.size(); ++i) {
        double e = model.predict(X[i]) - y[i];
        s += e * e;
    }
    return sqrt(s / X.size());
}

Noisy Nonlinear Dataset

A quick test generates a two-feature regression problem with a smooth nonlinear target:

int main() {
    mt19937 rng(42);
    uniform_real_distribution<double> ux(-3.0, 3.0);
    normal_distribution<double> noise(0.0, 0.15);

    Matrix X;
    vector<double> y;
    for (int i = 0; i < 320; ++i) {
        double a = ux(rng), b = ux(rng);
        X.push_back({a, b});
        y.push_back(sin(a) + 0.35 * b * b + noise(rng));
    }

    vector<int> order(X.size());
    iota(order.begin(), order.end(), 0);
    shuffle(order.begin(), order.end(), rng);

    Matrix Xtr, Xte;
    vector<double> ytr, yte;
    for (size_t k = 0; k < order.size(); ++k) {
        int i = order[k];
        if (k < 240) { Xtr.push_back(X[i]); ytr.push_back(y[i]); }
        else         { Xte.push_back(X[i]); yte.push_back(y[i]); }
    }

    for (int rounds : {1, 5, 20, 80}) {
        RegressionGBM gbm(rounds, 3, 6, 0.08);
        gbm.fit(Xtr, ytr);
        printf("trees=%2d  train RMSE=%.3f  test RMSE=%.3f\\n",
               rounds, rmse(gbm, Xtr, ytr), rmse(gbm, Xte, yte));
    }
}

Compiled with g++ -O2 -std=c++17, a representative run prints:

trees= 1 train RMSE=0.960 test RMSE=0.928 trees= 5 train RMSE=0.610 test RMSE=0.632 trees=20 train RMSE=0.258 test RMSE=0.291 trees=80 train RMSE=0.119 test RMSE=0.190

Reading the Numbers

One tree is too weak to capture the nonlinear target. More rounds keep lowering training error, and test error follows until the model begins spending capacity on noise. The gap between train and test at 80 trees is the first hint that regularization and early stopping matter, which the final implementation post adds explicitly.

Next Step

The same structure can classify by changing only the prediction scale and the gradient target. Instead of a raw numeric prediction, the model stores a logit. Instead of residuals $y-F$, it fits pseudo-residuals $y-p$. That is the binary classification GBM in the next post.