Gradient Boosting

Open In Colab

The foundations of gradient boosting were laid by Jerome Friedman in 2001 in his paper Greedy Function Approximation: A Gradient Boosting Machine, published in the Annals of Statistics. His work united two ideas. The first was boosting, developed in the early 1990s as a way to combine many weak learners so that each new learner focuses on correcting the errors of the previous ones. Early boosting methods, such as AdaBoost, were mainly designed for classification and relied on reweighting misclassified samples. The second idea was gradient descent. Friedman showed that boosting can be viewed as numerical optimization, allowing the method to work with any differentiable loss and to handle both regression and classification.

For decision trees this results in a clear contrast with random forests. A random forest trains many trees independently and averages their predictions. Gradient boosting trains trees one after another, where each new tree attempts to fit the residuals of the current forest.

The widespread adoption of gradient boosting accelerated after the release of the XGBoost library in 2014 by Tianqi Chen, who was a PhD student at the University of Washington at the time. Its speed on large datasets was the key factor that made it the dominant approach for tabular machine learning competitions. In 2016 Chen and his Ph.D. advisor published the paper XGBoost: A Scalable Tree Boosting System, which described the algorithm and documented its success across many tasks.

“The most important factor behind the success of XGBoost is its scalability in all scenarios. The system runs more than ten times faster than existing popular solutions on a single machine and scales to billions of examples in distributed or memory-limited settings.”

Later, other high performing implementations appeared, including LightGBM from Microsoft Research Asia and CatBoost from Yandex, each introducing their own innovations.

In 2025, researchers from various institutes introduced TabArena, a large benchmark suite for evaluating algorithms on tabular data. Results from TabArena show that gradient boosting methods are one of the best methods on tabular datasets.

Today, you can even find job postings that list specific gradient boosting frameworks as required skills. For example, a machine learning engineer vacancy from December 2025 asked for:

Familiarity with ML libraries: TensorFlow, XGBoost, scikit-learn, etc.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

1. Gradient Boosting from Decision Trees

Like we did in the random forests notebook, let’s look at how gradient boosting combines multiple decision trees.

  • We will start with 3 decision trees
  • This time we will use a toy dataset from scikit-learn because visualizing individual splits is not the goal.
  • We will also use regression instead of classification because the gradient boosting formulation is easier to illustrate when we work directly with numerical residuals and gradient boosting classifiers rely on regression trees internally.
from sklearn.datasets import load_diabetes
X, y = load_diabetes(as_frame=True, return_X_y=True)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
# we will save the trees here
boosted_forest = []

# let's keep it simple
max_depth = 3
# nothing new here
from sklearn.tree import DecisionTreeRegressor

dt1 = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
dt1.fit(X_train, y_train) # just fit the tree

y_pred = dt1.predict(X_train) # get preds

boosted_forest.append(dt1) # save the first tree to the forest

Let’s check the errors (i.e. residuals) of the tree. Remember we covered residuals in 007_modeling_and_evaluation.ipynb section 9. Evaluating Regression Models.

y_residual = y_train - y_pred
pred_df = pd.DataFrame({
    'y_train': y_train,
    'tree1_pred': y_pred,
    'residual_1': y_residual
    }).reset_index(drop=True)
pred_df

So far everything looks straightforward. The next step is to fit the second tree on the residuals in order to correct the remaining errors. The residuals represent how far the first tree’s predictions are from the true values, so training the next tree on these residuals means we want it to learn the mistakes of the first tree and reduce them.

dt2 = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
dt2.fit(X_train, y_residual) # the second tree is fit on the y_residual instead of y_train
y_pred = dt2.predict(X_train)
boosted_forest.append(dt2)

When we moved from the first tree to the second tree, we computed the residuals directly as y_train - y_pred because no corrections had been made yet. From this point onward the situation changes. Each new tree is meant to fix whatever error remains after all previous trees, so we no longer compare to the original targets. Instead, we update the residuals by subtracting the predictions of the newest tree from the current residuals. This way

y_residual = y_residual - y_pred

ensures that the third tree (and every tree after it) focuses only on the part of the error that previous trees have not yet corrected.

