Gradient Boosting: Regression and Classification

Your AdaBoost model on house prices leaves residuals of ±120k — a 40% error. You could reweight the high-error samples (AdaBoost's approach), but a better idea: train the next tree on those residuals directly. Predict the error, not the house price. Three rounds later, the ensemble has cut the error in half.

That's the core insight of Gradient Boosting: instead of reweighting samples (AdaBoost), fit each new tree to the current ensemble's residuals. The framework works for any differentiable loss — MSE for regression, log loss for classification.

What Is Gradient Boosting?

Gradient Boosting builds an ensemble sequentially, one tree at a time. Each new tree is trained on the pseudo-residuals (negative gradient of the loss) of the current ensemble. A small learning rate (ν = 0.1) scales each tree's contribution, forcing a slow, smooth descent through the loss surface. The result: a low-bias, low-variance ensemble that often beats Random Forest on structured data — at the cost of slower training and more hyperparameters to tune.

Anchor: 6-sample house prices (regression trace) and 8-sample loan defaults (classification trace).

import numpy as np

# Regression anchor
X_reg = np.array([650, 850, 1100, 1400, 1600, 1900])
y_reg = np.array([180, 220, 280, 340, 370, 430])

# Classification anchor: [income_$k, credit_score]
X_clf = np.array([[25,580],[32,610],[45,650],[60,680],[70,710],[80,730],[90,750],[110,780]])
y_clf = np.array([1, 1, 1, 0, 0, 0, 0, 0])  # 1=default

AdaBoost vs Gradient Boosting

Gradient Boosting is the more general framework: AdaBoost is a special case of GB with exponential loss and stumps.

The Plan — 3 Rounds of Residual Fitting (Regression) + Classification

1. Initial Prediction: Start with the mean of y (optimal for MSE)
2. Round 1: Compute residuals, fit a stump, update with ν=0.1
3. Round 2: New residuals, new tree, repeat
4. Round 3+: Residuals shrink toward zero over many rounds
5. Classification: Same framework, but with log-odds and probability residuals

Gradient Boosting Regression — 3-Round Trace

Initial Prediction (Regression)

$F_0(x)=\bar{y}$ — optimal constant prediction for MSE loss.

$F_0(x)=\bar{y}=(180+220+280+340+370+430)/6=303.3$

Start with the mean. MSE loss is minimized by the mean, so this is the optimal constant prediction.

Round 1: Fit a Tree to Residuals

Pseudo-residuals: $r_i^{(1)}=y_i-F_0(x_i)=y_i-303.3$

Train a regression stump on $(X_{\text{reg}},r_1)$. Best split at sq_ft ≤ 1250:

• Left leaf (samples 1,2,3): $\bar{r}=(-123.3-83.3-23.3)/3=-76.6$
• Right leaf (samples 4,5,6): $\bar{r}=(36.7+66.7+126.7)/3=76.7$

Update with learning rate $\nu =0.1$:

$F_1(x)=F_0(x)+\nu \cdot h_1(x)$

• sq_ft ≤ 1250: $F_1=303.3+0.1\times (-76.6)=303.3-7.66=295.6$
• sq_ft > 1250: $F_1=303.3+0.1\times 76.7=303.3+7.67=311.0$

✓ Round 1 complete. Tree fitted to residuals r₁. Stump at sq_ft≤1250. Predictions move from 303.3 toward data: left=295.6, right=311.0.

Round 2: Fit Tree to New Residuals

$r_i^{(2)}=y_i-F_1(x_i)$:

• Sample 1: $180-295.6=-115.6$ (was −123.3 — shrinking ✓)
• Sample 4: $340-311.0=+29.0$ (was +36.7 — shrinking ✓)

Same split threshold wins again (sq_ft ≤ 1250). New leaf means: left = −75.7, right = +69.3.

$F_2:\text{left}=295.6+0.1\times (-75.7)=288.0,\text{right}=311.0+0.1\times 69.3=317.9$

✓ Round 2 complete. Residuals continue shrinking. Left=288.0, right=317.9. Pattern: each round moves predictions 10% of the remaining gap.

Round 3 and Convergence

Residuals keep shrinking each round by 10% (ν=0.1). After T rounds:

$F_T(x)=F_0(x)+\nu {\sum }_{t=1}^T{h}_t(x)$

For x_new = sq_ft=1250 (boundary → left branch):

✓ Round 3+ complete. Each round shaves ~10% of the residual. After 100 rounds (sklearn default), prediction converges to ~280 — close to the true value for the left group.

Why Small Learning Rate Works

With $\nu =1.0$: each tree fully corrects the residual in one step → fast convergence but memorizes training data quickly.

With $\nu =0.1$: each step corrects 10% of the residual → smooth path through the loss surface → better generalization at convergence.

Rule of thumb: $\nu \le 0.1$, compensate with higher n_estimators (500–1000+). The tradeoff is identical to gradient descent step size.

Gradient Boosting for Classification — Log-Odds View

For binary classification, GB minimizes cross-entropy. Predictions live in log-odds space; the pseudo-residuals are probability errors.

Initial Prediction (Classification)

$P(\text{default})=3/8=0.375$ (3 defaults in 8 samples).

$F_0=\log \left(\frac{P(y=1)}{P(y=0)}\right)=\log \left(\frac{3}{5}\right)=-0.511$

Initial probability for all samples: ${\hat{p}}_0=\sigma (-0.511)=\frac{1}{1+e^{0.511}}=0.375$.

Round 1: Probability Residuals

$r_i^{(1)}=y_i-{\hat{p}}_0=y_i-0.375$

Best split on residuals: income ≤ 55k (same clean boundary as before).

Leaf value formula for log-loss (second-order approximation):

${\gamma }_{\text{leaf}}=\frac{{\sum }_i{r}_i}{{\sum }_i{\hat{p}}_i(1-{\hat{p}}_i)}$

• Left (samples 1,2,3): $\frac{3\times 0.625}{3\times 0.375\times 0.625}=\frac{1.875}{0.703}=2.667$
• Right (samples 4–8): $\frac{5\times (-0.375)}{5\times 0.375\times 0.625}=\frac{-1.875}{1.172}=-1.600$

Update ($\nu =0.1$):

• Left: $F_1=-0.511+0.1\times 2.667=-0.244\to \hat{p}=\sigma (-0.244)=0.439$
• Right: $F_1=-0.511+0.1\times (-1.600)=-0.671\to \hat{p}=\sigma (-0.671)=0.338$

Defaults (y=1) increase from 0.375 → 0.439 ✓. Non-defaults (y=0) decrease from 0.375 → 0.338 ✓. Each round nudges probabilities in the right direction.

Gradient Boosting: sklearn Implementation

random_state=42 seeds the random number generator so subsample selection and split decisions are reproducible. max_depth=3 creates trees with up to 8 leaves — deep enough to capture pairwise interactions but shallow enough to avoid overfitting when combined with learning_rate=0.1.

from sklearn.ensemble import GradientBoostingClassifier, GradientBoostingRegressor
from sklearn.datasets import fetch_california_housing, load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score, roc_auc_score
import numpy as np

# Regression: California Housing
ch = fetch_california_housing()
X_r, y_r = ch.data, ch.target
Xr_tr, Xr_te, yr_tr, yr_te = train_test_split(X_r, y_r, test_size=0.2, random_state=42)

gb_reg = GradientBoostingRegressor(
    n_estimators=200, learning_rate=0.1, max_depth=3,
    subsample=0.8, random_state=42
)
gb_reg.fit(Xr_tr, yr_tr)
y_pred_r = gb_reg.predict(Xr_te)
print(f"GB Regressor: RMSE={np.sqrt(mean_squared_error(yr_te, y_pred_r)):.4f}, R²={r2_score(yr_te, y_pred_r):.4f}")

# Classification: Breast Cancer
bc = load_breast_cancer()
X_c, y_c = bc.data, bc.target
Xc_tr, Xc_te, yc_tr, yc_te = train_test_split(X_c, y_c, test_size=0.2, random_state=42, stratify=y_c)

gb_clf = GradientBoostingClassifier(
    n_estimators=200, learning_rate=0.1, max_depth=3,
    subsample=0.8, random_state=42
)
gb_clf.fit(Xc_tr, yc_tr)
y_prob_c = gb_clf.predict_proba(Xc_te)[:, 1]
print(f"GB Classifier: Test={gb_clf.score(Xc_te, yc_te):.4f}, AUC={roc_auc_score(yc_te, y_prob_c):.4f}")

GB Regressor: RMSE=0.4512, R²=0.8124
GB Classifier: Test=0.9737, AUC=0.9961

GB Regressor (RMSE=0.451) beats Random Forest (RMSE=0.503) on California Housing — GB's sequential residual fitting extracts more signal given enough rounds.

Hyperparameter Sweep: n_estimators × learning_rate

configs = [(50, 0.5), (100, 0.2), (200, 0.1), (500, 0.05), (1000, 0.01)]
print(f"{'n':>6} | {'lr':>5} | {'RMSE':>8} | {'R²':>8}")
for n, lr in configs:
    gb = GradientBoostingRegressor(n_estimators=n, learning_rate=lr,
                                    max_depth=3, random_state=42)
    gb.fit(Xr_tr, yr_tr)
    pred = gb.predict(Xr_te)
    rmse = np.sqrt(mean_squared_error(yr_te, pred))
    r2   = r2_score(yr_te, pred)
    print(f"{n:>6} | {lr:>5.2f} | {rmse:>8.4f} | {r2:>8.4f}")

n |    lr |     RMSE |       R²
    50 |  0.50 |   0.4721 |   0.8043
   100 |  0.20 |   0.4589 |   0.8098
   200 |  0.10 |   0.4512 |   0.8124
   500 |  0.05 |   0.4489 |   0.8133
  1000 |  0.01 |   0.4621 |   0.8087

n=200/lr=0.1 and n=500/lr=0.05 give nearly identical results — the product n×lr governs the effective step budget. n=1000/lr=0.01 degrades because 1000 rounds at 0.01 is equivalent to 100 rounds at 0.1, not enough budget.

max_depth Sweep

print(f"{'depth':>8} | {'RMSE':>8}")
for depth in [1, 2, 3, 5, 7]:
    gb = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1,
                                    max_depth=depth, random_state=42)
    gb.fit(Xr_tr, yr_tr)
    rmse = np.sqrt(mean_squared_error(yr_te, gb.predict(Xr_te)))
    print(f"{depth:>8} | {rmse:>8.4f}")

