Gradient Boosting

Imagine you're trying to predict house prices. Your first model predicts $200k for every house (just using the mean). Obviously this is terrible—some houses should be $100k, others $500k.

Instead of starting over with a completely different model, what if you built a second model that predicts the errors of the first model? Then your combined prediction would be:

The second model doesn't try to predict house prices directly. It predicts “how much did the first model miss by?” If the first model predicted $200k for a house actually worth $300k, the second model learns to predict +$100k for that type of house.

Now you have a better model. But it's still not perfect. So you add a third model that predicts the errors of (Model₁ + Model₂). And a fourth that predicts the errors of (Model₁ + Model₂ + Model₃). Each model incrementally improves the ensemble by focusing on what previous models got wrong.

Gradient boosting can work with any base learner—linear models, neural networks, even other ensembles—but decision trees are the standard choice. Why?

In gradient boosting, we typically use shallow regression trees(even for classification, as we'll see). These are stumps or small trees that output a single numeric value at each leaf. Deep trees would overfit and wouldn't ensemble well.

The algorithm proceeds in steps. First, we initialize with a constant prediction:

1def fit(self, X, y):
2    self.initial_prediction = np.mean(y)
3    F = np.full(len(y), self.initial_prediction)  # Current predictions
4    self.trees = []
5    
6    for m in range(self.n_estimators):
7        # Step 1: Compute residuals
8        residuals = y - F
9        
10        # Step 2: Fit tree to residuals
11        tree = DecisionTreeRegressor(max_depth=self.max_depth)
12        tree.fit(X, residuals)
13        self.trees.append(tree)
14        
15        # Step 3: Update predictions
16        predictions = tree.predict(X)
17        F += self.learning_rate * predictions

Let's break down what happens at each iteration:

After M iterations, the final ensemble prediction is:

1def predict(self, X):
2    # Start with the initial mean prediction
3    F = np.full(len(X), self.initial_prediction)
4    
5    # Add each tree's contribution
6    for tree in self.trees:
7        F += self.learning_rate * tree.predict(X)
8    
9    return F

Here's where the “gradient” comes from. Gradient boosting is gradient descent in function space. In normal gradient descent, we minimize a loss L(w) by updating parameters:

In gradient boosting, we minimize a loss L(F) by updating the function F:

∂L/∂F isn't a number—it's a vector, one value for each training sample. This is the negative gradient of the loss with respect to the predictions. For MSE loss L = ½(y − F)²:

The beauty of gradient boosting is you can use any differentiable loss function. Just compute its gradient with respect to predictions.

Each loss function gives different residuals to fit. MSE focuses on all errors equally. MAE is robust to outliers because large errors don't dominate. Huber smoothly transitions between the two.

Classification is trickier because labels are discrete (0 or 1), but we need continuous predictions to fit trees to. The trick: predict log-odds instead of probabilities, and use log-loss (binary cross-entropy).

1class GradientBoostingClassifier:
2    def _sigmoid(self, x):
3        return 1 / (1 + np.exp(-np.clip(x, -500, 500)))
4    
5    def fit(self, X, y):
6        # Initialize with log-odds
7        p_positive = np.mean(y)
8        self.initial_prediction = np.log(p_positive / (1 - p_positive))
9        F = np.full(len(y), self.initial_prediction)
10        
11        for m in range(self.n_estimators):
12            # Convert log-odds to probabilities
13            probabilities = self._sigmoid(F)
14            
15            # Negative gradient of log-loss = y - p
16            residuals = y - probabilities
17            
18            # Fit tree to negative gradient
19            tree = DecisionTreeRegressor(max_depth=self.max_depth)
20            tree.fit(X, residuals)
21            self.trees.append(tree)
22            
23            # Update log-odds
24            F += self.learning_rate * tree.predict(X)
25    
26    def predict_proba(self, X):
27        F = np.full(len(X), self.initial_prediction)
28        for tree in self.trees:
29            F += self.learning_rate * tree.predict(X)
30        return self._sigmoid(F)
31    
32    def predict(self, X):
33        return (self.predict_proba(X) > 0.5).astype(int)

For log-loss L = −[y log(p) + (1−y) log(1−p)], the negative gradient is simply y − p, the difference between true labels and predicted probabilities. This is remarkably similar to MSE residuals!

Gradient boosting dominates machine learning competitions and real-world applications for several fundamental reasons:

Gradient boosting has many hyperparameters, but understanding their interactions makes tuning manageable.

Our implementation is the foundation, but production libraries add many optimizations that make them 10-100x faster and more accurate: