1. The Intuition: The Relay Putting Team

1.1 The Golf Putting Team Analogy

Imagine a relay golf team trying to sink a ball into a hole located exactly 100 meters away. Instead of having one star golfer hit the ball repeatedly from the starting tee, the team works sequentially in a cooperative chain:

• Golfer 1 (Initial Model $F_0$): Hits the ball 80 meters down the fairway. The remaining distance (residual error) is $+20$ meters.
• Golfer 2 (Weak Learner $h_1$): Walks to where Golfer 1's ball landed. Instead of aiming at the original tee, Golfer 2 aims specifically at the $+20$ meter residual error, hitting the ball 16 meters. The new remaining residual error is $+4$ meters.
• Golfer 3 (Weak Learner $h_2$): Aims at the $+4$ meter error, hitting the ball 3.8 meters, leaving a final error of $+0.2$ meters.

By scaling each golfer's stroke by a learning rate $\eta = 0.1$ to prevent overshooting, the team's combined prediction $F_M(x) = F_0(x) + \eta h_1(x) + \eta h_2(x) + \dots + \eta h_M(x)$ achieves pinpoint accuracy. Gradient Boosting operates on this exact principle: each new decision tree is trained to predict the negative gradient (the residual errors) of the ensemble's previous predictions.

Diagram: Additive sequential boosting workflow training decision trees on pseudo-residuals.

Pitfall — Confusing Bagging (Random Forests) with Boosting (GBDT): Random Forests build deep, independent trees in parallel using bootstrap samples (bagging) to reduce variance. Gradient Boosting builds shallow, sequential trees serially on residual errors (boosting) to reduce bias. Training deep trees in GBDT causes immediate, catastrophic overfitting.

2. Mathematical Formulation of Gradient Boosting

2.1 Loss Minimization & Pseudo-Residuals

Given a dataset $\{(x_i, y_i)\}_{i=1}^n$ and a differentiable loss function $L(y, F(x))$, our goal is to find an additive ensemble model $F(x)$ that minimizes total empirical risk:

At iteration $m$, the model is updated as $F_m(x) = F_{m-1}(x) + \eta \cdot h_m(x)$. To determine the optimal direction for $h_m(x)$, we calculate the negative gradient of the loss function with respect to the current model's predictions, termed the **pseudo-residuals** $r_{im}$:

2.2 Mean Squared Error (MSE) vs Log-Loss Derivations

For Mean Squared Error loss $L(y, F) = \frac{1}{2}(y - F)^2$:

Under MSE loss, pseudo-residuals simplify exactly to the standard residual difference $y_i - F_{m-1}(x_i)$. For binary classification using Log-Loss $L(y, F) = -y \log(p) - (1-y) \log(1-p)$ where $p = \frac{1}{1 + e^{-F}}$, the pseudo-residual is $r_{im} = y_i - p_i$.

2.3 Robust Regression with Huber Loss

When training datasets contain extreme outliers, Mean Squared Error ($L_2$) loss penalizes outliers quadratically, pulling decision tree splits far off-target. Mean Absolute Error ($L_1$) loss is robust to outliers but has a non-differentiable discontinuity at 0. **Huber Loss** combines the best of both worlds by acting quadratically for small errors and linearly for large errors:

The corresponding pseudo-residuals $r_{im}$ for Huber loss switch dynamically based on threshold $\delta$:

By clipping large pseudo-residuals to $\pm \delta$, Huber loss prevents outliers from dominating tree node splits while maintaining smooth convergence for small residuals.

Pitfall — Setting learning rate (η) too high: Setting $\eta = 1.0$ turns off shrinkage, causing the model to over-correct on early noisy samples. Always use a small learning rate ($\eta \in [0.01, 0.1]$) combined with early stopping based on validation loss.

3. Second-Order Taylor Expansion & XGBoost Objectives

3.1 The Second-Order Taylor Series Approximation

Standard GBDT uses only first-order gradients. XGBoost improves optimization speed by taking a **second-order Taylor expansion** of the objective function around the previous prediction $F_{m-1}(x_i)$:

Where $g_i$ is the first-order gradient and $h_i$ is the second-order Hessian:

3.2 Leaf Weight Optimization with Regularization

Removing constant terms $L(y_i, F_{m-1}(x_i))$ and adding L2 tree regularization $\Omega(f) = \gamma T + \frac{1}{2}\lambda \sum_{j=1}^T w_j^2$ yields the simplified objective for leaf node $j$ containing sample set $I_j$:

