A Binary Classification GBM in C++

July 2026 · 18 min read

The regression GBM predicted a number directly. A classification GBM should output a probability. The standard approach is to make the additive model predict a logit $F(x)$, then convert it to a probability with a sigmoid:

The trees still add corrections to $F(x)$. The probability is derived only when evaluating the loss or making a final prediction.

Pseudo-Residuals for Log Loss

With labels $y_i \in \{0,1\}$, binary log loss is

The derivative with respect to the logit is $p_i-y_i$, so the negative gradient is $y_i-p_i$. This becomes the target for the next regression tree.

Sigmoid and Metrics

double sigmoid(double z) {
    if (z >= 0) {
        double e = exp(-z);
        return 1.0 / (1.0 + e);
    }
    double e = exp(z);
    return e / (1.0 + e);
}

double clip01(double p) {
    return min(1.0 - 1e-12, max(1e-12, p));
}

The numerically stable sigmoid avoids overflow for large positive or negative logits. Clipping probabilities keeps log loss finite during evaluation.

The Logistic GBM Class

This implementation uses a first-order logistic GBM: it fits shallow regression trees to $y-p$ and adds their predictions to the logit with shrinkage.

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

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

    void fit(const Matrix& X, const vector<int>& y) {
        int n = (int)X.size();
        double pos = accumulate(y.begin(), y.end(), 0.0);
        double rate = min(1.0 - 1e-6, max(1e-6, pos / n));
        F0 = log(rate / (1.0 - rate));              // best constant logit

        vector<double> logit(n, F0);
        trees.clear();
        trees.reserve(rounds);

        for (int m = 0; m < rounds; ++m) {
            vector<double> pseudo(n);
            for (int i = 0; i < n; ++i) {
                double p = sigmoid(logit[i]);
                pseudo[i] = y[i] - p;               // negative gradient
            }

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

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

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

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

    double predict_proba(const vector<double>& x) const {
        return sigmoid(predict_logit(x));
    }

    int predict(const vector<double>& x) const {
        return predict_proba(x) >= 0.5 ? 1 : 0;
    }
};

Accuracy and Log Loss

double accuracy(const BinaryGBM& model,
                const Matrix& X,
                const vector<int>& y) {
    int ok = 0;
    for (size_t i = 0; i < X.size(); ++i)
        ok += model.predict(X[i]) == y[i];
    return 100.0 * ok / X.size();
}

double log_loss(const BinaryGBM& model,
                const Matrix& X,
                const vector<int>& y) {
    double s = 0.0;
    for (size_t i = 0; i < X.size(); ++i) {
        double p = clip01(model.predict_proba(X[i]));
        s += -(y[i] * log(p) + (1 - y[i]) * log(1 - p));
    }
    return s / X.size();
}

Ring Versus Center Dataset

The test dataset labels points near the origin as class 0 and points in a surrounding ring as class 1. Linear models struggle with this boundary; shallow trees can carve it in pieces, and boosting improves the boundary over rounds.

int label(double x, double y) {
    double r2 = x * x + y * y;
    return r2 > 2.0 ? 1 : 0;
}

int main() {
    mt19937 rng(7);
    uniform_real_distribution<double> u(-2.2, 2.2);
    normal_distribution<double> noise(0.0, 0.08);

    Matrix X;
    vector<int> y;
    for (int i = 0; i < 500; ++i) {
        double a = u(rng), b = u(rng);
        X.push_back({a + noise(rng), b + noise(rng)});
        y.push_back(label(a, b));
    }

    Matrix Xtr(X.begin(), X.begin() + 380), Xte(X.begin() + 380, X.end());
    vector<int> ytr(y.begin(), y.begin() + 380), yte(y.begin() + 380, y.end());

    for (int rounds : {1, 5, 20, 80}) {
        BinaryGBM gbm(rounds, 3, 8, 0.12);
        gbm.fit(Xtr, ytr);
        printf("trees=%2d  test acc=%.1f%%  test logloss=%.3f\\n",
               rounds, accuracy(gbm, Xte, yte), log_loss(gbm, Xte, yte));
    }
}

Implementation Notes

This first-order version is intentionally small. Production libraries usually use Newton leaf values with both gradient and Hessian sums. That improves step sizes, especially for logistic loss. The conceptual StatQuest classification derivation post derives the second-order formula; the important implementation point here is that the additive-tree skeleton is unchanged.