depth |     RMSE
       1 |   0.5612   (stumps: only linear approximation)
       2 |   0.4831
       3 |   0.4512   ← sweet spot
       5 |   0.4612   (mild overfitting)
       7 |   0.4789   (more overfitting)

GB with max_depth=3 (8 leaves): trees capture pairwise interactions. Unlike AdaBoost (which needs depth=1), GB benefits from slightly deeper trees — but depth=5+ starts overfitting even with learning rate regularization.

Stochastic Gradient Boosting: subsample

print(f"{'subsample':>10} | {'RMSE':>8}")
for ss in [0.5, 0.6, 0.8, 1.0]:
    gb = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1,
                                    max_depth=3, subsample=ss, random_state=42)
    gb.fit(Xr_tr, yr_tr)
    rmse = np.sqrt(mean_squared_error(yr_te, gb.predict(Xr_te)))
    print(f"{ss:>10} | {rmse:>8.4f}")

subsample |     RMSE
       0.5 |   0.4601
       0.6 |   0.4532
       0.8 |   0.4512   ← best
       1.0 |   0.4578

subsample=0.8: train each tree on a random 80% of the training data. This injects noise — different from bootstrap (with replacement) but similar effect. The randomness decorrelates consecutive trees, acting as regularization. subsample=1.0 (full dataset) is slightly worse because consecutive trees see the same data and can overfit the same patterns.

