Gradient Boosting — How XGBoost and LightGBM Work

You are learning to throw darts. First throw: you miss the bullseye by 8cm to the right. A coach watches and says "next throw, aim 8cm to the left of wherever you aimed before." Second throw: miss by 3cm upward. Coach: "aim 3cm down from last time." Each throw corrects the accumulated error of all previous throws.

Gradient Boosting trains each new tree to hit where the previous ensemble missed. The final prediction is the sum of all trees — each one having corrected the previous collection's errors.

import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.metrics import mean_absolute_error

np.random.seed(42)
n = 1000

distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep])
y = delivery

# ── Gradient Boosting from scratch ────────────────────────────────────
class GradientBoostingScratch:
    """
    Gradient Boosting Regressor built from scratch.
    Every step is explicit — no hidden magic.
    """
    def __init__(self, n_estimators=50, learning_rate=0.1, max_depth=3):
        self.n_estimators   = n_estimators
        self.learning_rate  = learning_rate
        self.max_depth      = max_depth
        self.trees          = []
        self.initial_pred   = None

    def fit(self, X, y):
        # Step 1: Initial prediction = mean of target
        self.initial_pred = y.mean()
        current_pred = np.full(len(y), self.initial_pred)

        print(f"Initial prediction (mean): {self.initial_pred:.2f}")
        print(f"Initial MAE: {mean_absolute_error(y, current_pred):.4f}
")

        for i in range(self.n_estimators):
            # Step 2: Compute residuals — what the ensemble got wrong
            residuals = y - current_pred   # the "negative gradient"

            # Step 3: Train a shallow tree on the RESIDUALS (not y!)
            tree = DecisionTreeRegressor(max_depth=self.max_depth)
            tree.fit(X, residuals)
            self.trees.append(tree)

            # Step 4: Update predictions — small step in residual direction
            update = self.learning_rate * tree.predict(X)
            current_pred = current_pred + update

            if (i + 1) % 10 == 0:
                mae = mean_absolute_error(y, current_pred)
                print(f"  After tree {i+1:3d}: MAE = {mae:.4f}  "
                      f"mean residual = {np.abs(residuals).mean():.4f}")

        return self

    def predict(self, X):
        # Start with initial prediction
        pred = np.full(len(X), self.initial_pred)
        # Add each tree's contribution (scaled by learning rate)
        for tree in self.trees:
            pred += self.learning_rate * tree.predict(X)
        return pred


# Train from scratch
from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

print("=== Gradient Boosting From Scratch ===")
gb_scratch = GradientBoostingScratch(
    n_estimators=50, learning_rate=0.1, max_depth=3
)
gb_scratch.fit(X_tr, y_tr)
scratch_mae = mean_absolute_error(y_te, gb_scratch.predict(X_te))
print(f"
Final test MAE: {scratch_mae:.4f} min")

# Verify with sklearn
from sklearn.ensemble import GradientBoostingRegressor
gb_sk = GradientBoostingRegressor(
    n_estimators=50, learning_rate=0.1, max_depth=3, random_state=42
)
gb_sk.fit(X_tr, y_tr)
sk_mae = mean_absolute_error(y_te, gb_sk.predict(X_te))
print(f"sklearn GBR test MAE: {sk_mae:.4f} min  (should be similar)")

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 3000
distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value     = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep, value])
y = delivery

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── Learning rate × n_estimators trade-off ───────────────────────────
print(f"{'learning_rate':<16} {'n_estimators':<15} {'Train MAE':<12} {'Test MAE':<12} {'CV MAE'}")
print("─" * 70)

configs = [
    (0.3,   50),
    (0.1,  100),
    (0.1,  300),
    (0.05, 300),
    (0.05, 500),
    (0.01, 1000),
    (0.01, 2000),
]
for lr, n_est in configs:
    gb = GradientBoostingRegressor(
        learning_rate=lr, n_estimators=n_est,
        max_depth=4, subsample=0.8,
        random_state=42,
    )
    gb.fit(X_tr, y_tr)
    tr_mae = mean_absolute_error(y_tr, gb.predict(X_tr))
    te_mae = mean_absolute_error(y_te, gb.predict(X_te))
    cv     = -cross_val_score(gb, X_tr, y_tr, cv=3,
                               scoring='neg_mean_absolute_error').mean()
    flag   = ' ← overfit' if tr_mae < te_mae * 0.75 else ''
    print(f"  {lr:<14}  {n_est:<13}  {tr_mae:<12.4f} {te_mae:<12.4f} {cv:.4f}{flag}")

