Gradient Boosting, in depth: functional gradient descent, the second-order split gain, histograms, and XGBoost vs LightGBM vs CatBoost

Boosting as gradient descent in function space

Boosting builds a model as a sum of weak learners, added one at a time. After m rounds the prediction for a row x is F_m(x) = F_0 + eta x (f_1(x) + ... + f_m(x)), where each f is a small tree and eta, the learning rate or shrinkage, is typically between 0.01 and 0.3. Each new tree is chosen to reduce the training loss of the ensemble so far.

Here is the key step. Treat the vector of current predictions, one per training row, as the parameters. The loss L = sum of l(y_i, F(x_i)) has a gradient with respect to each prediction, g_i = dl/dF at F(x_i). Gradient descent would move each prediction by minus eta times g_i. We cannot do that directly, because the model must generalise to rows it has not seen, so instead we fit a tree to the negative gradients and step along the tree. That is why it is called gradient boosting: it is gradient descent in function space, with each step constrained to be a tree.

For squared error the negative gradient is the residual, which gives the familiar story: each tree fits what the previous trees got wrong. For other losses the residual intuition breaks, and the gradient view is the one that generalises. Absolute error shows why the Hessian matters next: it is zero, so second-order methods need a workaround for it.

The second-order objective and the gain formula

XGBoost's contribution was to take a second-order Taylor expansion of the loss around the current predictions and add explicit regularisation on the tree. For a tree with leaves j, leaf weights w_j and leaf row sets I_j, write G_j and H_j for the sums of g_i and h_i over rows in leaf j. The objective for the new tree is approximately the sum over leaves of G_j w_j + (H_j + lambda) w_j^2 / 2, plus gamma per leaf. Minimising each leaf independently gives a closed form:

optimal leaf weight:   w_j* = -G_j / (H_j + lambda)
leaf score:            -1/2 x G_j^2 / (H_j + lambda)
split gain:            1/2 x [ G_L^2/(H_L+lambda) + G_R^2/(H_R+lambda) - (G_L+G_R)^2/(H_L+H_R+lambda) ] - gamma

Every hyperparameter in that formula has a direct meaning. lambda (L2 on leaf weights) shrinks weights toward zero, more strongly in leaves with little Hessian mass. gamma, called min_split_loss in XGBoost, is the minimum gain a split must deliver, so it prunes weak splits. min_child_weight is a lower bound on H in each child; for squared error H is just a row count, while for logistic loss it is the sum of p(1-p), so leaves full of confident predictions carry little weight. Split finding is the same scan for every loss: sort or bin a feature, sweep left to right accumulating G_L and H_L, and evaluate the gain at each threshold.

Worked example: one split by hand

Four rows, one feature, squared error, lambda = 1, gamma = 0, eta = 0.3. Targets y = 1, 2, 6 and 7 at feature values x = 1, 2, 3 and 4. Start from the mean, F_0 = 4. The gradients are g = F - y = 3, 2, -2, -3 and every h = 1. At the root G = 0 and H = 4, so the root score term is 0.

The middle split wins. Leaf weights are -5/3 = -1.667 on the left and +1.667 on the right. After shrinkage by 0.3 the left rows move from 4 to 3.5 and the right rows to 4.5. The squared loss falls from 13 to 8.5, a deliberate under-step: the leaf weights are already shrunk by lambda, from the unregularised -2.5 toward zero, and eta damps them again so later trees can correct the structure instead of fighting an over-committed first tree.

Code: a gradient booster in NumPy

The booster below implements exactly those formulas, for squared and logistic loss, with exact greedy splits. It is slow and has no missing-value handling, but every line maps to the derivation, and on a small dataset its predictions track a library with the same settings closely.

import numpy as np

def grad_hess(y, F, loss):
    if loss == "squared":
        return F - y, np.ones_like(F)
    p = 1.0 / (1.0 + np.exp(-F))                     # logistic
    return p - y, p * (1.0 - p)

def build(X, g, h, depth, lam, gamma, min_h):
    G, H = g.sum(), h.sum()
    leaf = {"w": -G / (H + lam)}
    if depth == 0 or H < 2 * min_h:
        return leaf
    best = (0.0, None)
    parent = G * G / (H + lam)
    for f in range(X.shape[1]):
        order = np.argsort(X[:, f])
        gl, hl = np.cumsum(g[order])[:-1], np.cumsum(h[order])[:-1]
        gr, hr = G - gl, H - hl
        gain = 0.5 * (gl**2 / (hl + lam) + gr**2 / (hr + lam) - parent) - gamma
        gain[(hl < min_h) | (hr < min_h)] = -np.inf
        xs = X[order, f]
        gain[xs[:-1] == xs[1:]] = -np.inf             # no split between equal values
        i = int(np.argmax(gain))
        if gain[i] > best[0]:
            best = (gain[i], (f, (xs[i] + xs[i + 1]) / 2))
    if best[1] is None:
        return leaf
    f, thr = best[1]
    m = X[:, f] < thr
    return {"f": f, "thr": thr,
            "l": build(X[m], g[m], h[m], depth - 1, lam, gamma, min_h),
            "r": build(X[~m], g[~m], h[~m], depth - 1, lam, gamma, min_h)}