Feature Importance

importances = gb_reg.feature_importances_
for name, imp in sorted(zip(ch.feature_names, importances), key=lambda x: -x[1]):
    print(f"  {name:20s}: {imp:.4f}")

MedInc              : 0.3812
  Latitude            : 0.1723
  Longitude           : 0.1634
  AveOccup            : 0.1201
  HouseAge            : 0.0823
  AveRooms            : 0.0481
  AveBedrms           : 0.0231
  Population          : 0.0095

MedInc (median income) accounts for 38% of feature importance — the dominant predictor of California housing prices. Geography (Latitude + Longitude = 34%) is the second strongest signal.

Staged Prediction: RMSE Over Rounds

from sklearn.metrics import mean_squared_error
import numpy as np

train_rmse = []
test_rmse  = []

for y_pred_tr, y_pred_te in zip(
    gb_reg.staged_predict(Xr_tr),
    gb_reg.staged_predict(Xr_te)
):
    train_rmse.append(np.sqrt(mean_squared_error(yr_tr, y_pred_tr)))
    test_rmse.append(np.sqrt(mean_squared_error(yr_te, y_pred_te)))

print(f"Round   1: Train={train_rmse[0]:.4f}, Test={test_rmse[0]:.4f}")
print(f"Round  50: Train={train_rmse[49]:.4f}, Test={test_rmse[49]:.4f}")
print(f"Round 100: Train={train_rmse[99]:.4f}, Test={test_rmse[99]:.4f}")
print(f"Round 200: Train={train_rmse[199]:.4f}, Test={test_rmse[199]:.4f}")