y_residual = y_residual - y_pred # instead of y_residual = y_train - y_pred
pred_df['tree2_pred'] = y_pred
pred_df['residual_2'] = y_residual.values
pred_df
# repeat
dt3 = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
dt3.fit(X_train, y_residual)

y_pred = dt3.predict(X_train)
boosted_forest.append(dt3)
y_residual = y_residual - y_pred
pred_df['tree3_pred'] = y_pred
pred_df['residual_3'] = y_residual.values
pred_df

We now have a gradient boosting regressor made of three decision trees. To make a prediction, the model adds together the contributions of all three trees, so the final output is tree_1 + tree_2 + tree_3. Let’s inspect this process on a single datapoint to see how each tree contributes to the final prediction.

test_datapoint = X_test.iloc[0:1]
test_datapoint
# prediction of the first tree
p1 = dt1.predict(test_datapoint)
p1
# prediction of the second tree
p2 = dt2.predict(test_datapoint)
p2
# prediction of the third tree
p3 = dt3.predict(test_datapoint)
p3
# the final prediction
p1+p2+p3
# the label
y_test.iloc[0]
# for all X_test
y_pred = boosted_forest[0].predict(X_test) + boosted_forest[1].predict(X_test) + boosted_forest[2].predict(X_test)
y_pred
# we can make this work for any number of trees in the forest
y_pred = sum(dt.predict(X_test) for dt in boosted_forest)
y_pred

and calculate the MAE:

from sklearn.metrics import mean_absolute_error
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

let’s put this into a loop where we can change the number of trees easily:

n_estimators = 3
boosted_forest = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    # for the next iteration
    y_pred = dt.predict(X_train)
    y_residual = y_residual - y_pred
# the same mae is expected
y_pred = sum(dt.predict(X_test) for dt in boosted_forest)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

Now we can add more trees:

n_estimators = 50 # the only change
boosted_forest = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    # next iteration
    y_pred = dt.predict(X_train)
    y_residual = y_residual - y_pred

y_pred = sum(dt.predict(X_test) for dt in boosted_forest)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

It got worse 😨😨😨 But we learned in random forests that adding more trees should generally improve the ensemble, so why is this happening? The short answer is random forests are independent trees whose errors average out, whereas in gradient boosting trees grow sequentially and each one modifies the previous ensemble. To answer in depth, we need to understand the role of the learning rate and update our implementation accordingly. Before doing that, let’s take a look at how to use gradient boosting in scikit-learn.

2. Gradient Boosting in scikit-learn

Gradient boosting in scikit-learn is available through two main estimators: GradientBoostingClassifier and GradientBoostingRegressor. These models implement Friedman’s gradient boosting framework.

For larger datasets (n_samples >= 10000), scikit-learn provides much faster histogram-based implementations. These are HistGradientBoostingClassifier and HistGradientBoostingRegressor. They are not modified versions of the two estimators above; instead, they are separate estimators inspired by LightGBM’s histogram-based approach. This design makes them considerably more scalable and efficient, and they also bring additional functionality such as native handling of missing values.

from sklearn.ensemble import GradientBoostingRegressor
gbr = GradientBoostingRegressor(random_state=0, n_estimators=3)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

As you can see, using a new algorithm in scikit-learn works just like any other estimator. However, becoming proficient with gradient boosting requires understanding the ideas behind it and gaining hands-on experience. The more you engage with it and experiment on your own, the faster the concepts will click and the better your models will become.

Now let’s look at what happens when we add more estimators. In our implementation the model got worse, because it was missing a key piece.

# n_estimators 10
gbr = GradientBoostingRegressor(random_state=0, n_estimators=10)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)
# n_estimators 30
gbr = GradientBoostingRegressor(random_state=0, n_estimators=30)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)
# n_estimators 50
gbr = GradientBoostingRegressor(random_state=0, n_estimators=50)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

much better! Now let’s repeat this with a for loop, ranging n_estimators from 1 to 100 and plot:

n_estimators_list = []
test_error = []
for n_estimators in range(1, 100, 3):
    gbr = GradientBoostingRegressor(random_state=0, n_estimators=n_estimators)
    gbr.fit(X_train, y_train)
    y_pred = gbr.predict(X_test)
    mae = mean_absolute_error(y_test, y_pred)
    n_estimators_list.append(n_estimators)
    test_error.append(mae)
plt.plot(n_estimators_list, test_error)
plt.xlabel('n_estimators')
plt.ylabel('Test MAE')

As you can see, the test error decreases up to about 20–30 estimators and then begins to rise again. This behavior must be familiar to you at this point: overfitting! Let’s plot this together with the training error:

n_estimators_list = []
test_error = []
training_error = []
for n_estimators in range(1, 100, 3):
    gbr = GradientBoostingRegressor(random_state=0, n_estimators=n_estimators)
    gbr.fit(X_train, y_train)
    y_pred = gbr.predict(X_test)
    mae = mean_absolute_error(y_test, y_pred)
    n_estimators_list.append(n_estimators)
    test_error.append(mae)
    training_error.append(mean_absolute_error(y_train, gbr.predict(X_train)))
plt.plot(n_estimators_list, test_error, label='test_error')
plt.plot(n_estimators_list, training_error, label='training_error')
plt.xlabel('n_estimators')
plt.ylabel('Test MAE')
plt.legend()

A masterpiece of a plot. We need to understand the hyperparameters next to prevent overfitting.

YOUR TURN

Try GradientBoostingClassifier on the breast cancer dataset:

from sklearn.datasets import load_breast_cancer
# YOUR CODE HERE

3. Hyperparameters

Most hyperparameters will look familiar from decision trees and random forests. Gradient boosting adds a few parameters that control how each tree corrects the errors of the ensemble. Different libraries implement these ideas in slightly different ways, so aim to understand their role conceptually rather than tuning them mechanically. We will cover:

  • Learning Rate
  • Early Stopping
  • Loss

3.1 Learning Rate

Let’s get back to the problem we encountered previously, where adding trees increased the test error in our initial implementation:

n_estimators = 50
boosted_forest = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    # next iteration
    y_pred = dt.predict(X_train)
    y_residual = y_residual - y_pred

y_pred = sum(dt.predict(X_test) for dt in boosted_forest)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

The reason this implementation fails is that we have not included a learning rate, so the model updates in large jumps rather than taking the small, controlled steps needed to avoid overshooting the optimum. The key idea is to reduce the influence of each tree by scaling its contribution with a step size. This works exactly like gradient descent and can be implemented as:

y_residual = y_residual - learning_rate * y_pred
learning_rate = 0.1 # ADD lr
n_estimators = 50
boosted_forest = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    y_pred = dt.predict(X_train)
    y_residual = y_residual - learning_rate * y_pred # ADD lr
y_pred = sum(learning_rate * dt.predict(X_test) for dt in boosted_forest) # ADD lr
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

Suddenly the results are much closer to the scikit-learn implementation 🥳🥳 This shows that the core idea behind boosting is conceptually simple, even if production implementations add many refinements.

# learning rate = 1
# this is equal to our initial non-learning rate implementation
learning_rate = 1
n_estimators = 50
boosted_forest = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    y_pred = dt.predict(X_train)
    y_residual = y_residual - learning_rate * y_pred
y_pred = sum(learning_rate * dt.predict(X_test) for dt in boosted_forest)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

Now let’s understand the influence of changing the learning rate, by repeating the test error vs n_estimators plot at different learning rates:

plot_dict = {}

for learning_rate in [0.01, 0.03, 0.1, 0.3, 1]:

    x_for_plotting = []
    y_for_plotting = []
    for n_estimators in range(1, 100, 3):
        gbr = GradientBoostingRegressor(random_state=0, n_estimators=n_estimators, learning_rate=learning_rate)
        gbr.fit(X_train, y_train)
        y_pred = gbr.predict(X_test)
        mae = mean_absolute_error(y_test, y_pred)
        x_for_plotting.append(n_estimators)
        y_for_plotting.append(mae)

    plot_dict[learning_rate] = (x_for_plotting, y_for_plotting)
annotation_location = {0.01: 15, 0.03: -5, 0.1: -5, 0.3: -5, 1: -5}
for lr in plot_dict:
    x_vals, y_vals = plot_dict[lr]
    plt.plot(x_vals, y_vals)
    # annotate lr on the curves
    plt.text(x_vals[annotation_location[lr]], y_vals[annotation_location[lr]], f'lr={lr}', fontsize=10, bbox=dict(facecolor='white', alpha=0.6, edgecolor='black'))
    
plt.xlabel('n_estimators')
plt.ylabel('Test MAE')

The best way to read this plot is to track how the behavior changes as the learning rate decreases from 1 to 0.01.

  • For lr = 1 the error shoots up immediately, which tells us the step size is far too large.
  • For lr = 0.3 the curve improves slightly but still overshoots, so it remains too high.
  • For lr = 0.1 we finally get a smooth trajectory with a clear minimum around 20–30 trees, although the learning rate may still be a bit high.
  • For lr = 0.03 the error decreases more slowly and appears to reach an even lower value than the lr = 0.1 optimum.
  • For lr = 0.01 the curve looks poor only because we stopped at 100 trees; with more boosting rounds it might reach a better solution.

This illustrates the tradeoff between the learning rate and the number of estimators. Smaller learning rates often lead to better generalization but require more trees and longer training times. These two hyperparameters must be tuned together. Choosing 100 trees for lr = 0.01 leads to underfitting because the model needs many more steps, whereas choosing 100 trees for lr = 0.1 leads to overfitting because the steps are too large.

3.2 Early Stopping

The key idea to tune lr and n_estimators together is to use early stopping. Set n_estimators to a large value, monitor the validation error after each tree, and stop when the error starts to rise. However, it is not sensible to stop the moment a single tree increases the error. The curve is not perfectly smooth and small temporary increases can occur even when the overall trend is improving. A better approach is to use patience: allow the model to continue for several more iterations, for example 20, and stop only if the error has not improved during that window. This is called early stopping because training ends before reaching the maximum number of estimators. For this reason it is best to set n_estimators high and rely on early stopping to find the right point. If training actually reaches n_estimators without stopping, the model is likely underfitting, as in the case we have seen above with lr=0.01 and n_estimators=100.

# we need a validation set
# instead of tracking the test error
# you should know why at this point
X_dev, X_test, y_dev, y_test = train_test_split(X, y, test_size=1/10, random_state=0)
X_train, X_val, y_train, y_val = train_test_split(X_dev, y_dev, test_size=1/9, random_state=0)
# let's repeat what we have done
# only replace the test data with the validation data
learning_rate = 0.1
n_estimators = 100
boosted_forest = []
mae_values = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    y_pred = dt.predict(X_train)
    y_residual = y_residual - learning_rate * y_pred

    y_pred = sum(learning_rate * dt.predict(X_val) for dt in boosted_forest)
    mae = mean_absolute_error(y_val, y_pred)
    mae_values.append(mae)
n_estimators_list = list(np.arange(len(mae_values))+1)
plt.plot(n_estimators_list, mae_values)
plt.ylabel('Validation MAE')
plt.xlabel('n_estimators')
# This is where we want to stop
best_val_mae = min(mae_values)
best_n_estimators = n_estimators_list[mae_values.index(best_val_mae)]

print(f"Best validation MAE: {best_val_mae:.2f}")
print(f"Best n_estimators: {best_n_estimators}")

a closer look to the vicinity:

plt.plot(n_estimators_list, mae_values)
plt.xlim(best_n_estimators-15,best_n_estimators+5)
plt.ylim(best_val_mae-1,best_val_mae+5)
plt.scatter(best_n_estimators, best_val_mae, color='red')

here is how we can add early stopping, changes marked by comments:

learning_rate = 0.1
n_estimators = 500
boosted_forest = []
mae_values = []
y_residual = y_train.copy()

# Addition for early stopping
patience = 20  # num iterations to wait without improvement
best_mae = float("inf") # keep track of the best mae so far
best_iter = -1 # keep track of the best iter so far
no_improve_count = 0 # if this value hits patience we break

for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    y_pred_train = dt.predict(X_train)
    y_residual = y_residual - learning_rate * y_pred_train

    y_pred_val = sum(learning_rate * tree.predict(X_val) for tree in boosted_forest)
    mae = mean_absolute_error(y_val, y_pred_val)
    mae_values.append(mae)

    # Addition for early stopping
    if mae < best_mae:
        best_mae = mae
        best_iter = i
        no_improve_count = 0
    else:
        no_improve_count += 1
        if no_improve_count >= patience:
            print(f"Early stopping at iteration {i} with best iteration {best_iter} and MAE={best_mae:.4f}")
            break

Keep in mind that when we use early stopping with patience, we do not actually stop at the minimum. Because we use patience, training continues past the optimum and only stops after the model fails to improve for the amount of rounds equal to patience.

# MAE will be higher
y_pred = sum(learning_rate*dt.predict(X_val) for dt in boosted_forest)
mae = mean_absolute_error(y_val, y_pred)
print("MAE:", mae)

so we need to slice the forest until the optimum point:

# keep only the trees up to the best iteration
best_n_estimators = best_iter + 1
boosted_forest = boosted_forest[:best_n_estimators] # slice until the optimum
mae_values = mae_values[:best_n_estimators]

print(f"Using {best_n_estimators} trees with best validation MAE = {best_mae:.4f}")
y_pred = sum(learning_rate*dt.predict(X_val) for dt in boosted_forest)
mae = mean_absolute_error(y_val, y_pred)
print("MAE:", mae)

final step would be to report the test error.

# Test error
y_pred = sum(learning_rate*dt.predict(X_test) for dt in boosted_forest)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

3.3 Loss

We mentioned above that the foundations of gradient boosting were laid by Jerome Friedman in 2001 in his paper Greedy Function Approximation: A Gradient Boosting Machine. In the abstract he says:

Function estimation/approximation is viewed from the perspective of numerical optimization in function space, rather than parameter space. A connection is made between stagewise additive expansions and steepest descent minimization. A general gradient descent “boosting” paradigm is developed for additive expansions based on any fitting criterion. Specific algorithms are presented for least-squares, least absolute deviation, and Huber-M loss functions for regression, and multiclass logistic likelihood for classification.

Gradient descent, the optimization method behind many machine learning models, tries to find the best parameters (a search in the parameter space).

\[new\_parameters = old\_parameters − lr × gradient\_of\_the\_loss\_at\_the\_old\_parameters\]

The analogous rule for gradient boosting (the update in function space) becomes:

\[new\_function = old\_function − lr × tree\_that\_approximates\_the\_gradient\_of\_the\_loss\]

How is loss a hyperparameter in our implementation?

learning_rate = 0.1
n_estimators = 100
boosted_forest = []
mae_values = []
y_residual = y_train.copy()
for i in range(n_estimators):
    dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth)
    dt.fit(X_train, y_residual)
    boosted_forest.append(dt)
    y_pred = dt.predict(X_train)
    
    # HERE IT IS
    y_residual = y_residual - learning_rate * y_pred # GRADIENT OF THE LOSS

    y_pred = sum(learning_rate * dt.predict(X_val) for dt in boosted_forest)
    mae = mean_absolute_error(y_val, y_pred)
    mae_values.append(mae)

the loss we implemented is the squared error (you are seeing the gradient in the code). In scikit-learn this is the default loss:

gbr = GradientBoostingRegressor(random_state=0, n_estimators=50, loss='squared_error')
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)

Check the documentation for available loss functions:

loss{‘squared_error’, ‘absolute_error’, ‘huber’, ‘quantile’}, default=’squared_error’

Loss function to be optimized. ‘squared_error’ refers to the squared error for regression. ‘absolute_error’ refers to the absolute error of regression and is a robust loss function. ‘huber’ is a combination of the two. ‘quantile’ allows quantile regression (use alpha to specify the quantile). See Prediction Intervals for Gradient Boosting Regression for an example that demonstrates quantile regression for creating prediction intervals with loss=‘quantile’.

gbr = GradientBoostingRegressor(random_state=0, n_estimators=50, loss='huber')
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", mae)