print("
Key observation:")
print("  Lower lr + more trees generally gives better test MAE")
print("  Very low lr (0.01) with enough trees usually wins but trains slowly")

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score, RandomizedSearchCV
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

# ── Effect of regularisation parameters ───────────────────────────────
base_params = dict(
    n_estimators=200, learning_rate=0.1,
    random_state=42,
)

print("Regularisation parameter effects (CV MAE):")
print(f"{'Config':<45} {'CV MAE'}")
print("─" * 55)

configs = [
    ('No regularisation (default)',
     dict(max_depth=3, subsample=1.0, min_samples_leaf=1, max_features=None)),
    ('Shallow trees (max_depth=2)',
     dict(max_depth=2, subsample=1.0, min_samples_leaf=1, max_features=None)),
    ('Deep trees (max_depth=6)',
     dict(max_depth=6, subsample=1.0, min_samples_leaf=1, max_features=None)),
    ('Stochastic GB (subsample=0.8)',
     dict(max_depth=3, subsample=0.8, min_samples_leaf=1, max_features=None)),
    ('Min leaf samples=10',
     dict(max_depth=3, subsample=1.0, min_samples_leaf=10, max_features=None)),
    ('Random features (max_features=sqrt)',
     dict(max_depth=3, subsample=1.0, min_samples_leaf=1, max_features='sqrt')),
    ('All regularisation combined',
     dict(max_depth=3, subsample=0.8, min_samples_leaf=5, max_features='sqrt')),
]

for name, reg_params in configs:
    model = GradientBoostingRegressor(**base_params, **reg_params)
    cv    = -cross_val_score(model, X_tr, y_tr, cv=5,
                              scoring='neg_mean_absolute_error').mean()
    print(f"  {name:<43} {cv:.4f}")

# ── RandomizedSearchCV — full hyperparameter search ───────────────────
param_dist = {
    'n_estimators':     [100, 200, 300, 500],
    'learning_rate':    [0.01, 0.05, 0.1, 0.2],
    'max_depth':        [2, 3, 4, 5],
    'subsample':        [0.6, 0.7, 0.8, 0.9, 1.0],
    'min_samples_leaf': [1, 5, 10, 20],
    'max_features':     ['sqrt', 'log2', None],
}

search = RandomizedSearchCV(
    GradientBoostingRegressor(random_state=42),
    param_dist,
    n_iter=30, cv=5,
    scoring='neg_mean_absolute_error',
    random_state=42, n_jobs=-1,
)
search.fit(X_tr, y_tr)
print(f"
Best params: {search.best_params_}")
print(f"Best CV MAE: {-search.best_score_:.4f} min")
best_test = mean_absolute_error(y_te, search.predict(X_te))
print(f"Test MAE:    {best_test:.4f} min")

import numpy as np
from sklearn.ensemble import GradientBoostingClassifier, GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import roc_auc_score, log_loss, mean_absolute_error

np.random.seed(42)
n = 3000
distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value     = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep, value])
y_reg = delivery
y_cls = (delivery > 45).astype(int)

X_tr, X_te, y_tr_r, y_te_r = train_test_split(X, y_reg, test_size=0.2, random_state=42)
_, _,  y_tr_c, y_te_c       = train_test_split(X, y_cls, test_size=0.2,
                                                 stratify=y_cls, random_state=42)

# ── Different loss functions for regression ────────────────────────────
print("Regression loss functions:")
for loss in ['squared_error', 'absolute_error', 'huber']:
    gb = GradientBoostingRegressor(
        loss=loss, n_estimators=200, learning_rate=0.1,
        max_depth=3, random_state=42,
    )
    gb.fit(X_tr, y_tr_r)
    mae = mean_absolute_error(y_te_r, gb.predict(X_te))
    print(f"  loss='{loss}': MAE={mae:.4f}")

# huber loss: squared_error for small residuals, absolute_error for large
# More robust to outliers than squared_error, faster to converge than absolute_error

# ── Classification — uses log-loss internally ─────────────────────────
gb_cls = GradientBoostingClassifier(
    n_estimators=200, learning_rate=0.1,
    max_depth=3, subsample=0.8,
    random_state=42,
)
gb_cls.fit(X_tr, y_tr_c)

