Gradient Boosting Regressor

Introduction

Gradient Boosting Regressor is a potent machine learning algorithm that is part of the ensemble methods family, particularly within the boosting techniques category. It operates by sequentially incorporating predictors (commonly decision trees), with each one amending the errors of its predecessor. Unlike other boosting approaches that modify the weights of observations, Gradient Boosting concentrates on minimizing the loss function, a differentiable function that quantifies the model's error, making it highly effective for regression tasks.

Background and Theory

Gradient Boosting combines multiple weak learners (models that perform slightly better than random guessing) to create a strong learner in a sequential manner. The core principle is to construct new models that predict the residuals or errors of prior models and then add these new models to the ensemble.

Mathematical Foundations

Let's denote our data set as $(x_i, y_i)$, where $x_i$ represents the features and $y_i$ the target variable for the $i^\text{th}$ observation. The goal is to find a function $F(x)$ that maps input features $x$ to predictions as close as possible to the target values $y$.

The algorithm starts with a model that makes constant predictions:

where $L(y_i, \gamma)$ is the loss function evaluated at the initial prediction $\gamma$ and the true target value $y_i$.

At each stage $m = 1$ to $M$, where $M$ is the total number of trees, the algorithm fits a new model $h_m(x)$ to the negative gradient of the loss function evaluated at the previous model's predictions:

Then, it finds the best-fit model $h_m(x)$ to these residuals. The model $h_m(x)$ is typically a decision tree.

Next, it computes a multiplier $\gamma_m$ for the model $h_m(x)$ that minimizes the loss when added to the existing model $F_{m-1}(x)$:

The model is then updated:

where $\nu$ is the learning rate, a parameter that scales the contribution of each $h_m(x)$.

The process repeats for $M$ iterations or until a stopping condition is reached, resulting in a final model:

Algorithmic Contemplation

Let’s try generalizing the problem Gradient Boosting algorithm is about to solve. We want to find the function $\hat{F}\in\mathcal{H}$ which minimizes the expected loss function in terms of $p$-dimensional loss function $L$, where $\mathcal{H}=\{F\mid F:\mathbb{R}^p\rightarrow\mathbb{R}\}$.

Unfortunately, $\mathcal{H}$ is an enormous set for the algorithm to find optimal $\hat{F}$. Therefore, some constraints are needed, which constrains $\mathcal{H}$ to:

where $M$ is a natural number.

In short, a structure which can be expressed through a linear combination of a constant and several functions is added, rather than a simple generic function.

Moreover, since the distributions of $x,y$ are unknown, it is highly demanding to figure out the expected loss. Therefore, we have to find $\hat{F}$ that minimizes the experienced loss function, which can be derived by utilizing a finite amount of data $(x_i,y_i),~i=1,\ldots,n$.

This can be found by following several steps:

At this point, Gradient Boosting reforms the last equation by substituting the minimizing part with the gradient descent method as following:

$F_m$ can be approximated with respect to the positive integer $\gamma$. This is the reason why the term Gradient is included in the name of this boosting algorithm.

Here, the gradient can be computed as:

which in most cases(except for MSE loss), are very complitated to compute. Therefore, we take the 1st-order Taylor approximation with respect to $f_m$, which necessarily requires the loss function to be differentiable.

If we differentiate the right-hand term of the above equation with respect to $f_m(x_i)$, we can approximate the $\nabla L(y_i,F_{m-1}(X_i)+f_m(x_i))$ by:

And therefore:

Additionaly, we can optimize $\gamma$ in order to minimize the residuals:

As a result, we can define the pseudo-residual $r_{im}$ as following:

For the next step, we train a new function $g_m$, which accepts $x_i$ as an input variable and $r_{im}$ as an output variable, by approximating $r_{im}$ with respect to $x_i$. Then we can redefine $\gamma_m$ as following:

Consequently, $F_m$ can be finally set as $F_m(x)=F_{m-1}(x)+l\cdot\gamma_m\cdot g_m(x)$ by adopting $l$, which is for the prevention of overfitting.

From previous explanations, $h_i(x),~i=1,\ldots,M$ are in fact $\gamma_i\cdot g_i(x)$.

Loss Functions

In the Gradient Boosting Regressor, the choice of the loss function is crucial as it directly influences the model's performance by determining how well the predictions fit the actual data. Here are three popular loss functions and their mathematical integration into the algorithm:

Mean Squared Error (MSE) Loss

For regression tasks where the target variable $y$ represents the true value, and the model's prediction is denoted by $\hat{y}$, the Mean Squared Error (MSE) loss for a single observation is defined as:

When applying Gradient Boosting for regression, the gradient (negative gradient of the loss with respect to the model's prediction $\hat{y}$ used to fit the next model is:

where $\hat{y}_i$ is the prediction for the $i$-th observation based on the model up to iteration $m-1$.

Mean Absolute Error (MAE) Loss

For the same regression context, the Mean Absolute Error (MAE) loss, which measures the average magnitude of errors between predicted and actual values without considering their direction, for a single observation is:

The gradient for the MAE loss with respect to the model's prediction $\hat{y}$ is:

where $\text{sign}(x)$ is a function that returns 1 if $x > 0$, -1 if $x < 0$, and 0 if $x = 0$, and $\hat{y}_i$ is as previously defined.

Huber Loss

Huber Loss is used in regression and is less sensitive to outliers in data than the squared error loss. It combines the properties of both MSE and MAE. The loss function for a single observation can be defined piecewise as:

where $\delta$ is a hyperparameter that dictates the transition between the squared loss and the absolute loss.

The gradient for the Huber loss with respect to $\hat{y}$ is:

where $\hat{y}_i$ and $\delta$ are as defined above. This formulation allows for adjusting the model's sensitivity to outliers in the gradient boosting process.

Procedural Steps

1. Initialization: Start with a constant model $F_0(x)$.
2. For $m = 1$ to $M$:
   • Compute the residuals $r_{im}$, the negative gradient of the loss function.
   • Fit a model $h_m(x)$ to these residuals.
   • Find the optimal multiplier $\gamma_m$.
   • Update the model $F_m(x) = F_{m-1}(x) + \nu \gamma_m h_m(x)$.
3. Output: The final model $F_M(x)$.

Implementation

Parameters

• base_estimator: Estimator, default = DecisionTreeRegressor()
  Base estimator for training multiple models
  
• n_estimators: int, default = 100
  Number of base estimators to fit
  
• learning_rate: float, default = 0.1
  Shrinking factor of the contribution of each estimator
  
• subsampling: float, default = 1.0
  Fraction of samples to be used for fitting each estimator
  
• loss: Literal['mse', 'mae', 'huber'], default = ‘mse’
  Type of loss function
  
• delta: float, default = 1.0
  Balancing factor between MSE loss and MAE loss
  

Examples

Test on the synthesized dataset of the curve $y=\sin{3x}\cdot\cos{e^{x/2}}+\epsilon$:

from luma.ensemble.boost import GradientBoostingRegressor
from luma.preprocessing.scaler import StandardScaler
from luma.metric.regression import RootMeanSquaredError
from luma.visual.evaluation import ResidualPlot

import matplotlib.pyplot as plt
import numpy as np

np.random.seed(42)

X = np.linspace(-4, 4, 400).reshape(-1, 1)
y = (np.sin(3 * X) * np.cos(np.exp(X / 2))).flatten()
y += 0.15 * np.random.randn(400)

sc = StandardScaler()
y_trans = sc.fit_transform(y)

gb = GradientBoostingRegressor(n_estimators=50,
                               learning_rate=0.1,
                               subsample=1.0,
                               loss='mae',
                               max_depth=3)

gb.fit(X, y_trans)

y_pred = gb.predict(X)
score = gb.score(X, y_trans, metric=RootMeanSquaredError)

fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.scatter(X, y_trans, 
            s=10, c='black', alpha=0.3, 
            label=r'$y=\sin{3x}\cdot\cos{e^{x/2}}+\epsilon$')

for tree in gb:
    ax1.plot(X, tree.predict(X), c='pink', alpha=0.1)

ax1.plot(X, y_pred, lw=2, c='crimson', label='Predicted Plot')
ax1.legend(loc='upper right')
ax1.set_xlabel('x')
ax1.set_ylabel('y (Standardized)')
ax1.set_title(f'GradientBoostingRegressor [RMSE: {score:.4f}]')

res = ResidualPlot(gb, X, y_trans)
res.plot(ax=ax2)

plt.tight_layout()
plt.show()

Applications

Gradient Boosting Classifier is used in a variety of domains due to its flexibility and accuracy. Common applications include but are not limited to:

• Fraud detection in banking.
• Customer churn prediction.
• Disease outbreak prediction.
• Demand forecasting in retail.

Strengths and Limitations

Strengths

• Highly Accurate: Often provides very high accuracy across a wide range of applications.
• Flexibility: Can handle different types of data and is adaptable to various loss functions.

Limitations

• Prone to Overfitting: Especially with noisy data and without proper regularization.
• Computationally Intensive: Can be slower to train due to the sequential nature of boosting.
• Parameter Tuning: Requires careful tuning of parameters and stopping criteria to avoid overfitting and underfitting.

Advanced Topics

• Regularization Techniques: Techniques like shrinkage (learning rate) and stochastic gradient boosting can help prevent overfitting.
• Loss Functions: Beyond the typical use with classification and regression tasks, custom loss functions can be implemented for specific applications.

References

1. Friedman, J. H. (2001). "Greedy Function Approximation: A Gradient Boosting Machine." Annals of Statistics.
2. Natekin, A., & Knoll, A. (2013). "Gradient boosting machines, a tutorial." Frontiers in Neurorobotics.