Gradient Boosting from First Principles

Bagging fits each member without using the other members’ predictions, then averages them. Boosting builds an additive model in sequence: each new tree is fitted using the current ensemble’s predictions and the chosen loss. That lets later trees address patterns the current model has missed. Both methods can affect bias and variance; which predicts better depends on the data and fitting choices.

Boosting is gradient descent in function space

Ordinary gradient descent updates parameters: \(\theta\leftarrow\theta-\eta\nabla_\theta L\). Gradient boosting instead builds a prediction function \(F\) by adding small trees. For the empirical loss \(L=\sum_i\ell(y_i,F(x_i))\), compute the negative derivative with respect to each training prediction: \(r_i=-\partial\ell(y_i,F(x_i))/\partial F(x_i)\). These \(n\) numbers indicate which way to change the current predictions to reduce loss locally.

Fit a regression tree \(h\) to the targets \(r_i\), then update \(F(x)\leftarrow F(x)+\eta h(x)\). The tree approximates the negative-gradient vector at the training inputs and supplies predictions at new inputs. This is a restricted descent direction, not the unconstrained gradient itself. The learning rate \(\eta\) scales the tree’s contribution; the actual change also depends on the value predicted by that tree.

For half-squared loss the curvature is one, so gradient and Newton targets and weights coincide when the tree-fitting and step conventions match. Other quadratic losses can differ by a constant scale. Non-quadratic losses can have varying, zero, or undefined curvature; not every loss admits this weighted least-squares Newton construction.

import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split

X, y = make_friedman1(n_samples=1000, noise=1.0, random_state=0)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)
lr, M, depth = 0.1, 100, 3

F = np.full(len(y_tr), y_tr.mean())            # squared loss: the optimal constant
F_te = np.full(len(y_te), y_tr.mean())
tree_rng = np.random.RandomState(0)
for _ in range(M):
    residual = y_tr - F                        # the negative gradient
    t = DecisionTreeRegressor(max_depth=depth, random_state=tree_rng).fit(X_tr, residual)
    F += lr * t.predict(X_tr)
    F_te += lr * t.predict(X_te)

sk = GradientBoostingRegressor(n_estimators=M, learning_rate=lr,
                               max_depth=depth, random_state=0).fit(X_tr, y_tr)
print(f"from scratch test MSE {np.mean((F_te - y_te) ** 2):.6f}")
print(f"sklearn      test MSE {np.mean((sk.predict(X_te) - y_te) ** 2):.6f}")
np.testing.assert_allclose(F_te, sk.predict(X_te), rtol=1e-10, atol=1e-10)
# from scratch test MSE 2.284947
# sklearn      test MSE 2.284947

The two test MSE values agree to the displayed precision, and the code checks the prediction vectors directly. Using one shared random-number generator across the hand-built trees matches the library’s continuing random stream in this example; resetting the tree seed to zero at every round would reuse the same initial random stream instead. For squared loss, the residual mean already minimizes loss within a fixed leaf. It does not require a different loss-specific leaf correction. Matching this example is a check of the implementation here, not a promise of identical numerical results across all versions or settings.

Shrinkage and the number of trees

The learning rate scales each added tree. A smaller value often needs more rounds to achieve a comparable training fit, and can change the entire sequence of fitted trees. It can help generalization, but can also underfit within a fixed round budget. Learning rate and the number of rounds should therefore be considered together.

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

print(f"{'lr':>6} {'M=50':>9} {'M=200':>9} {'M=800':>9} {'M=3200':>9}")
for lr in (0.5, 0.1, 0.03, 0.01):
    row = []
    for M in (50, 200, 800, 3200):
        g = GradientBoostingRegressor(n_estimators=M, learning_rate=lr,
                                      max_depth=3, random_state=0).fit(X_tr, y_tr)
        row.append(np.mean((g.predict(X_te) - y_te) ** 2))
    print(f"{lr:6.2f} " + " ".join(f"{v:9.3f}" for v in row))
