📉 Gradient Boosting from Scratch
Gradient Boosting from Scratch
[[adaboost]] reweights the training rows. [[xgboost]], [[lightgbm]] and [[catboost]] all describe themselves as gradient boosting. This page is the step in between — the idea those three libraries are optimised implementations of.
The idea in one sentence: fit each new tree to the negative gradient of the loss with respect to the current predictions.
That sounds abstract until you compute the gradient for squared error, at which point it becomes something much more familiar:
$$L = \tfrac{1}{2}(y - F)^2 \qquad \frac{\partial L}{\partial F} = -(y - F) \qquad -\frac{\partial L}{\partial F} = y - F = \text{the residual}$$
So "fit the next tree to the residuals" is not a separate algorithm — it is what gradient boosting becomes when the loss is squared error. Change the loss and the same machinery gives you something else, which is the whole point of the generalisation.
Every number below is computed by hand on 12 EEG recordings and checked against scikit-learn.
Display helpers ready: banner, step, check, expand, show_matrix, show_matmul, show_elementwise.
12 recordings x 3 features
id region line_length spike_rate age y_class y_reg
-------------------------------------------------------------------------------
rec00 frontal 12.0 0.5 24.0 0 0.2
rec01 temporal 15.0 1.2 31.0 0 0.5
rec02 temporal 42.0 6.8 58.0 1 6.4
rec03 occipital 38.0 5.1 47.0 1 4.8
rec04 frontal 11.0 0.3 19.0 0 0.1
rec05 temporal 55.0 8.2 63.0 1 9.1
rec06 occipital 40.0 1.9 35.0 1 3.5
rec07 temporal 47.0 7.4 52.0 1 7.2
rec08 frontal 14.0 0.8 28.0 0 0.3
rec09 occipital 51.0 6.1 60.0 1 6.9
rec10 frontal 45.0 2.4 41.0 0 1.6
rec11 temporal 35.0 4.7 30.0 1 4.1
-------------------------------------------------------------------------------
class balance: 5 negative, 7 positive
y_reg mean = 3.7250, std = 3.0518
1. Three rounds, by hand
The algorithm (Friedman, 2001):
- Start with a constant prediction $F_0$ — the value minimising the loss with no features. For squared error that is the mean of $y$.
- For $m = 1 \dots M$:
- compute the pseudo-residuals $r_i = -\partial L/\partial F(x_i)$ evaluated at the current $F_{m-1}$
- fit a shallow regression tree to those residuals
- update $F_m = F_{m-1} + \eta\, h_m$, where $\eta$ is the learning rate
The trees are deliberately weak — depth 1 or 2 — and $\eta$ is small. Each round corrects a little of what remains.
==============================================================================
GRADIENT BOOSTING, THREE ROUNDS ON 12 RECORDINGS
==============================================================================
Round 0: the constant that minimises squared error
==============================================================================
F_0 = mean(y) = (0.2 + 0.5 + 6.4 + 4.8 + ...) / 12 = 3.7250
SSE = 111.7625
Round 1
==============================================================================
pseudo-residuals r = y - F (this IS -dL/dF for squared error)
id y F r = y-F
rec00 0.20 3.7250 -3.5250
rec01 0.50 3.7250 -3.2250
rec02 6.40 3.7250 2.6750
rec03 4.80 3.7250 1.0750
rec04 0.10 3.7250 -3.6250
rec05 9.10 3.7250 5.3750
rec06 3.50 3.7250 -0.2250
rec07 7.20 3.7250 3.4750
rec08 0.30 3.7250 -3.4250
rec09 6.90 3.7250 3.1750
rec10 1.60 3.7250 -2.1250
rec11 4.10 3.7250 0.3750
tree 1: split on spike_rate <= 3.550
left leaf (<= 3.550) : -2.6917
right leaf (> 3.550) : +2.6917
check: mean(r | left) = -2.6917 MATCH: True
mean(r | right) = +2.6917 MATCH: True
update F <- F + 0.5 * h
SSE now = 46.5569
Round 2
==============================================================================
pseudo-residuals r = y - F (this IS -dL/dF for squared error)
id y F r = y-F
rec00 0.20 2.3792 -2.1792
rec01 0.50 2.3792 -1.8792
rec02 6.40 5.0708 1.3292
rec03 4.80 5.0708 -0.2708
rec04 0.10 2.3792 -2.2792
rec05 9.10 5.0708 4.0292
rec06 3.50 2.3792 1.1208
rec07 7.20 5.0708 2.1292
rec08 0.30 2.3792 -2.0792
rec09 6.90 5.0708 1.8292
rec10 1.60 2.3792 -0.7792
rec11 4.10 5.0708 -0.9708
tree 2: split on age <= 49.500
left leaf (<= 49.500) : -1.1646
right leaf (> 49.500) : +2.3292
check: mean(r | left) = -1.1646 MATCH: True
mean(r | right) = +2.3292 MATCH: True
update F <- F + 0.5 * h
SSE now = 22.1443
Round 3
==============================================================================
pseudo-residuals r = y - F (this IS -dL/dF for squared error)
id y F r = y-F
rec00 0.20 1.7969 -1.5969
rec01 0.50 1.7969 -1.2969
rec02 6.40 6.2354 0.1646
rec03 4.80 4.4885 0.3115
rec04 0.10 1.7969 -1.6969
rec05 9.10 6.2354 2.8646
rec06 3.50 1.7969 1.7031
rec07 7.20 6.2354 0.9646
rec08 0.30 1.7969 -1.4969
rec09 6.90 6.2354 0.6646
rec10 1.60 1.7969 -0.1969
rec11 4.10 4.4885 -0.3885
tree 3: split on line_length <= 36.500
left leaf (<= 36.500) : -1.2952
right leaf (> 36.500) : +0.9251
check: mean(r | left) = -1.2952 MATCH: True
mean(r | right) = +0.9251 MATCH: True
update F <- F + 0.5 * h
SSE now = 11.3600
Three rounds of progress
==============================================================================
round SSE reduction
------------------------------
0 111.7625
1 46.5569 65.2056
2 22.1443 24.4126
3 11.3600 10.7843
------------------------------
2. Checking it against scikit-learn
The hand-built model above should be exactly what GradientBoostingRegressor
produces with the same settings. Worth confirming, because it means the mental
model is the implementation, not an analogy for it.
==============================================================================
HAND-BUILT vs SCIKIT-LEARN
==============================================================================
id y by hand sklearn diff
------------------------------------------------
rec00 0.20 1.1493 1.1493 0.00e+00
rec01 0.50 1.1493 1.1493 0.00e+00
rec02 6.40 6.6980 6.6980 0.00e+00
rec03 4.80 4.9511 4.9511 0.00e+00
rec04 0.10 1.1493 1.1493 0.00e+00
rec05 9.10 6.6980 6.6980 0.00e+00
rec06 3.50 2.2594 2.2594 -4.44e-16
rec07 7.20 6.6980 6.6980 0.00e+00
rec08 0.30 1.1493 1.1493 0.00e+00
rec09 6.90 6.6980 6.6980 0.00e+00
rec10 1.60 2.2594 2.2594 -4.44e-16
rec11 4.10 3.8409 3.8409 0.00e+00
------------------------------------------------
max|diff| = 4.44e-16 MATCH: True
init prediction: mine 3.725000 | sklearn 3.725000 | MATCH: True
3. The learning rate: why shrink at all?
Nothing forces $\eta < 1$. With $\eta = 1$ each tree is applied at full strength and the training error falls fastest. So why shrink?
Because each tree was fitted to residuals computed from the current model, and it is only a good correction near that point. Taking the full step overshoots into territory the tree knows nothing about. Small steps, many of them, track the true function more carefully — the same argument as a learning rate in [[optimisers-sgd-to-adam]], and it produces the same trade: lower $\eta$ needs more rounds.
The table below shows both halves — training error and held-out error — because they say different things.
==============================================================================
LEARNING RATE AND NUMBER OF ROUNDS
==============================================================================
Training error always prefers the bigger step
==============================================================================
eta 5 trees 20 trees 100 trees
----------------------------------------
0.05 20.4023 10.4408 1.8296
0.1 15.7170 5.5012 0.8946
0.5 4.5781 1.0351 0.0302
1.0 3.9095 0.9189 0.0016
----------------------------------------
Held-out error does not
==============================================================================
eta 5 trees 20 trees 100 trees 400 trees
---------------------------------------------------
0.05 21.6237 13.0418 5.8020 4.8544
0.1 17.5169 8.8375 5.0782 4.8759
0.5 8.5494 6.8702 6.9667 6.9842
1.0 10.3849 10.8073 11.2455 11.2466
---------------------------------------------------
best held-out: eta=0.05, 400 trees -> 4.8544
eta=1.0 reaches its best early and then gets WORSE with more trees.
Small eta keeps improving for longer. Shrinkage and number of rounds
are one decision, not two -- which is why tuning them independently
wastes most of a grid ([[hyperparameter-search]]).
4. Changing the loss changes the algorithm
This is the payoff for phrasing it as gradients rather than residuals. Swap the loss and everything else stays put:
| loss | $-\partial L/\partial F$ | behaviour |
|---|---|---|
| squared error $\tfrac12(y-F)^2$ | $y - F$ | the plain residual; sensitive to outliers |
| absolute error $\lvert y-F \rvert$ | $\operatorname{sign}(y - F)$ | every row pulls with equal force — robust |
| Huber | residual, clipped | squared near zero, absolute in the tail |
| log loss (binary) | $y - p$, with $p = \sigma(F)$ | classification, on the log-odds scale |
Look at the absolute-error row. The gradient is $\pm 1$ regardless of how wrong the prediction is, so a wildly outlying point pulls no harder than a slightly wrong one. That is the entire robustness mechanism, and it falls straight out of the derivative.
==============================================================================
THE SAME MACHINERY, DIFFERENT LOSSES
==============================================================================
Pseudo-residuals for one badly outlying point
==============================================================================
recording rec03 corrupted: 4.8 -> 60.0
id y y - F squared-err grad abs-err grad
----------------------------------------------------------
rec00 0.2 -3.6000 -3.6000 -1.0
rec02 6.4 2.6000 2.6000 1.0
rec03 60.0 56.2000 56.2000 1.0
rec05 9.1 5.3000 5.3000 1.0
----------------------------------------------------------
The outlier's squared-error gradient is 56.2 -- it will
dominate the next tree entirely. Its absolute-error gradient is 1,
exactly the same pull as any other under-predicted row.
What that does to the fitted model
==============================================================================
12 rows is too few to show this: with 60 trees the model simply
memorises every row, outlier included, and nothing is damaged.
So this needs a real dataset and a held-out CLEAN test set.
240 training rows, 14 of them (6%) corrupted by +/-60
test set is 160 CLEAN rows the model never sees
loss test MAE
----------------------------
squared error 5.1901
huber 2.8239
absolute error 1.7613
----------------------------
(no corruption) 1.5127 <- the target to recover
Squared error is 3.4x the clean baseline: 14 corrupted rows
out of 240 dominate every tree, because a residual of 60
contributes a gradient of 60.
Absolute error gets within 0.25 of that
baseline, because those rows pull with force 1, exactly like
everything else. Huber sits between the two, by construction.
ordering as predicted (abs < huber < squared): True
Binary log loss: the gradient is y - p
==============================================================================
base rate = 0.5833 -> F_0 = log-odds = 0.3365
p = sigmoid(F_0) = 0.5833 (recovers the base rate: True)
id y p grad = y - p
------------------------------------
rec00 0 0.5833 -0.5833
rec01 0 0.5833 -0.5833
rec02 1 0.5833 0.4167
rec03 1 0.5833 0.4167
rec04 0 0.5833 -0.5833
rec05 1 0.5833 0.4167
------------------------------------
Positives pull up by 1-p, negatives pull down by p. Boosting then
proceeds identically -- fit a tree to that column, add eta * h.
The trees live on the LOG-ODDS scale; sigmoid is applied at the end.
5. Where XGBoost goes further
Everything above uses only the first derivative — it is gradient descent in function space. XGBoost's central change is to use the second derivative too, making it Newton's method in function space.
Expand the loss to second order around the current prediction, with $g_i = \partial L/\partial F$ and $h_i = \partial^2 L/\partial F^2$:
$$L \approx \sum_i \left[ g_i f(x_i) + \tfrac{1}{2} h_i f(x_i)^2 \right] + \Omega(f)$$
Minimising over a leaf containing rows $I$, with $L_2$ regularisation $\lambda$, gives the optimal leaf weight directly:
$$w^*_I = -\frac{\sum_{i \in I} g_i}{\sum_{i \in I} h_i + \lambda}$$
Compare with plain gradient boosting, where the leaf value is just the mean of the gradients in that leaf. The denominator is the difference: leaves whose rows have high curvature get smaller steps. For squared error $h_i = 1$, so $w^* = -\bar{g}\cdot n/(n+\lambda)$ and the two coincide apart from shrinkage — which is why the distinction is invisible on regression and matters for classification, where $h_i = p_i(1-p_i)$ varies a great deal across rows.
==============================================================================
FIRST-ORDER vs SECOND-ORDER LEAF VALUES
==============================================================================
For squared error the two agree (h = 1 everywhere)
==============================================================================
leaf with 6 rows
gradient-boosting leaf = mean(residual) = -0.208333
XGBoost leaf = -sum(g)/(sum(h)+lambda) = -0.208333
MATCH: True
For log loss they do NOT (h = p(1-p) varies)
==============================================================================
id p g = -(y-p) h = p(1-p)
----------------------------------------
rec00 0.0759 0.0759 0.0701
rec01 0.1419 0.1419 0.1217
rec02 0.9002 -0.0998 0.0898
rec03 0.8022 -0.1978 0.1587
rec04 0.0474 0.0474 0.0452
rec05 0.9427 -0.0573 0.0540
rec06 0.7109 -0.2891 0.2055
rec07 0.8699 -0.1301 0.1132
rec08 0.1091 0.1091 0.0972
rec09 0.9168 -0.0832 0.0763
rec10 0.3543 0.3543 0.2288
rec11 0.7503 -0.2497 0.1874
----------------------------------------
h ranges from 0.0452 to 0.2288 -- a 5x spread.
Confident rows (p near 0 or 1) have tiny curvature and barely move.
lambda=0.0: first-order leaf = +0.014959 second-order leaf = +0.166358 ratio 11.121
lambda=1.0: first-order leaf = +0.014959 second-order leaf = +0.058301 ratio 3.897
This is the step from this page to [[xgboost]]: same loop, same
trees, a better leaf value -- plus the regularisation term that
lambda is doing in that denominator.
6. What to take away
"Fit the residual" is a special case, not the definition. Gradient boosting fits the negative gradient; for squared error that happens to be the residual. Phrasing it as a gradient is what lets the same loop do robust regression, classification, ranking, and survival analysis.
The loss determines the behaviour. Squared error chases outliers because its gradient grows with the error; absolute error ignores them because its gradient is $\pm 1$. Neither is a special robustness feature — both are just derivatives.
Shrinkage and rounds are one decision. A large $\eta$ minimises training error fastest and generalises worst. Tune the pair together.
XGBoost is this plus a second derivative. Leaf value $-\sum g / (\sum h + \lambda)$ rather than the mean gradient. Identical for squared error, materially different for log loss.
Verified in this notebook
- Three hand-built rounds reproduce
GradientBoostingRegressorto machine precision, including the initial constant - Each stump's leaf values equal the mean of the residuals falling in that leaf
- $\eta = 1.0$ gets worse with more trees on held-out data while $\eta = 0.05$ keeps improving
- The first- and second-order leaf formulas agree exactly for squared error and diverge for log loss, where curvature spans a large range
Where this connects
[[adaboost]] is the same idea with exponential loss and row reweighting. [[xgboost]] adds the second-order step and explicit regularisation; [[lightgbm]] and [[catboost]] change how the trees are grown and how splits are found. [[optimisers-sgd-to-adam]] is the same learning-rate argument in parameter space rather than function space.
References
- Greedy Function Approximation: A Gradient Boosting Machine
- Stochastic Gradient Boosting