Differentiating with respect to leaf weight $w_j$ and setting to 0 yields the **Optimal Leaf Weight $w_j^*$**:

3.3 The Split Gain Equation

To evaluate whether a candidate split divides leaf node $I$ into left child $I_L$ and right child $I_R$, XGBoost calculates the **Split Gain $G$**:

If the Gain $G < 0$, the split is rejected and the node is pruned, directly controlled by the minimum gain complexity parameter $\gamma$.

Pitfall — Setting lambda (L2 regularization) to zero: Without L2 regularization ($\lambda = 0$), leaf nodes with very few samples and near-zero Hessians ($\sum h_i \to 0$) produce enormous leaf weights ($w_j^* \to \infty$), destabilizing training. Keep $\lambda \ge 1.0$.

4. Complete Python Implementation from Scratch

Below is a complete, object-oriented Python implementation of a Gradient Boosting Regressor built from scratch using NumPy and a custom DecisionTreeRegressor:

5. Handling Categorical Features & Missing Values

5.1 Default Direction for Missing Data in XGBoost

In traditional decision trees, missing values must be imputed before training. XGBoost handles missing values automatically during split-finding. For each candidate split, missing values are routed to the left child and the gain is calculated, then routed to the right child and calculated. XGBoost permanently assigns the default direction that yields the maximum gain.

5.2 Target Encoding & CatBoost Ordered Boosting

One-hot encoding high-cardinality categorical features creates sparse matrices that severely slow down tree building. CatBoost solves this using **Ordered Target Encoding**, calculating mean target values on historical permutations of data to eliminate target leakage.

5.3 CatBoost Ordered Target Encoding Mathematics

When converting high-cardinality categorical features into numerical values, naive target encoding replaces category $k$ with the mean target value of all rows matching $k$. This introduces severe target leakage and overfitting. CatBoost resolves this by calculating an **Ordered Target Encoder** over randomized permutations $\sigma$ of the training dataset:

Where $D_i$ is the set of historical samples appearing *before* instance $i$ in permutation $\sigma$, $P$ is the global target mean prior, and $a > 0$ is a smoothing parameter (usually $a = 1$). Because sample $i$'s own target value $y_i$ is excluded from its encoder calculation, target leakage is completely eliminated.

Pitfall — Overfitting with One-Hot Encoding on Tree Models: Applying One-Hot Encoding to a categorical feature with 100 levels creates 100 sparse binary columns. Decision trees struggle to split on sparse 0/1 features because individual binary splits isolate tiny fractions of data. Use target encoding or CatBoost categorical handling instead.

6. Advanced: Histogram-Based Boosting, GOSS, and EFB (LightGBM vs XGBoost)

6.1 Binning Continuous Features into 255 Bins

Evaluating every exact unique threshold value for continuous features costs $O(n \cdot d)$ time per split. **LightGBM** and `HistGradientBoostingRegressor` bin continuous floating-point features into 256 discrete integer bins (fitting inside 8-bit `uint8` RAM).

Instead of sorting raw floats, split-finding builds feature histograms in $O(n)$ time. Furthermore, histograms for a right child node can be computed instantaneously by subtracting the left child histogram from the parent histogram ($\text{Hist}_{\text{Right}} = \text{Hist}_{\text{Parent}} - \text{Hist}_{\text{Left}}$), speeding up training by **10x to 20x**.

6.2 Gradient-based One-Side Sampling (GOSS)

In standard GBDT, all data instances contribute equally to split gain computation. However, instances with small gradients (small pseudo-residuals) are already well-trained by previous trees and provide minimal information gain. **Gradient-based One-Side Sampling (GOSS)** downsamples instances intelligently:

1. Sort instances by absolute value of their gradients $|g_i|$ in descending order.
2. Keep the top $a \times 100\%$ instances with the largest gradients (e.g. top 20%).
3. Randomly sample $b \times 100\%$ instances from the remaining small-gradient instances (e.g. 10% of the remaining 80%).
4. Multiply the small-gradient samples by a constant weight factor $\frac{1 - a}{b}$ when computing split gain to restore the original data distribution without biasing gradient estimates.

6.3 Exclusive Feature Bundling (EFB)

High-dimensional tabular datasets are frequently sparse (e.g. one-hot encoded vectors where features are mutually exclusive, rarely taking non-zero values simultaneously). **Exclusive Feature Bundling (EFB)** reduces feature dimensionality $d$ by bundling mutually exclusive features into single composite feature bins using a graph coloring algorithm, slashing RAM overhead and split search complexity.