#     lr      M=50     M=200     M=800    M=3200
#   0.50     3.185     3.281     3.285     3.285
#   0.10     2.938     1.985     1.956     1.983
#   0.03     8.388     2.630     1.986     1.964
#   0.01    16.054     6.635     2.448     2.048

At learning rate 0.5, test MSE stays near 3.285 at the larger round counts shown. At 0.1, the best recorded value is about 1.956. These observations do not establish that steps overshot an optimum or that recovery is impossible. The algorithms follow different sequences of trees, and test error is not their optimization objective. At 0.01, error continues to improve through 3,200 rounds, but this does not show that it must eventually reach the same value as another rate.

The settings (0.1, 200), (0.03, 800), and (0.01, 3200) produce similar test errors, but their learning-rate–round products are 20, 24, and 32, not a conserved quantity. That product is a rough way to compare fitting budgets, not an equivalence between models. One practical search tests a few learning rates with generous round limits and selects the stopping point using validation data. The displayed test set should not be reused as the final assessment after selecting among these settings.

At learning rate 0.1, test MSE rises from 1.956 at 800 trees to 1.983 at 3,200. This small increase on one test set is consistent with diminishing generalization benefit; it does not establish that only noise remains in the residuals. Boosting rounds change the fitted function by sequential optimization. Forest tree counts primarily change how accurately an average is approximated, though an individual forest test score can also fluctuate as trees are added.

Second-order boosting

For binary log loss, \(F\) is a raw score interpreted as log odds and \(p=1/(1+e^{-F})\) is the class-one probability. Define \(g=y-p\), the negative gradient with respect to \(F\), and \(h=p(1-p)\), the second derivative. A second-order expansion for a correction \(u\) is \(\ell(F+u)\approx\ell(F)-gu+\tfrac12hu^2\). With \(h>0\), minimizing the summed approximation is equivalent to fitting targets \(g/h\) with weights \(h\). A leaf’s unregularized correction is therefore \(\sum_i g_i/\sum_i h_i\), not the unweighted mean of the ratios. This is the principle behind Newton-style tree updates; production libraries add further split rules and regularization.

import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from scipy.special import expit

Xc, yc = make_classification(n_samples=3000, n_features=10, n_informative=6,
                             random_state=0)
A_tr, A_te, b_tr, b_te = train_test_split(Xc, yc, test_size=0.3, random_state=0)

def boost(newton, M=60, lr=0.1, depth=3):
    F, F_te, curve = np.zeros(len(b_tr)), np.zeros(len(b_te)), []
    for _ in range(M):
        p = expit(F)
        g = b_tr - p                           # negative gradient
        if newton:
            h = np.clip(p * (1 - p), 1e-6, None)   # second derivative
            target, weight = g / h, h
        else:
            target, weight = g, None
        t = DecisionTreeRegressor(max_depth=depth, random_state=0).fit(
            A_tr, target, sample_weight=weight)
        F += lr * t.predict(A_tr)
        F_te += lr * t.predict(A_te)
        curve.append(np.mean(np.logaddexp(0.0, F_te) - b_te * F_te))
    return curve

first, newton = boost(False), boost(True)
for m in (10, 20, 40, 60):
    print(f"after {m:3d} trees   first-order {first[m-1]:.4f}   Newton {newton[m-1]:.4f}")
# after  10 trees   first-order 0.5862   Newton 0.4210
# after  20 trees   first-order 0.5180   Newton 0.3301
# after  40 trees   first-order 0.4386   Newton 0.2650
# after  60 trees   first-order 0.3934   Newton 0.2420

At 60 rounds, the Newton-style version reaches log loss about 0.242 and the first-order version about 0.393. They use the same numerical learning rate, but their step magnitudes differ, so this is not a comparison of separately tuned methods or a measurement of how many times faster one converges. Both start at \(F=0\), or probability 0.5. At that starting point \(h=0.25\), making each Newton target four times the first-order target.