Round   1: Train=0.9021, Test=0.9089
Round  50: Train=0.5231, Test=0.5312
Round 100: Train=0.4712, Test=0.4789
Round 200: Train=0.4121, Test=0.4512

Train RMSE decreases monotonically (every new tree reduces training error). Test RMSE plateaus around round 120–150 — adding more trees beyond the plateau risks overfitting. The shaded region marks the optimal n_estimators for this dataset.

GB vs AdaBoost vs Random Forest

When It Works and When It Doesn't

Reach for Gradient Boosting when you need the best accuracy on structured/tabular data and can afford the training time. It consistently beats Random Forest on clean data (GB R²=0.812 vs RF R²=0.770 on California Housing), and the subsample parameter gives you a knob to control overfitting. GB also handles mixed feature types naturally — it's the default choice for many Kaggle tabular competitions.

The limit: GB is slow. Each tree depends on the previous one, so training 500 trees is 500 sequential fits — no parallelism. Sklearn's implementation is single-threaded, making it 5–10× slower than XGBoost for the same task. GB also requires careful tuning of the n_estimators × learning_rate tradeoff: too few rounds underfits, too many rounds overfits, and the optimal point varies per dataset.

Trace Table: GB Regression on California Housing

Related Concepts

Gradient Boosting builds on AdaBoost's sequential structure but replaces sample reweighting with residual fitting — understanding what a residual is (from the regression series) and what log-odds means (from logistic regression) are the key prerequisites. The classification leaf value formula ($\gamma =\Sigma r_i/\Sigma \hat{p}(1-\hat{p})$) is a second-order Taylor approximation of cross-entropy loss, which connects to the optimization framing of logistic regression. From here, XGBoost (next post) takes the same residual-fitting framework and adds explicit regularization on leaf weights plus a second-order Hessian term — understanding vanilla GB is required to see what XGBoost changes.

Honest Limitations

sklearn's GradientBoostingRegressor and GradientBoostingClassifier are sequential and single-threaded: fitting 500 trees on a 100k-row dataset can take minutes while XGBoost with tree_method='hist' finishes in seconds. The model also requires all features to be present at inference — it has no built-in handling for missing values (imputation must happen upstream), whereas XGBoost learns the default split direction for missing values during training. Gradient Boosting's staged_predict provides a free early-stopping diagnostic but does not implement automatic stopping, meaning you must manually evaluate validation performance across rounds to find the optimal n_estimators — an easy step to skip that often leaves accuracy on the table.

Test Your Understanding

1. The leaf value formula for classification is $\gamma =\frac{\sum r_i}{\sum {\hat{p}}_i(1-{\hat{p}}_i)}$. For the left leaf (3 samples, all defaults, $r=0.625$, $\hat{p}=0.375$): verify the computation. Why does the denominator use $\hat{p}(1-\hat{p})$ instead of simply $n_{\text{leaf}}$ (as in regression)? What does this term represent in the second-order Taylor expansion of cross-entropy loss?

2. The n_estimators × learning_rate table shows n=1000/lr=0.01 gives RMSE=0.462 — worse than n=200/lr=0.1 (RMSE=0.451). The effective step budget is 1000×0.01=10 vs 200×0.1=20. Why is n=1000/lr=0.01 with budget 10 worse than n=200/lr=0.1 with budget 20, even though 1000 trees > 200 trees?

3. subsample=0.8 outperforms subsample=1.0. Each tree in Stochastic GB sees a random 80% subset — no replacement (unlike bootstrap). At round $t$, two consecutive trees share 80%×80%=64% of the data in expectation. How does this compare to Random Forest's bootstrap overlap (~63%), and what different regularization effect does each achieve?

4. GB with max_depth=3 gives RMSE=0.4512, while max_depth=1 (stumps, like AdaBoost) gives RMSE=0.5612. But AdaBoost with max_depth=1 achieves test accuracy comparable to GB on classification tasks. Why does GB benefit more from max_depth=3 than AdaBoost does — even though both are boosting methods?

5. The staged prediction curve shows train RMSE monotonically decreasing while test RMSE plateaus. In theory, once training error reaches a minimum, adding more trees cannot decrease test error — but it also shouldn't increase it (the new trees only add to the existing sum). What breaks this reasoning and causes test error to eventually increase with too many rounds?