def predict_tree(t, X):
    if "f" not in t:
        return np.full(len(X), t["w"])
    m = X[:, t["f"]] < t["thr"]
    out = np.empty(len(X))
    out[m], out[~m] = predict_tree(t["l"], X[m]), predict_tree(t["r"], X[~m])
    return out

def boost(X, y, loss="squared", rounds=200, eta=0.1, depth=3, lam=1.0, gamma=0.0, min_h=1.0):
    F0 = y.mean() if loss == "squared" else np.log(y.mean() / (1 - y.mean()))
    F, trees = np.full(len(y), F0), []
    for _ in range(rounds):
        g, h = grad_hess(y, F, loss)
        t = build(X, g, h, depth, lam, gamma, min_h)
        F += eta * predict_tree(t, X)
        trees.append(t)
    return F0, trees

Run it on the four-row example with rounds=1, eta=0.3, depth=1 and you get predictions of 3.5 and 4.5, matching the hand calculation.

Histograms, tree growth and categories at scale

Exact greedy splitting sorts every feature at every node, which is too slow for millions of rows. All three libraries instead bin each feature into a fixed number of quantile buckets once (XGBoost's hist method uses 256 by default, LightGBM 255), then accumulate G and H per bin. Finding a split becomes a scan over a few hundred bins instead of millions of rows, and a node's histogram can be computed as its parent's minus its sibling's, halving the work. Histogram building is also what GPUs accelerate: in XGBoost, device="cuda" with tree_method="hist" moves it onto the GPU.

The libraries differ mainly in how they grow trees and handle categories. XGBoost grows level-wise by default, splitting every node at a depth, which gives balanced trees and predictable cost. LightGBM grows leaf-wise, always splitting the single leaf with the highest gain, capped by num_leaves; it reaches lower loss with fewer leaves but overfits small data unless min_data_in_leaf is raised. LightGBM also adds GOSS, which keeps rows with large gradients and samples the rest, and exclusive feature bundling for sparse data. CatBoost uses oblivious trees, where every node at a depth shares one split, which regularises strongly and makes inference very fast. It encodes categorical features with ordered target statistics, computing each row's encoding only from rows before it in a random permutation, which prevents the target leakage that naive mean encoding causes.

Using XGBoost and LightGBM, and how to tune them

In practice you call a library. The pattern is the same in each: hold out a validation set, set a low learning rate and a generous round limit, and let early stopping pick the number of trees.

import xgboost as xgb, lightgbm as lgb

model = xgb.XGBClassifier(
    n_estimators=5000, learning_rate=0.05, max_depth=6,
    min_child_weight=5, subsample=0.8, colsample_bytree=0.8,
    reg_lambda=1.0, tree_method="hist", device="cuda",
    eval_metric="logloss", early_stopping_rounds=100,
    enable_categorical=True,          # pandas category dtype columns
)
model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=False)
print(model.best_iteration)

lgbm = lgb.LGBMClassifier(n_estimators=5000, learning_rate=0.05, num_leaves=63,
                          min_child_samples=50, subsample=0.8, subsample_freq=1,
                          colsample_bytree=0.8, reg_lambda=1.0)
lgbm.fit(X_train, y_train, eval_set=[(X_valid, y_valid)],
         callbacks=[lgb.early_stopping(100)])

Tune in this order. Fix the learning rate at 0.05 to 0.1 and let early stopping choose the rounds. Then tune tree size, max_depth or num_leaves, together with min_child_weight or min_child_samples, since these interact most with overfitting. Then row and column subsampling, then lambda and gamma. Finally lower the learning rate, raise the round limit, and refit for a small extra gain. Use cross-validation for the search, and keep a test set that no early-stopping or tuning decision ever sees. Domain knowledge can be encoded too: monotone_constraints forces predictions to rise or fall with a feature, which regulators and users often expect for features such as income in credit models.

Failure modes

Overfitting in boosting is gradual: each tree is small, so validation loss usually bottoms out and rises slowly. That makes early stopping effective, but only if the validation data resembles production; see overfitting for the general diagnosis.

Trade-offs

Boosting also has a structural cost: training is sequential across rounds, so it parallelises within a tree, not across trees. Inference is a sum of hundreds of tree traversals, which is fast on CPU but worth benchmarking at high request rates.

What to do next

1. Derive the gain formula yourself from the second-order objective, then verify the worked example numbers with the NumPy booster.
2. Train XGBoost and LightGBM on one of your own tabular datasets with a time-based split, early stopping on validation only, and an untouched test set.
3. Tune in the order given: learning rate and rounds, tree size with leaf minimums, sampling, then regularisation.
4. Audit every target-derived feature for leakage, and compare against CatBoost with the raw categorical columns.
5. Check calibration and per-segment error before shipping, and add monotone constraints where the domain demands them.