y_proba = gb_cls.predict_proba(X_te)[:, 1]
print(f"
Classification (log-loss internally):")
print(f"  ROC-AUC:  {roc_auc_score(y_te_c, y_proba):.4f}")
print(f"  Log-loss: {log_loss(y_te_c, y_proba):.4f}")

# ── Staged predictions — watch the model improve over iterations ───────
print("
Gradient Boosting convergence (staged predictions):")
print(f"{'n_trees':<10} {'Train MAE':<12} {'Test MAE'}")
print("─" * 35)
gb_staged = GradientBoostingRegressor(
    n_estimators=500, learning_rate=0.1, max_depth=3, random_state=42
)
gb_staged.fit(X_tr, y_tr_r)

# staged_predict gives predictions at each stage of boosting
staged = list(gb_staged.staged_predict(X_te))
staged_tr = list(gb_staged.staged_predict(X_tr))

for n_trees in [1, 5, 10, 25, 50, 100, 200, 500]:
    tr_mae = mean_absolute_error(y_tr_r, staged_tr[n_trees - 1])
    te_mae = mean_absolute_error(y_te_r, staged[n_trees - 1])
    print(f"  {n_trees:<8}  {tr_mae:<12.4f} {te_mae:.4f}")

import numpy as np
import time
from sklearn.ensemble import GradientBoostingRegressor, HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 10_000
distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value     = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep, value])
y = delivery
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ── sklearn HistGradientBoostingRegressor — faster than GBR ──────────
# sklearn's own histogram-based GB (similar speedup to LightGBM)
# Available since sklearn 0.21, stable since 1.0
models = {
    'sklearn GBR':      GradientBoostingRegressor(
                            n_estimators=200, learning_rate=0.1,
                            max_depth=3, random_state=42),
    'sklearn HistGBR':  HistGradientBoostingRegressor(
                            max_iter=200, learning_rate=0.1,
                            max_depth=3, random_state=42),
}

# Try XGBoost and LightGBM if installed
try:
    import xgboost as xgb
    models['XGBoost'] = xgb.XGBRegressor(
        n_estimators=200, learning_rate=0.1, max_depth=3,
        subsample=0.8, colsample_bytree=0.8,
        random_state=42, verbosity=0,
    )
except ImportError:
    print("XGBoost not installed: pip install xgboost")

try:
    import lightgbm as lgb
    models['LightGBM'] = lgb.LGBMRegressor(
        n_estimators=200, learning_rate=0.1, max_depth=3,
        subsample=0.8, colsample_bytree=0.8,
        random_state=42, verbose=-1,
    )
except ImportError:
    print("LightGBM not installed: pip install lightgbm")

print(f"{'Model':<22} {'Train time':>12} {'Test MAE':>10}")
print("─" * 48)
for name, model in models.items():
    t0 = time.time()
    model.fit(X_tr, y_tr)
    train_time = time.time() - t0
    mae = mean_absolute_error(y_te, model.predict(X_te))
    print(f"  {name:<20}  {train_time:>10.3f}s  {mae:>10.4f}")

# ── HistGradientBoostingRegressor handles NaN natively ────────────────
import pandas as pd
X_with_nan = X_tr.copy().astype(float)
X_with_nan[np.random.choice(len(X_with_nan), 200, replace=False), 0] = np.nan

# HistGBR handles NaN — no imputation needed
hgbr = HistGradientBoostingRegressor(max_iter=100, random_state=42)
hgbr.fit(X_with_nan, y_tr)   # NaN handled internally
print(f"
HistGBR with NaN input — no error, MAE: "
      f"{mean_absolute_error(y_te, hgbr.predict(X_te)):.4f}")

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder
from sklearn.model_selection import cross_validate, StratifiedKFold, KFold
from sklearn.metrics import mean_absolute_error
import joblib, warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 8_000

restaurants = ['Pizza Hut','KFC','Dominos','Chipotle',"McDonald's",'Subway']
cities      = ['Seattle','New York','Denver','Austin','Boston','Chicago']
time_slots  = ['breakfast','lunch','evening','dinner']

df = pd.DataFrame({
    'distance_km':     np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15),
    'traffic_score':   np.random.randint(1, 11, n).astype(float),
    'restaurant_prep': np.abs(np.random.normal(15, 5, n)).clip(5, 35),
    'order_value':     np.abs(np.random.normal(350, 150, n)).clip(50, 1200),
    'restaurant':      np.random.choice(restaurants, n),
    'city':            np.random.choice(cities, n),
    'time_slot':       np.random.choice(time_slots, n),
    'is_weekend':      np.random.randint(0, 2, n).astype(float),
})
# Introduce missing values — real data always has them
df.loc[np.random.choice(n, 200, replace=False), 'restaurant_prep'] = np.nan
df.loc[np.random.choice(n, 100, replace=False), 'traffic_score']   = np.nan

