LEARN · GRADIENT BOOSTING WITH XGBOOST
When hyperparameter tuning stops paying
Once learning rate, depth, regularization, subsampling, and early stopping are sensible, another search rarely creates a step-change. The next gains usually come from changing one of three things:
-
The information the trees receive
-
The errors the booster optimizes
-
The behavior the model is permitted to learn
The code below targets XGBoost 3.3.0, the current stable release as of June 17, 2026. It uses the current scikit-learn interface and hist tree method.
We will use synthetic property listings: a neutral tabular problem with nonlinear interactions, uneven neighborhood activity, heteroscedastic noise, and a few exceptional sales. Run the snippets in order.
import numpy as np import pandas as pd import xgboost as xgb from sklearn.metrics import mean_absolute_error, root_mean_squared_error from sklearn.model_selection import train_test_split SEED = 42 REFERENCE_YEAR = 2026 rng = np.random.default_rng(SEED) n = 14_000 zone_probabilities = rng.dirichlet(np.linspace(0.3, 2.0, 60)) zone_id = rng.choice(60, size=n, p=zone_probabilities) living_area = np.clip( rng.lognormal(np.log(1_350), 0.38, n), 450, 4_200, ) bedrooms = np.clip( np.rint(living_area / 520 + rng.normal(0, 0.8, n)), 1, 7, ).astype(int) bathrooms = np.clip( np.round((0.55 * bedrooms + rng.normal(0.8, 0.55, n)) * 2) / 2, 1, 5, ) X = pd.DataFrame( { "living_area": living_area, "bedrooms": bedrooms, "bathrooms": bathrooms, "year_built": rng.integers(1940, 2027, n), "distance_to_center": np.clip( rng.gamma(2.2, 4.5, n), 0.2, 45, ), "lot_size": np.clip( living_area * rng.uniform(1.1, 4.5, n) + rng.normal(0, 600, n), 500, 18_000, ), "garage_spaces": np.clip(rng.poisson(1.2, n), 0, 3), "pool": rng.binomial(1, 0.12, n), "solar": rng.binomial(1, 0.16, n), "renovated": rng.binomial(1, 0.28, n), "listing_month": rng.integers(1, 13, n), "zone_id": zone_id, } ) X["school_score"] = np.clip( 10 - 0.11 * X["distance_to_center"] + rng.normal(0, 1.1, n), 1, 10, ) X["transit_score"] = np.clip( 95 - 2 * X["distance_to_center"] + rng.normal(0, 12, n), 0, 100, ) age = REFERENCE_YEAR - X["year_built"] area_per_room = living_area / (bedrooms + bathrooms) amenity_count = ( X["pool"] + X["solar"] + X["renovated"] + (X["garage_spaces"] > 0).astype(int) ) zone_rarity = -np.log(zone_probabilities[zone_id] + 1e-8) seasonality = 12_000 * np.sin( 2 * np.pi * (X["listing_month"] - 3) / 12 ) expected_price = ( 85_000 + 205 * living_area + 36 * X["lot_size"] + 14_500 * X["school_score"] - 3_900 * X["distance_to_center"] + 420 * X["transit_score"] + 17_000 * X["garage_spaces"] + 28_000 * X["pool"] + 18_000 * X["solar"] + 13_000 * X["renovated"] - 950 * age + 90 * area_per_room + 14_000 * amenity_count + 24_000 * zone_rarity + seasonality + 42 * living_area * X["pool"] + 230 * age * X["renovated"] ) noise_scale = 18_000 + 0.035 * expected_price y = expected_price.to_numpy() + rng.normal( 0, noise_scale.to_numpy(), ) exceptional_sale = rng.random(n) < 0.015 y[exceptional_sale] += rng.normal( 0, 220_000, exceptional_sale.sum(), ) X_train, X_temp, y_train, y_temp = train_test_split( X, y, test_size=0.30, random_state=SEED, ) X_valid, X_test, y_valid, y_test = train_test_split( X_temp, y_temp, test_size=0.50, random_state=SEED, )
Freeze a serious baseline
Do not change the split while testing ideas. A favorable resample can look like model progress.
common_params = { "n_estimators": 1_200, "learning_rate": 0.035, "max_depth": 4, "min_child_weight": 8, "subsample": 0.85, "colsample_bytree": 0.85, "reg_lambda": 10.0, "reg_alpha": 0.05, "tree_method": "hist", "eval_metric": "rmse", "early_stopping_rounds": 70, "random_state": SEED, "n_jobs": -1, } baseline = xgb.XGBRegressor( objective="reg:squarederror", **common_params, ) baseline.fit( X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=False, ) baseline_prediction = baseline.predict(X_test) print( f"RMSE: " f"{root_mean_squared_error(y_test, baseline_prediction):,.0f}" ) print( f"MAE: " f"{mean_absolute_error(y_test, baseline_prediction):,.0f}" )
The seeded run produced 57,918 RMSE and 41,384 MAE. That is the number every later experiment must beat—or deliberately trade away for better behavior elsewhere.
Track both metrics. RMSE reacts strongly to exceptional sales, while MAE describes the typical absolute miss more directly. Also record best_iteration: when an experiment improves only because it trains many more trees, compare latency and model size before declaring victory.
For a serious project, repeat the ablation across several fixed folds, but never change the protocol midway through the comparison.
Engineer operations a shallow tree must approximate
Trees discover threshold interactions, but they do not directly divide two columns, calculate a periodic transform, or know how rare a group is. A depth-limited model must spend several splits approximating those operations.
High-value candidates include:
-
Ratios with a domain interpretation
-
Counts assembled from sparse indicators
-
Cyclical encodings for month, hour, or weekday
-
Group frequencies fitted on training data only
-
Gated interactions such as age × renovation status
def add_features( reference: pd.DataFrame, frame: pd.DataFrame, ) -> pd.DataFrame: zone_counts = reference["zone_id"].value_counts() result = frame.copy() result["property_age"] = ( REFERENCE_YEAR - result["year_built"] ) result["area_per_room"] = ( result["living_area"] / (result["bedrooms"] + result["bathrooms"]) ) result["bathroom_density"] = ( result["bathrooms"] / result["bedrooms"] ) result["amenity_count"] = ( result[["pool", "solar", "renovated"]].sum(axis=1) + (result["garage_spaces"] > 0).astype(int) ) angle = 2 * np.pi * (result["listing_month"] - 1) / 12 result["month_sin"] = np.sin(angle) result["month_cos"] = np.cos(angle) result["zone_listing_count"] = ( result["zone_id"] .map(zone_counts) .fillna(0) .astype(float) ) return result.drop(columns="year_built") Xe_train = add_features(X_train, X_train) Xe_valid = add_features(X_train, X_valid) Xe_test = add_features(X_train, X_test) engineered = xgb.XGBRegressor( objective="reg:squarederror", **common_params, ) engineered.fit( Xe_train, y_train, eval_set=[(Xe_valid, y_valid)], verbose=False, ) engineered_prediction = engineered.predict(Xe_test)
RMSE fell to 57,584, a modest 0.58% improvement. Small gains still matter when they come from stable information rather than search noise.
Fit frequency features on the training partition. Using validation or test rows would give the model knowledge unavailable at production time.
Ablate feature families separately. Ratios, counts, seasonality, and group frequency answer different hypotheses; adding them all at once hides which one helped. Check missing-value behavior too: an unseen zone maps to zero here, which is explicit and reproducible, but another application may need an unknown_zone indicator.
Constrain behavior you must defend
Monotonic constraints encode directional rules, while interaction constraints control which features may appear together along a decision path. Current XGBoost accepts feature-name dictionaries and nested feature-name lists.
With hist, monotonic rules can eliminate candidate splits; increasing max_bin gives the learner more candidates.
constrained = xgb.XGBRegressor( objective="reg:squarederror", monotone_constraints={ "living_area": 1, "school_score": 1, "distance_to_center": -1, }, interaction_constraints=[ [ "living_area", "bedrooms", "bathrooms", "lot_size", "area_per_room", "bathroom_density", "pool", ], [ "property_age", "renovated", "amenity_count", ], [ "distance_to_center", "school_score", "transit_score", "zone_id", "zone_listing_count", ], [ "garage_spaces", "pool", "solar", "amenity_count", ], [ "listing_month", "month_sin", "month_cos", ], ], max_bin=512, **common_params, ) constrained.fit( Xe_train, y_train, eval_set=[(Xe_valid, y_valid)], verbose=False, )
The constrained model reached 56,986 RMSE, another 1.04% improvement over the engineered model. Treat that accuracy gain as a bonus. The primary value is a model whose local direction and permitted interactions are easier to justify.
Do not encode a rule merely because it sounds plausible. More floor area may not be monotonic when the dataset mixes apartments, warehouses, and land parcels.
After fitting, probe each constrained feature over a realistic grid while holding the others fixed. The constraint guarantees direction, but it does not guarantee that the resulting curve is economically sensible or well supported by data.
Use sample weights for cost-sensitive errors
Suppose errors on the top quarter of listings cost three times as much. Sample weights express that business preference without inventing a new loss.
high_value_cutoff = np.quantile(y_train, 0.75) train_weights = np.where( y_train >= high_value_cutoff, 3.0, 1.0, ) valid_weights = np.where( y_valid >= high_value_cutoff, 3.0, 1.0, ) weighted_model = xgb.XGBRegressor( objective="reg:squarederror", **common_params, ) weighted_model.fit( Xe_train, y_train, sample_weight=train_weights, eval_set=[(Xe_valid, y_valid)], sample_weight_eval_set=[valid_weights], verbose=False, ) weighted_prediction = weighted_model.predict(Xe_test) high_value_test = y_test >= high_value_cutoff print( mean_absolute_error( y_test[high_value_test], engineered_prediction[high_value_test], ) ) print( mean_absolute_error( y_test[high_value_test], weighted_prediction[high_value_test], ) )
Overall RMSE worsened from 57,584 to 58,127, but high-value MAE improved from 52,792 to 50,172—a 4.96% gain in the segment that mattered.
That is not a failed model; it is a different optimization target.
Derive the cutoff from training labels only, then freeze it. Recomputing the threshold on each evaluation set changes the business segment being measured. Also report unweighted metrics beside weighted ones so stakeholders can see where performance was purchased and where it was surrendered.
Write a custom objective with explicit derivatives
A custom objective for XGBRegressor returns one gradient and Hessian per row. XGBoost expects objectives that are smooth, twice differentiable, row-additive, and preferably convex; negative Hessians are clipped. Callable metrics in the scikit-learn interface are costs to minimize.
For residual r = prediction − target, use log-cosh:
-
Loss = δ² × log(cosh(r / δ))
-
Gradient = δ × tanh(r / δ)
-
Hessian = 1 − tanh²(r / δ)
DELTA = 50_000.0
def log_cosh_objective(
y_true: np.ndarray,
y_pred: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
scaled_residual = (y_pred - y_true) / DELTA
tanh_residual = np.tanh(scaled_residual)
gradient = DELTA * tanh_residual
hessian = 1.0 - tanh_residual**2
return gradient, np.maximum(hessian, 1e-6)
def mean_log_cosh(
y_true: np.ndarray,
y_pred: np.ndarray,
) -> float:
scaled_residual = (y_pred - y_true) / DELTA
stable_log_cosh = (
np.logaddexp(
scaled_residual,
-scaled_residual,
)
- np.log(2.0)
)
return float(
np.mean(DELTA**2 * stable_log_cosh)
)
custom_params = {
key: value
for key, value in common_params.items()
if key != "eval_metric"
}
custom_model = xgb.XGBRegressor(
objective=log_cosh_objective,
eval_metric=mean_log_cosh,
base_score=float(np.median(y_train)),
**custom_params,
)
custom_model.fit(
Xe_train,
y_train,
eval_set=[(Xe_valid, y_valid)],
verbose=False,
)
The custom loss reached 57,454 RMSE and 40,743 MAE: only a 0.23% RMSE improvement over ordinary squared error on engineered features.
The lesson is not “custom losses win.” It is “custom losses are worth maintaining only when their error shape matches the decision cost.”
Because this objective uses an identity link, raw predictions already live on the target scale. For non-identity links, XGBoost cannot infer the transformation for a custom objective; you must apply it consistently in both the objective and metric.
Before a long training run, compare the analytic gradient with a finite-difference approximation on a handful of rows. Most custom-loss failures are sign, scaling, or argument-order mistakes.
Setting base_score near the target median also gives the first trees a sensible starting point. XGBoost can estimate intercepts for selected built-in objectives, but it cannot infer the optimum for an arbitrary callable.
Try native pseudo-Huber before owning custom code
XGBoost includes reg:pseudohubererror, with huber_slope controlling δ. The current default is 1.0, which is usually far too small for targets measured in hundreds of thousands.
pseudo_huber = xgb.XGBRegressor( objective="reg:pseudohubererror", eval_metric="mphe", huber_slope=25_000.0, base_score=float(np.median(y_train)), **custom_params, ) pseudo_huber.fit( Xe_train, y_train, eval_set=[(Xe_valid, y_valid)], verbose=False, )
It produced 57,797 RMSE and 41,024 MAE—worse than squared error here.
Robust loss is not an automatic upgrade simply because outliers exist. Keep it when repeated validation or operational costs support it.
Cherry on the cake: uncertainty bands for one extra fit
The quantile objective accepts an array of quantiles. One booster can therefore learn the 5th, 50th, and 95th conditional quantiles and return three columns per row.
XGBoost recommends hist rather than exact for this objective and warns that quantile crossing can occur.
quantile_train = xgb.QuantileDMatrix( Xe_train, label=y_train, ) quantile_valid = xgb.QuantileDMatrix( Xe_valid, label=y_valid, ref=quantile_train, ) quantile_test = xgb.QuantileDMatrix( Xe_test, ref=quantile_train, ) quantile_model = xgb.train( { "objective": "reg:quantileerror", "quantile_alpha": np.array( [0.05, 0.50, 0.95] ), "tree_method": "hist", "learning_rate": 0.03, "max_depth": 5, "min_child_weight": 8, "subsample": 0.90, "colsample_bytree": 0.90, "seed": SEED, }, quantile_train, num_boost_round=1_400, evals=[(quantile_valid, "validation")], early_stopping_rounds=100, verbose_eval=False, ) quantiles = quantile_model.predict(quantile_test) lower = quantiles[:, 0] median = quantiles[:, 1] upper = quantiles[:, 2] coverage = np.mean( (y_test >= lower) & (y_test <= upper) ) crossing_rate = np.mean(lower > upper) median_width = np.median(upper - lower) print(f"Coverage: {coverage:.1%}") print(f"Crossing rate: {crossing_rate:.1%}") print(f"Median width: {median_width:,.0f}")
The nominal 90% band achieved 82.1% empirical coverage, zero crossings on this test split, and a median width of about 138,106.
That undercoverage is important: conditional quantiles are useful uncertainty bands, not automatically calibrated prediction intervals. Validate coverage by segment and time period before using them for escalation or pricing guardrails.
Coverage and width must be read together. A trivially wide band can cover almost everything while being operationally useless; a narrow band can look attractive while missing expensive tail cases.
Monitor lower-tail misses, upper-tail misses, and crossing rate separately rather than collapsing uncertainty quality into one average.
Finish with an honest stack
Stacking helps only when base models make complementary errors. Train the meta-model on out-of-fold predictions and compare every component on the same final training rows.
from lightgbm import LGBMRegressor from sklearn.compose import ColumnTransformer from sklearn.ensemble import StackingRegressor from sklearn.linear_model import RidgeCV from sklearn.pipeline import make_pipeline from sklearn.preprocessing import OneHotEncoder, StandardScaler best_rounds = engineered.best_iteration + 1 xgb_base = xgb.XGBRegressor( objective="reg:squarederror", n_estimators=best_rounds, learning_rate=0.035, max_depth=4, min_child_weight=8, subsample=0.85, colsample_bytree=0.85, reg_lambda=10.0, reg_alpha=0.05, tree_method="hist", random_state=SEED, n_jobs=-1, ) lgbm_base = LGBMRegressor( objective="regression", n_estimators=650, learning_rate=0.035, num_leaves=31, min_child_samples=40, colsample_bytree=0.85, reg_lambda=10.0, random_state=SEED, n_jobs=-1, verbosity=-1, ) categorical_columns = ["zone_id"] numeric_columns = [ column for column in Xe_train.columns if column not in categorical_columns ] linear_base = make_pipeline( ColumnTransformer( [ ( "numeric", StandardScaler(), numeric_columns, ), ( "categorical", OneHotEncoder( handle_unknown="ignore" ), categorical_columns, ), ] ), RidgeCV( alphas=np.logspace(-2, 4, 13) ), ) stack = StackingRegressor( estimators=[ ("xgb", xgb_base), ("lgbm", lgbm_base), ("ridge", linear_base), ], final_estimator=RidgeCV( alphas=np.logspace(-3, 3, 13) ), cv=3, passthrough=False, n_jobs=1, ) X_fit = pd.concat( [Xe_train, Xe_valid], ignore_index=True, ) y_fit = np.concatenate([y_train, y_valid]) models = { "XGBoost": xgb_base, "LightGBM": lgbm_base, "Ridge": linear_base, "Stack": stack, } for name, model in models.items(): model.fit(X_fit, y_fit) prediction = model.predict(Xe_test) print( name, f"{root_mean_squared_error(y_test, prediction):,.0f}", )
The seeded reference run produced:
| Model | Test RMSE |
|---|---|
| XGBoost | 57,299 |
| LightGBM | 57,243 |
| Ridge | 55,831 |
| Stack | 55,630 |
The stack beat XGBoost by 2.91%, but beat the strongest base model by only 0.36%.
That is the honest delta. The linear model was unusually strong because much of the synthetic target was additive; the stack retained that strength and corrected a small amount of residual nonlinearity.
The ablation order that saves time
After ordinary tuning plateaus:
-
Freeze the split, metric, and baseline.
-
Add one leakage-safe feature family at a time.
-
Add monotonic and interaction rules only where defensible.
-
Use sample weights when segment errors have different costs.
-
Test native robust objectives before maintaining derivatives.
-
Measure quantile coverage, width, and crossings—not just pinball loss.
-
Stack only after out-of-fold residuals show complementary errors.
Take your current best model and run these as seven separate ablations. Keep a change only when it improves the metric, the behavior, or the business cost you actually care about—and record the real delta, including when it is nearly zero.