The label matters when interpreting curvature. If \(y=1\), then \(p=0.99\) gives \(g=0.01\) and \(g/h\approx1.01\), while \(p=0.6\) gives \(g=0.4\) and \(g/h\approx1.67\). For a confident wrong prediction with \(y=0,p=0.99\), the ratio is −100. Small curvature can produce large corrections, so safeguards matter. The code floors curvature at \(10^{-6}\), changing the quadratic approximation when that floor is active, and uses a learning rate below one. Those choices alone are not a universal stability guarantee.

For half-squared loss, \(h=1\) everywhere, so the two constructions coincide when tree fitting and step conventions match. The logistic example instead has prediction-dependent curvature, which changes both the targets and the fitting weights. Other losses require their own curvature analysis; zero or undefined curvature does not support this Newton weighted least-squares construction.

Early stopping

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split

Xf, yf = make_friedman1(n_samples=800, noise=2.0, random_state=0)
P_tr, P_te, q_tr, q_te = train_test_split(Xf, yf, test_size=0.3, random_state=0)

stopped = GradientBoostingRegressor(n_estimators=2000, learning_rate=0.05,
                                    max_depth=3, validation_fraction=0.2,
                                    n_iter_no_change=20, random_state=0).fit(P_tr, q_tr)
full = GradientBoostingRegressor(n_estimators=2000, learning_rate=0.05,
                                 max_depth=3, random_state=0).fit(P_tr, q_tr)

print(f"early stopping: {stopped.n_estimators_:4d} trees"
      f"   test MSE {np.mean((stopped.predict(P_te) - q_te) ** 2):.4f}")
print(f"all trees:      {full.n_estimators_:4d} trees"
      f"   test MSE {np.mean((full.predict(P_te) - q_te) ** 2):.4f}")
refit = GradientBoostingRegressor(n_estimators=stopped.n_estimators_,
                                  learning_rate=0.05, max_depth=3,
                                  random_state=0).fit(P_tr, q_tr)
print(f"refit selected: {refit.n_estimators_:4d} trees"
      f"   test MSE {np.mean((refit.predict(P_te) - q_te) ** 2):.4f}")
# early stopping:  237 trees   test MSE 6.1170
# all trees:      2000 trees   test MSE 6.1469
# refit selected:  237 trees   test MSE 6.0851

The stopped estimator reserves 20% of the 560 outer training rows for its internal validation, leaving 448 rows to fit its trees. The 2,000-round estimator uses all 560. Their test-score comparison therefore changes both the stopping rule and the fitting sample size. The additional refit row uses the selected stopping count on all 560 training rows. Here test MSE is 6.1170 for the stopped model, 6.1469 for the full 2,000-round model, and 6.0851 for the 237-round refit. This is a comparison of those procedures on one test set, not a guarantee that refitting improves performance.

Early stopping can avoid building unnecessary trees, but the reduction in elapsed time is not exactly the ratio of tree counts: validation prediction and other work also cost time. If the stopping rule reaches the round limit, it may save no trees. The stored stopping count also need not be the round with the smallest validation loss; it reflects the estimator’s tolerance and patience rule.

The stopping split must match the prediction task. A random internal split is unsuitable for evaluating future-time or new-group generalization. This GradientBoostingRegressor API does not accept an eval_set argument. For an explicit split, fit a model without internal stopping and evaluate staged_predict on your validation set, then choose a round count and refit. That approach selects model size but does not save the initial full fit; online stopping needs an appropriate training loop or estimator API. Keep the final test set outside all stopping and hyperparameter choices.

References

• Friedman, J. H. (2001). Greedy Function Approximation: A Gradient Boosting Machine. Annals of Statistics.
• Friedman, J. H. (2002). Stochastic Gradient Boosting. Computational Statistics & Data Analysis.
• scikit-learn documentation. sklearn.ensemble.GradientBoostingRegressor
• InterpretML. Explainable Boosting Machine
• XGBoost documentation. Parameters: column subsampling