y = (8.6 + 7.3*df['distance_km'] + 0.8*df['restaurant_prep'].fillna(15)
     + 1.5*df['traffic_score'].fillna(5) + np.random.normal(0, 4, n)).clip(10, 120)

from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(df, y, test_size=0.2, random_state=42)

NUM_COLS = ['distance_km', 'traffic_score', 'restaurant_prep',
            'order_value', 'is_weekend']
CAT_COLS = ['restaurant', 'city', 'time_slot']

# ── HistGBR handles NaN natively — only encode categoricals ──────────
preprocessor = ColumnTransformer([
    ('num', 'passthrough', NUM_COLS),   # pass through — HistGBR handles NaN
    ('cat', OrdinalEncoder(
        handle_unknown='use_encoded_value', unknown_value=-1
    ), CAT_COLS),
], remainder='drop')

pipeline = Pipeline([
    ('prep',  preprocessor),
    ('model', HistGradientBoostingRegressor(
        max_iter=300,
        learning_rate=0.05,
        max_depth=4,
        min_samples_leaf=10,
        l2_regularization=0.1,    # L2 regularisation (like Ridge penalty on leaf values)
        max_bins=255,              # histogram bins — higher = more accurate but slower
        early_stopping=True,       # stop if validation score stops improving
        validation_fraction=0.1,   # fraction of training data for early stopping
        n_iter_no_change=20,       # patience — stop after 20 rounds no improvement
        random_state=42,
    )),
])

# ── Cross-validation ──────────────────────────────────────────────────
cv_results = cross_validate(
    pipeline, df, y,
    cv=KFold(n_splits=5, shuffle=True, random_state=42),
    scoring='neg_mean_absolute_error',
    return_train_score=True,
)
print(f"Train MAE: {-cv_results['train_score'].mean():.4f} ± {cv_results['train_score'].std():.4f}")
print(f"Val MAE:   {-cv_results['test_score'].mean():.4f} ± {cv_results['test_score'].std():.4f}")

# ── Train final model and evaluate ────────────────────────────────────
pipeline.fit(X_tr, y_tr)
test_mae = mean_absolute_error(y_te, pipeline.predict(X_te))
print(f"Test MAE:  {test_mae:.4f} min")

# ── Feature importance ────────────────────────────────────────────────
model = pipeline.named_steps['model']
all_cols = NUM_COLS + CAT_COLS
importance = pd.Series(
    model.feature_importances_,
    index=all_cols,
).sort_values(ascending=False)

print("
Feature importance:")
for feat, imp in importance.items():
    bar = '█' * int(imp * 60)
    print(f"  {feat:<20}: {bar} {imp:.4f}")

# ── Save for production ───────────────────────────────────────────────
joblib.dump(pipeline, '/tmp/doordash_delivery_gbm.pkl')
print("
Model saved: /tmp/doordash_delivery_gbm.pkl")

# ── Score new orders ─────────────────────────────────────────────────
new_orders = pd.DataFrame([
    {'distance_km':5.2,'traffic_score':8,'restaurant_prep':22,
     'order_value':450,'restaurant':'KFC','city':'Seattle',
     'time_slot':'dinner','is_weekend':1},
    {'distance_km':1.5,'traffic_score':2,'restaurant_prep':10,
     'order_value':180,'restaurant':'Dominos','city':'New York',
     'time_slot':'lunch','is_weekend':0},
])
preds = pipeline.predict(new_orders)
for i, (_, row) in enumerate(new_orders.iterrows()):
    print(f"
Order {i+1}: {row['city']} | {row['restaurant']} | "
          f"{row['distance_km']}km | {row['time_slot']}")
    print(f"  Predicted delivery time: {preds[i]:.1f} min")

It is not just branding — "gradient" refers to a specific, precise quantity: the negative gradient of the loss function with respect to the current prediction. For MSE loss this negative gradient happens to equal the plain residual (y − ŷ), which is why Section 2 could get away with saying "each tree predicts the residual." But swap in log-loss for classification, and the tree is fitting the gradient of log-loss (a probability residual, y − sigmoid(ŷ)) — not a literal label difference. The "residual" framing is a special case of the gradient framing that only looks identical because MSE's gradient is unusually simple.

Lower learning rate paired with more trees does generally improve test performance, but it is not a free lunch — it is a straight trade of training compute for a small accuracy gain, with diminishing returns. The learning-rate sweep earlier in this module needed 2,000 trees at lr=0.01 to match what 50 trees did at lr=0.3 — a 40× increase in training cost for a modest MAE improvement. Past some point, additional trees stop helping regardless of how small the learning rate is, which is exactly why you tune both together with cross-validation and let early stopping — not intuition — decide where to stop.

It is usually the opposite. Random Forest averages many independent trees trained on bootstrap samples — that averaging cancels out noise by construction, whether you regularise it or not. Gradient boosting has no such built-in safety net: each new tree is explicitly optimising away whatever error remains, including error that is just noise in the training labels. Left unconstrained, later trees will happily fit that noise. This is exactly why gradient boosting needs the four regularisation levers from Section 4 (shallow trees, subsampling, min_samples_leaf, feature sampling) tuned deliberately — Random Forest gets comparable protection almost for free from its bagging mechanism.

It wins on most tabular benchmarks, which is why it dominates competitions — but "most" is not "always." On small or noisy-labelled datasets, an untuned Random Forest is often more robust out of the box, because bagging's variance-cancelling averaging degrades gracefully with poor hyperparameters, whereas an untuned gradient boosting model (high learning rate, no regularisation) can overfit badly and underperform a default Random Forest. Gradient boosting earns its accuracy advantage only when you invest the tuning time — an untuned comparison is not a fair one.

Random Forest trees are trained completely independently, so building 500 of them across 500 cores is trivial. Gradient boosting cannot do this: tree i+1 needs the residuals produced by the ensemble through tree i, so trees must be built one after another — there is a hard sequential dependency between boosting rounds. What XGBoost and LightGBM actually parallelise is the work *inside* each tree — histogram construction and split-finding across features run on multiple threads — but the 500 rounds themselves still happen in sequence. This is why gradient boosting training time scales roughly linearly with n_estimators almost regardless of core count, unlike Random Forest which genuinely speeds up with more parallel hardware.

Start with a single constant prediction — the mean of the target for regression. Compute the residual for every training point: actual value minus the current ensemble's prediction. Fit a new, shallow tree whose target is that residual, not the original label. Scale that tree's output by the learning rate and add it to the running ensemble prediction. Repeat for n_estimators rounds — each new tree is trained specifically to correct whatever error is left after every previous tree. The final prediction is the initial constant plus the learning rate times the sum of every tree's output.

For MSE loss, the residual (y − ŷ) is exactly the negative gradient of the loss with respect to the current prediction, so fitting a tree on the residual is mathematically the same as taking one step of gradient descent — except the "step" is an entire decision tree function instead of a fixed numeric update. This generalises beyond MSE: for any differentiable loss (log-loss for classification, a custom ranking loss), you compute the negative gradient — the pseudo-residual — at every point and fit the next tree on that. XGBoost goes one step further and also uses the second derivative (the Hessian) for a more accurate Newton-style update.

Learning rate scales down each tree's contribution; a smaller value means gentler, more conservative corrections that require more trees to converge but generalise better because no single tree can dominate the ensemble. A reliable rule of thumb is that halving the learning rate while doubling n_estimators tends to improve test performance. In practice you rarely hand-tune n_estimators directly — set it very high (thousands) with a small learning rate (0.01–0.05), then use early stopping against a validation set so training automatically halts at the point where validation error stops improving, regardless of the learning rate chosen.

Random Forest trains trees independently on bootstrap samples and averages them — that is inherently robust, needs very little hyperparameter tuning to get a strong baseline, trains in parallel across cores, and tends to tolerate noisy labels better because errors get averaged away rather than explicitly chased. Gradient boosting, especially with a high learning rate or too many trees, can overfit noisy data because later trees are built specifically to reduce whatever training residual remains, noise included. I would reach for Random Forest when I need a fast, low-maintenance baseline or I'm working with a small/noisy dataset, and reach for gradient boosting when I have the time budget to tune it properly and want the best possible accuracy on structured tabular data.

First symptom: a large gap between training and validation/test error — near-perfect training score with much worse held-out score is the classic signature. Fix it using the four regularisation levers together, not just one: lower the learning rate and correspondingly raise n_estimators with early stopping enabled; keep max_depth shallow (3–5, since gradient boosting trees are meant to be weak learners, unlike Random Forest's often-unconstrained trees); add subsample below 1.0 so each tree only sees a random fraction of rows (stochastic gradient boosting); and raise min_samples_leaf so trees can't carve out single-sample leaves. Tune all of this with cross-validation, never by reading training score alone.