Gradient Boosting Regression from Scratch
Build a gradient boosting regressor from scratch using NumPy.
Introduction
30 Second Summary
Every time you see a house price estimate online or a weather forecast on your phone, there is likely a gradient boosting model working behind the scenes. These models are the workhorses of modern prediction tasks, dominating competitions and production systems alike.
In this project, you will build a complete Gradient Boosting Regressor from scratch using Python and NumPy, then validate it against scikit-learn's implementation on real California housing data.
What You'll Build
A from-scratch implementation of Friedman's Gradient Boosting algorithm, including a custom Decision Tree Regressor, evaluation metrics, and a comparison pipeline against scikit-learn.
By the end of this project, you'll have:
- A from-scratch Decision Tree Regressor with recursive MSE-based splitting that serves as the weak learner for boosting.
- A complete Gradient Boosting Regressor that iteratively fits pseudo-residuals to make predictions on 20,640 real housing samples.
- A side-by-side comparison with scikit-learn's implementation, validating your model's correctness through feature importance rankings and R-squared scores.
- Secret Mission: Implement early stopping with validation monitoring and visualize learning curves showing where overfitting begins.
Are there any prerequisites?
This project assumes solid Python and NumPy proficiency, a working understanding of decision trees, and familiarity with train/test evaluation. You should be comfortable reading and writing NumPy array operations.
Set Up Data and Metrics
Before building any machine learning model, you need data to train on and metrics to measure performance.
In this step, you'll load a real housing dataset, build a reproducible train/test split from scratch, and write two evaluation functions that will score every model in this project.
In this step, get ready to:
- Load the California Housing dataset and implement a train/test split using NumPy.
- Implement custom evaluation metrics for measuring prediction quality.
Add imports and write evaluation metrics
The top of your file needs all the libraries you'll use throughout the project. You'll also define two metric functions that quantify how good (or bad) a model's predictions are.
- In gradient_boosting.py, add the imports and evaluation metric functions by pasting the following code:
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor as SklearnGBR
import matplotlib.pyplot as plt
# ============================================================
# Evaluation Metrics
# ============================================================
def mean_squared_error(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
return 1 - (ss_res / ss_tot)What does this code do?
The imports pull in four libraries you'll use across the project:
- numpy provides the array math for building trees and computing metrics.
- fetch_california_housing loads the dataset you'll train on.
- SklearnGBR is scikit-learn's implementation that you'll compare against later.
- matplotlib will render comparison charts in a later step.
mean_squared_error computes the average squared difference between predictions and actual values. r_squared measures the proportion of variance your model explains, where 1.0 is a perfect fit and 0.0 means the model is no better than guessing the mean.
You won't see a change yet. The next chunks use these imports and metrics.
Implement the train/test split
A train/test split divides your dataset into two groups: one the model learns from, and one it has never seen before. You'll build this yourself using NumPy's random number generator for reproducible shuffling.
- Below the r_squared function, add the train/test split function by pasting this code:
# ============================================================
# Train/Test Split
# ============================================================
def train_test_split(X, y, test_size=0.2, random_state=42):
rng = np.random.default_rng(random_state)
n_samples = X.shape[0]
indices = rng.permutation(n_samples)
test_count = int(n_samples * test_size)
test_idx = indices[:test_count]
train_idx = indices[test_count:]
return X[train_idx], X[test_idx], y[train_idx], y[test_idx]What does this code do?
- np.random.default_rng(random_state) creates a reproducible random number generator. Using the same seed (42) guarantees the same shuffle every time you run the script.
- rng.permutation(n_samples) shuffles the row indices so the split is random, not sequential.
- The first 20% of shuffled indices become the test set. The remaining 80% become training data.
Why build this from scratch instead of using sklearn's train_test_split?
The goal of this project is to understand every piece of the gradient boosting pipeline. Writing your own split shows you exactly how data gets shuffled and partitioned. It also keeps the from-scratch portions of the project independent of scikit-learn utilities.
You won't see a change yet. The next chunk calls this function and prints the results.
Load the dataset and run the script
The California Housing dataset contains 20,640 samples with 8 numeric features describing California census block groups. The target variable is the median house value expressed in hundreds of thousands of dollars.
- Below the train_test_split function, add the main block that loads the data and prints summary information:
# ============================================================
# Main
# ============================================================
if __name__ == "__main__":
# Load dataset
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names
print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Features: {feature_names}")
print(f"Target: Median house value (in $100,000s)")
print()
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print(f"Train set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print()What does this code do?
- if __name__ == "__main__": ensures this code only runs when you execute the file directly, not when it's imported as a module.
- fetch_california_housing() downloads the dataset from scikit-learn's bundled datasets. It returns a Bunch object with .data (features) and .target (house values).
- The 80/20 split gives you 16,512 training samples and 4,128 test samples.
- Save gradient_boosting.py.
- Run the script in your terminal by executing this command:
python gradient_boosting.py- You should see output showing 20640 samples, 8 features, the list of feature names, and the train/test split sizes (16512 and 4128).
Not seeing the expected output?
- If you see ModuleNotFoundError, confirm your packages are installed by running pip show numpy scikit-learn matplotlib in your terminal.
- If you see IndentationError, check that the lines inside the if __name__ block are indented with 4 spaces consistently.
- If nothing prints at all, make sure the if __name__ == "__main__": line is at the top level (no indentation) and the code below it is indented.
Help me debug my gradient_boosting.py script that won't print the dataset info.
✔️ Awesome, I've got everything!
Great. Double check you've saved your file before moving on.
ⓧ I'd like to double check the full code
Here is the complete gradient_boosting.py file so far:
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor as SklearnGBR
import matplotlib.pyplot as plt
# ============================================================
# Evaluation Metrics
# ============================================================
def mean_squared_error(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
return 1 - (ss_res / ss_tot)
# ============================================================
# Train/Test Split
# ============================================================
def train_test_split(X, y, test_size=0.2, random_state=42):
rng = np.random.default_rng(random_state)
n_samples = X.shape[0]
indices = rng.permutation(n_samples)
test_count = int(n_samples * test_size)
test_idx = indices[:test_count]
train_idx = indices[test_count:]
return X[train_idx], X[test_idx], y[train_idx], y[test_idx]
# ============================================================
# Main
# ============================================================
if __name__ == "__main__":
# Load dataset
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names
print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Features: {feature_names}")
print(f"Target: Median house value (in $100,000s)")
print()
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print(f"Train set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print()Your data pipeline is ready and your metrics are in place. Next up, you'll build the decision tree that serves as the weak learner inside gradient boosting.
Build the Decision Tree Regressor
In the previous step you loaded the California Housing dataset and set up evaluation metrics. Now you need a model that can actually learn from that data.
Gradient boosting works by combining many weak learners into a strong predictor. The most common weak learner is a shallow decision tree. In this step, you will build a complete decision tree regressor from scratch that will serve as the building block for the full gradient boosting ensemble.
In this step, get ready to:
- Create a Node class and DecisionTreeRegressor class with MSE-based splitting.
- Implement recursive tree building with stopping criteria and feature importance tracking.
- Add prediction via tree traversal and verify the tree produces valid outputs.
Create the Node and DecisionTreeRegressor classes
A decision tree is built from nodes. Each internal node splits data based on a feature and threshold. Each leaf node holds a prediction value (the mean of its samples).
- In gradient_boosting.py, find the comment that says # Main.
- Add the following section comment and Node class above it:
# ============================================================
# Decision Tree Regressor (From Scratch)
# ============================================================
class Node:
def __init__(self, feature_idx=None, threshold=None, left=None, right=None, value=None):
self.feature_idx = feature_idx
self.threshold = threshold
self.left = left
self.right = right
self.value = value
def is_leaf(self):
return self.value is not NoneWhat does this code do?
The Node class represents a single point in the decision tree. Internal nodes store the feature_idx and threshold used to split data left or right.
Leaf nodes store a value (the prediction). The is_leaf() helper tells the traversal algorithm when to stop and return the prediction.
- Save gradient_boosting.py.
- You won't see a visible change yet. The Node class is a data structure that the tree builder will use next.
- Directly below the Node class, add the DecisionTreeRegressor class with its constructor, fit, and _mse methods:
class DecisionTreeRegressor:
def __init__(self, max_depth=3, min_samples_split=2):
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.root = None
self.n_features = None
self.feature_importances_ = None
def fit(self, X, y):
self.n_features = X.shape[1]
self.feature_importances_ = np.zeros(self.n_features)
self.root = self._build_tree(X, y, depth=0)
def _mse(self, y):
if len(y) == 0:
return 0.0
return np.mean((y - np.mean(y)) ** 2)What does this code do?
- max_depth limits how deep the tree can grow. Shallow trees (depth 3) are intentionally weak learners, which is what gradient boosting needs.
- fit() initializes feature importance tracking and kicks off recursive tree building.
- _mse() computes the variance of targets in a node. This is the impurity metric that measures how "mixed" the values are. A pure node (all same values) has MSE of 0.
- Save gradient_boosting.py.
- You won't see a change yet. The next chunk adds the split-finding logic that makes the tree learn.
The heart of a decision tree is finding the best split. For each feature, you need to test candidate thresholds and pick the one that reduces MSE the most.
- Add the _best_split method inside DecisionTreeRegressor, directly below the _mse method:
def _best_split(self, X, y):
n_samples, n_features = X.shape
best_mse_reduction = -np.inf
best_feature = None
best_threshold = None
current_mse = self._mse(y)
for feature_idx in range(n_features):
feature_values = X[:, feature_idx]
percentiles = np.percentile(feature_values, np.linspace(10, 90, 20))
thresholds = np.unique(percentiles)
for threshold in thresholds:
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
n_left = np.sum(left_mask)
n_right = np.sum(right_mask)
if n_left < 1 or n_right < 1:
continue
left_mse = self._mse(y[left_mask])
right_mse = self._mse(y[right_mask])
weighted_mse = (n_left * left_mse + n_right * right_mse) / n_samples
mse_reduction = current_mse - weighted_mse
if mse_reduction > best_mse_reduction:
best_mse_reduction = mse_reduction
best_feature = feature_idx
best_threshold = threshold
return best_feature, best_threshold, best_mse_reductionWhy percentile-based thresholds?
Testing every unique value in a feature column would be extremely slow on 16,000+ samples. Instead, this method tests 20 evenly-spaced percentiles (10th through 90th) per feature.
For each candidate threshold, the method splits data into left and right groups, computes the weighted MSE of both groups, and tracks whichever split produces the largest MSE reduction.
- Save gradient_boosting.py.
- You won't see output yet. The next method ties everything together by recursively calling _best_split to build the full tree.
Implement the recursive tree-building method
The tree grows by recursively finding the best split at each node, then building left and right subtrees. Stopping criteria prevent the tree from growing too deep or splitting nodes with too few samples.
- Add the _build_tree method inside DecisionTreeRegressor, directly below _best_split:
def _build_tree(self, X, y, depth):
n_samples = X.shape[0]
# Stopping criteria: max depth, too few samples, or zero impurity
if (depth >= self.max_depth or
n_samples < self.min_samples_split or
self._mse(y) == 0):
return Node(value=np.mean(y))
feature_idx, threshold, mse_reduction = self._best_split(X, y)
# No valid split found
if feature_idx is None or mse_reduction <= 0:
return Node(value=np.mean(y))
# Track feature importance as cumulative weighted MSE reduction
self.feature_importances_[feature_idx] += n_samples * mse_reduction
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
left_child = self._build_tree(X[left_mask], y[left_mask], depth + 1)
right_child = self._build_tree(X[right_mask], y[right_mask], depth + 1)
return Node(
feature_idx=feature_idx,
threshold=threshold,
left=left_child,
right=right_child
)What does this code do?
- Three stopping criteria create leaf nodes: reaching max_depth, having fewer than min_samples_split samples, or encountering zero impurity (all targets are identical).
- Feature importance is tracked by accumulating n_samples * mse_reduction each time a feature is chosen for a split. Features that reduce error more often and on larger nodes get higher importance.
- The recursion creates left and right children by masking data on the chosen threshold, then calling _build_tree again with incremented depth.
- Save gradient_boosting.py.
- You won't see output yet. The tree can now build itself, but it still needs a predict method to produce outputs.
Implement prediction and verify the tree works
Prediction works by traversing the tree for each sample. At each internal node, the sample's feature value determines whether it goes left or right. When it reaches a leaf, the leaf's value is the prediction.
- Add the predict and _traverse methods inside DecisionTreeRegressor, directly below _build_tree:
def predict(self, X):
return np.array([self._traverse(x, self.root) for x in X])
def _traverse(self, x, node):
if node.is_leaf():
return node.value
if x[node.feature_idx] <= node.threshold:
return self._traverse(x, node.left)
else:
return self._traverse(x, node.right)What does this code do?
- predict() loops over every sample in X and runs each one through the tree individually.
- _traverse() walks a single sample from the root down to a leaf. At each node it checks the sample's feature value against the threshold to decide left or right.
- Save gradient_boosting.py.
Now verify that the tree produces valid predictions on the housing data. You will add temporary test code to the bottom of the main block.
- At the bottom of the if __name__ == "__main__": block, add these lines to test the tree:
# Temporary test - verify decision tree works
tree = DecisionTreeRegressor(max_depth=3, min_samples_split=5)
tree.fit(X_train, y_train)
tree_preds = tree.predict(X_test[:5])
print(f"Tree predictions (first 5): {tree_preds}")
print(f"Predictions are varied: {len(set(tree_preds)) > 1}")
print(f"Feature importance sum: {tree.feature_importances_.sum():.2f}")- Run the script by running:
python gradient_boosting.py- You should see the dataset info from earlier, followed by three new lines: an array of varied prediction values, True confirming predictions are not all identical, and a positive feature importance sum.
- Once verified, remove the temporary test lines (the 5 lines starting from # Temporary test) and save the file.
Seeing all identical predictions or an error?
- If you see a NameError for Node or DecisionTreeRegressor, make sure both classes are defined above the main block, not inside it.
- If predictions are all the same value, check that _best_split returns a tuple of three values and that your indentation places all methods inside the class.
- If feature importance sum is 0.0, verify that _build_tree has the line self.feature_importances_[feature_idx] += n_samples * mse_reduction before the recursive calls.
Why isn't my DecisionTreeRegressor producing varied predictions?
✔️ Awesome, I've got everything!
Great. Double check you removed the temporary test lines and saved your file.
ⓧ I'd like to double check the full code
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor as SklearnGBR
import matplotlib.pyplot as plt
# ============================================================
# Evaluation Metrics
# ============================================================
def mean_squared_error(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
return 1 - (ss_res / ss_tot)
# ============================================================
# Train/Test Split
# ============================================================
def train_test_split(X, y, test_size=0.2, random_state=42):
rng = np.random.default_rng(random_state)
n_samples = X.shape[0]
indices = rng.permutation(n_samples)
test_count = int(n_samples * test_size)
test_idx = indices[:test_count]
train_idx = indices[test_count:]
return X[train_idx], X[test_idx], y[train_idx], y[test_idx]
# ============================================================
# Decision Tree Regressor (From Scratch)
# ============================================================
class Node:
def __init__(self, feature_idx=None, threshold=None, left=None, right=None, value=None):
self.feature_idx = feature_idx
self.threshold = threshold
self.left = left
self.right = right
self.value = value
def is_leaf(self):
return self.value is not None
class DecisionTreeRegressor:
def __init__(self, max_depth=3, min_samples_split=2):
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.root = None
self.n_features = None
self.feature_importances_ = None
def fit(self, X, y):
self.n_features = X.shape[1]
self.feature_importances_ = np.zeros(self.n_features)
self.root = self._build_tree(X, y, depth=0)
def _mse(self, y):
if len(y) == 0:
return 0.0
return np.mean((y - np.mean(y)) ** 2)
def _best_split(self, X, y):
n_samples, n_features = X.shape
best_mse_reduction = -np.inf
best_feature = None
best_threshold = None
current_mse = self._mse(y)
for feature_idx in range(n_features):
feature_values = X[:, feature_idx]
percentiles = np.percentile(feature_values, np.linspace(10, 90, 20))
thresholds = np.unique(percentiles)
for threshold in thresholds:
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
n_left = np.sum(left_mask)
n_right = np.sum(right_mask)
if n_left < 1 or n_right < 1:
continue
left_mse = self._mse(y[left_mask])
right_mse = self._mse(y[right_mask])
weighted_mse = (n_left * left_mse + n_right * right_mse) / n_samples
mse_reduction = current_mse - weighted_mse
if mse_reduction > best_mse_reduction:
best_mse_reduction = mse_reduction
best_feature = feature_idx
best_threshold = threshold
return best_feature, best_threshold, best_mse_reduction
def _build_tree(self, X, y, depth):
n_samples = X.shape[0]
if (depth >= self.max_depth or
n_samples < self.min_samples_split or
self._mse(y) == 0):
return Node(value=np.mean(y))
feature_idx, threshold, mse_reduction = self._best_split(X, y)
if feature_idx is None or mse_reduction <= 0:
return Node(value=np.mean(y))
self.feature_importances_[feature_idx] += n_samples * mse_reduction
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
left_child = self._build_tree(X[left_mask], y[left_mask], depth + 1)
right_child = self._build_tree(X[right_mask], y[right_mask], depth + 1)
return Node(
feature_idx=feature_idx,
threshold=threshold,
left=left_child,
right=right_child
)
def predict(self, X):
return np.array([self._traverse(x, self.root) for x in X])
def _traverse(self, x, node):
if node.is_leaf():
return node.value
if x[node.feature_idx] <= node.threshold:
return self._traverse(x, node.left)
else:
return self._traverse(x, node.right)
# ============================================================
# Main
# ============================================================
if __name__ == "__main__":
# Load dataset
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names
print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Features: {feature_names}")
print(f"Target: Median house value (in $100,000s)")
print()
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print(f"Train set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print()Your decision tree is fully functional. It can split data, grow recursively, and produce predictions. Next up, you will build the gradient boosting algorithm that ensembles many of these shallow trees to create a powerful predictor.
Build the Gradient Boosting Regressor
You have a working decision tree that can fit data and make predictions. A single shallow tree is a weak learner. It captures broad patterns but misses fine-grained details.
Gradient boosting turns that weakness into a strength by chaining many weak trees together. Each tree focuses exclusively on the errors the previous trees made, progressively refining the overall prediction.
In this step, get ready to:
- Implement Friedman's gradient boosting algorithm in a GradientBoostingRegressor class.
- Add predict and feature importance methods to the ensemble.
- Train and evaluate the model on California Housing data.
Create the GradientBoostingRegressor class and implement fit
Friedman's gradient boosting algorithm follows a simple loop. Start with a constant prediction (the mean of all target values). Then repeatedly compute the pseudo-residuals (the gap between actual values and current predictions), fit a new tree to those residuals, and nudge the predictions in the direction the tree suggests.
- In gradient_boosting.py, find the comment line # ============================================================ that sits directly above # Main.
- Add the following class definition above that comment block:
# ============================================================
# Gradient Boosting Regressor (From Scratch)
# ============================================================
class GradientBoostingRegressor:
def __init__(self, n_estimators=50, learning_rate=0.1, max_depth=3,
min_samples_split=5):
self.n_estimators = n_estimators
self.learning_rate = learning_rate
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.trees = []
self.initial_prediction = None
def fit(self, X, y):
self.initial_prediction = np.mean(y)
predictions = np.full(len(y), self.initial_prediction)
self.trees = []
for i in range(self.n_estimators):
residuals = y - predictions
tree = DecisionTreeRegressor(
max_depth=self.max_depth,
min_samples_split=self.min_samples_split
)
tree.fit(X, residuals)
self.trees.append(tree)
update = tree.predict(X)
predictions += self.learning_rate * update
if (i + 1) % 10 == 0:
mse = mean_squared_error(y, predictions)
print(f" Iteration {i + 1}/{self.n_estimators} - Train MSE: {mse:.4f}")What does this code do?
This is the core of gradient boosting. The fit method implements Friedman's algorithm step by step.
- self.initial_prediction = np.mean(y) sets the starting point. Before any trees are added, our best guess for every sample is the average target value.
- residuals = y - predictions computes the pseudo-residuals. For MSE loss, these are simply the errors between actual values and current predictions.
- Each iteration fits a new DecisionTreeRegressor to those residuals, teaching it to predict the remaining errors.
- predictions += self.learning_rate * update applies shrinkage. The learning rate scales down each tree's contribution, preventing overfitting and giving later trees room to correct mistakes.
You won't see a change yet. The next chunk adds the methods needed to generate predictions from the trained ensemble.
Implement predict and feature importances
The ensemble needs two more capabilities. First, a way to generate predictions by replaying all trees in sequence. Second, a way to measure which features contributed most across the entire ensemble.
- Add the following two methods inside the GradientBoostingRegressor class, directly below the fit method:
def predict(self, X):
predictions = np.full(X.shape[0], self.initial_prediction)
for tree in self.trees:
predictions += self.learning_rate * tree.predict(X)
return predictions
def feature_importances(self):
total_importance = np.zeros(self.trees[0].n_features)
for tree in self.trees:
total_importance += tree.feature_importances_
total = total_importance.sum()
if total > 0:
total_importance /= total
return total_importanceWhat does this code do?
- predict reconstructs the full ensemble prediction. It starts from the initial mean, then adds each tree's contribution scaled by the learning rate. This mirrors exactly what fit accumulated during training.
- feature_importances aggregates the MSE-reduction scores from every tree in the ensemble. It then normalizes so the importances sum to 1, giving a relative ranking of how much each feature contributed to reducing prediction error.
You won't see output yet. The next chunk adds the training and evaluation code that exercises these methods.
Train and evaluate the model
The class is complete. Now you need to instantiate it, train on the housing data, and measure how well the ensemble performs on both the training and test sets.
- In the if __name__ == "__main__": block, find the last two lines:
print(f"Test set: {X_test.shape[0]} samples")
print()- Add the following training, evaluation, and feature importance code directly below those lines:
# Train from-scratch model
print("Training Gradient Boosting Regressor (from scratch)...")
gbr = GradientBoostingRegressor(
n_estimators=50,
learning_rate=0.1,
max_depth=3,
min_samples_split=5
)
gbr.fit(X_train, y_train)
print()
# Evaluate from-scratch model
y_pred_train = gbr.predict(X_train)
y_pred_test = gbr.predict(X_test)
print("From-Scratch Results:")
print(f" Train MSE: {mean_squared_error(y_train, y_pred_train):.4f}")
print(f" Test MSE: {mean_squared_error(y_test, y_pred_test):.4f}")
print(f" Train R^2: {r_squared(y_train, y_pred_train):.4f}")
print(f" Test R^2: {r_squared(y_test, y_pred_test):.4f}")
print()
# Feature importances
importances = gbr.feature_importances()
print("Feature Importances (from scratch):")
for name, imp in sorted(zip(feature_names, importances), key=lambda x: x[1], reverse=True):
print(f" {name}: {imp:.4f}")
print()What does this code do?
- The model is instantiated with 50 trees, a learning rate of 0.1, and max depth 3. These are conservative hyperparameters that balance accuracy with training speed.
- gbr.fit(X_train, y_train) runs the full boosting loop. You will see progress printed every 10 iterations showing the training MSE decreasing.
- Evaluating on both train and test sets reveals the generalization gap. Train R-squared should be higher than test R-squared, confirming the model fits training data more tightly than unseen data.
- Feature importances are sorted in descending order so you can see which housing features (like median income) drive predictions the most.
- Save gradient_boosting.py.
- Run the script by executing:
python gradient_boosting.py- You should see training progress printed every 10 iterations with decreasing MSE values, followed by the final results showing Train MSE, Test MSE, Train R-squared, and Test R-squared. The test R-squared should be above 0.60.
What should I see?
The training will take 1-3 minutes depending on your hardware. You will see MSE values printed at iterations 10, 20, 30, 40, and 50. The MSE should decrease from around 0.5-0.6 down to roughly 0.3-0.4 as more trees are added.
The final results should show a test R-squared above 0.60 and a train R-squared that is noticeably higher (perhaps 0.70-0.80). The top feature importance should be MedInc (median income), since income strongly predicts housing prices.
Seeing an error or unexpected output?
- If you see an IndentationError, make sure the predict and feature_importances methods are indented at the same level as fit (one level inside the class).
- If R-squared is negative or near zero, check that your fit method uses residuals = y - predictions (not the other way around) and that predictions += self.learning_rate * update adds (not subtracts).
- If you see AttributeError: 'GradientBoostingRegressor' has no attribute 'initial_prediction', make sure self.initial_prediction = None is in your __init__ method.
Why is my gradient boosting model producing bad predictions?
✔️ Awesome, I've got everything!
Great. Double check you have saved your file before moving on.
ⓧ I'd like to double check the full code
Here is the complete gradient_boosting.py file as it should look at the end of this step:
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor as SklearnGBR
import matplotlib.pyplot as plt
# ============================================================
# Evaluation Metrics
# ============================================================
def mean_squared_error(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
return 1 - (ss_res / ss_tot)
# ============================================================
# Train/Test Split
# ============================================================
def train_test_split(X, y, test_size=0.2, random_state=42):
rng = np.random.default_rng(random_state)
n_samples = X.shape[0]
indices = rng.permutation(n_samples)
test_count = int(n_samples * test_size)
test_idx = indices[:test_count]
train_idx = indices[test_count:]
return X[train_idx], X[test_idx], y[train_idx], y[test_idx]
# ============================================================
# Decision Tree Regressor (From Scratch)
# ============================================================
class Node:
def __init__(self, feature_idx=None, threshold=None, left=None, right=None, value=None):
self.feature_idx = feature_idx
self.threshold = threshold
self.left = left
self.right = right
self.value = value
def is_leaf(self):
return self.value is not None
class DecisionTreeRegressor:
def __init__(self, max_depth=3, min_samples_split=2):
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.root = None
self.n_features = None
self.feature_importances_ = None
def fit(self, X, y):
self.n_features = X.shape[1]
self.feature_importances_ = np.zeros(self.n_features)
self.root = self._build_tree(X, y, depth=0)
def _mse(self, y):
if len(y) == 0:
return 0.0
return np.mean((y - np.mean(y)) ** 2)
def _best_split(self, X, y):
n_samples, n_features = X.shape
best_mse_reduction = -np.inf
best_feature = None
best_threshold = None
current_mse = self._mse(y)
for feature_idx in range(n_features):
feature_values = X[:, feature_idx]
percentiles = np.percentile(feature_values, np.linspace(10, 90, 20))
thresholds = np.unique(percentiles)
for threshold in thresholds:
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
n_left = np.sum(left_mask)
n_right = np.sum(right_mask)
if n_left < 1 or n_right < 1:
continue
left_mse = self._mse(y[left_mask])
right_mse = self._mse(y[right_mask])
weighted_mse = (n_left * left_mse + n_right * right_mse) / n_samples
mse_reduction = current_mse - weighted_mse
if mse_reduction > best_mse_reduction:
best_mse_reduction = mse_reduction
best_feature = feature_idx
best_threshold = threshold
return best_feature, best_threshold, best_mse_reduction
def _build_tree(self, X, y, depth):
n_samples = X.shape[0]
if (depth >= self.max_depth or
n_samples < self.min_samples_split or
self._mse(y) == 0):
return Node(value=np.mean(y))
feature_idx, threshold, mse_reduction = self._best_split(X, y)
if feature_idx is None or mse_reduction <= 0:
return Node(value=np.mean(y))
self.feature_importances_[feature_idx] += n_samples * mse_reduction
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
left_child = self._build_tree(X[left_mask], y[left_mask], depth + 1)
right_child = self._build_tree(X[right_mask], y[right_mask], depth + 1)
return Node(
feature_idx=feature_idx,
threshold=threshold,
left=left_child,
right=right_child
)
def predict(self, X):
return np.array([self._traverse(x, self.root) for x in X])
def _traverse(self, x, node):
if node.is_leaf():
return node.value
if x[node.feature_idx] <= node.threshold:
return self._traverse(x, node.left)
else:
return self._traverse(x, node.right)
# ============================================================
# Gradient Boosting Regressor (From Scratch)
# ============================================================
class GradientBoostingRegressor:
def __init__(self, n_estimators=50, learning_rate=0.1, max_depth=3,
min_samples_split=5):
self.n_estimators = n_estimators
self.learning_rate = learning_rate
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.trees = []
self.initial_prediction = None
def fit(self, X, y):
self.initial_prediction = np.mean(y)
predictions = np.full(len(y), self.initial_prediction)
self.trees = []
for i in range(self.n_estimators):
residuals = y - predictions
tree = DecisionTreeRegressor(
max_depth=self.max_depth,
min_samples_split=self.min_samples_split
)
tree.fit(X, residuals)
self.trees.append(tree)
update = tree.predict(X)
predictions += self.learning_rate * update
if (i + 1) % 10 == 0:
mse = mean_squared_error(y, predictions)
print(f" Iteration {i + 1}/{self.n_estimators} - Train MSE: {mse:.4f}")
def predict(self, X):
predictions = np.full(X.shape[0], self.initial_prediction)
for tree in self.trees:
predictions += self.learning_rate * tree.predict(X)
return predictions
def feature_importances(self):
total_importance = np.zeros(self.trees[0].n_features)
for tree in self.trees:
total_importance += tree.feature_importances_
total = total_importance.sum()
if total > 0:
total_importance /= total
return total_importance
# ============================================================
# Main
# ============================================================
if __name__ == "__main__":
# Load dataset
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names
print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Features: {feature_names}")
print(f"Target: Median house value (in $100,000s)")
print()
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print(f"Train set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print()
# Train from-scratch model
print("Training Gradient Boosting Regressor (from scratch)...")
gbr = GradientBoostingRegressor(
n_estimators=50,
learning_rate=0.1,
max_depth=3,
min_samples_split=5
)
gbr.fit(X_train, y_train)
print()
# Evaluate from-scratch model
y_pred_train = gbr.predict(X_train)
y_pred_test = gbr.predict(X_test)
print("From-Scratch Results:")
print(f" Train MSE: {mean_squared_error(y_train, y_pred_train):.4f}")
print(f" Test MSE: {mean_squared_error(y_test, y_pred_test):.4f}")
print(f" Train R^2: {r_squared(y_train, y_pred_train):.4f}")
print(f" Test R^2: {r_squared(y_test, y_pred_test):.4f}")
print()
# Feature importances
importances = gbr.feature_importances()
print("Feature Importances (from scratch):")
for name, imp in sorted(zip(feature_names, importances), key=lambda x: x[1], reverse=True):
print(f" {name}: {imp:.4f}")
print()Your from-scratch gradient boosting model is producing real, meaningful predictions on housing data. Next, you will validate these results by comparing them side-by-side against scikit-learn's production-grade implementation.
Validate Against scikit-learn
Your from-scratch gradient boosting model is trained and producing predictions. The real question now is: does it actually work correctly?
The best way to validate a from-scratch implementation is to compare it against a battle-tested library. In this step, you will train scikit-learn's GradientBoostingRegressor with identical hyperparameters and compare the results side-by-side.
In this step, get ready to:
- Train sklearn's GradientBoostingRegressor and compare test metrics.
- Compare feature importance rankings from both models.
- Visualize feature importances and predicted vs. actual scatter plots.
Train sklearn's model and compare metrics
If your implementation follows Friedman's algorithm correctly, it should produce R-squared and MSE values close to sklearn's output when trained with the same hyperparameters. The scores won't be identical because sklearn uses exact threshold search and per-leaf optimization, but they should be within approximately 0.05-0.10 R-squared of each other.
- In gradient_boosting.py, find the line that says print() at the end of the feature importances loop (the last line in your file).
- Add the sklearn comparison code below it:
# Compare with sklearn
print("Training sklearn GradientBoostingRegressor for comparison...")
sklearn_gbr = SklearnGBR(
n_estimators=50,
learning_rate=0.1,
max_depth=3,
min_samples_split=5,
random_state=42
)
sklearn_gbr.fit(X_train, y_train)
sklearn_pred_test = sklearn_gbr.predict(X_test)
print(f" sklearn Test MSE: {mean_squared_error(y_test, sklearn_pred_test):.4f}")
print(f" sklearn Test R^2: {r_squared(y_test, sklearn_pred_test):.4f}")
print()What does this code do?
You create an sklearn GradientBoostingRegressor with the exact same hyperparameters your from-scratch model used. The random_state=42 ensures sklearn's internal randomness is reproducible.
By feeding both models the same training data and evaluating on the same test set, any difference in output comes from implementation details (threshold search strategy, leaf optimization) rather than data differences.
- Save gradient_boosting.py and run the script.
- You should see sklearn's test MSE and R-squared printed below your from-scratch results. Both R-squared values should be above 0.60, and the difference between them should be within approximately 0.05-0.10.
sklearn scores are much higher than yours?
A gap larger than 0.10 in R-squared is expected. Your from-scratch model uses percentile-based threshold candidates (20 per feature), while sklearn evaluates every unique value. This means sklearn finds slightly better splits at every node, compounding across 50 trees.
If the gap is extremely large (above 0.20), double-check that your GradientBoostingRegressor uses min_samples_split=5 (matching sklearn's parameter) and that the learning rate is 0.1.
Help me debug why my from-scratch gradient boosting R-squared is much lower than sklearn's
Compare feature importance rankings
Matching R-squared scores is one signal. Matching feature importance rankings is another. If both models agree that the same features drive predictions, your tree-splitting logic is working correctly.
- Add the sklearn feature importances display below the comparison code you just added:
# sklearn feature importances
print("Feature Importances (sklearn):")
for name, imp in sorted(zip(feature_names, sklearn_gbr.feature_importances_), key=lambda x: x[1], reverse=True):
print(f" {name}: {imp:.4f}")
print()What does this code do?
sklearn exposes feature importances through the feature_importances_ attribute (based on impurity reduction, just like your implementation). This loop sorts features by importance and prints them in descending order so you can compare the rankings against your from-scratch output above.
- Save gradient_boosting.py and run the script.
- Compare the two feature importance lists in the output. MedInc (median income) should dominate both rankings as the most important feature.
Why does MedInc dominate?
The California Housing dataset predicts median house values. Median income of a neighborhood is the single strongest predictor of house prices, so both models consistently pick it for the highest-impact splits across all 50 boosting iterations.
Visualize feature importances and predictions
Numbers in a terminal are useful, but visual comparisons make it immediately obvious whether the two models agree. You will create two figures: a feature importance bar chart and a predicted-vs-actual scatter plot.
- Add the feature importance visualization below the sklearn feature importances loop:
# Visualization: Feature Importance Comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sorted_idx = np.argsort(importances)
axes[0].barh(range(len(feature_names)), importances[sorted_idx])
axes[0].set_yticks(range(len(feature_names)))
axes[0].set_yticklabels(np.array(feature_names)[sorted_idx])
axes[0].set_xlabel("Importance")
axes[0].set_title("Feature Importances (From Scratch)")
sklearn_importances = sklearn_gbr.feature_importances_
sorted_idx_sk = np.argsort(sklearn_importances)
axes[1].barh(range(len(feature_names)), sklearn_importances[sorted_idx_sk])
axes[1].set_yticks(range(len(feature_names)))
axes[1].set_yticklabels(np.array(feature_names)[sorted_idx_sk])
axes[1].set_xlabel("Importance")
axes[1].set_title("Feature Importances (sklearn)")
plt.tight_layout()
plt.savefig("feature_importances.png", dpi=100)
plt.show()What does this code do?
This creates a side-by-side horizontal bar chart. Each subplot sorts features by importance (lowest to highest, bottom to top) so the most important feature appears at the top of each chart.
The figure is saved as feature_importances.png and also displayed interactively. If both charts show the same ranking pattern, your splitting logic matches sklearn's behavior.
- Now add the predicted vs. actual scatter plot below the feature importance visualization:
# Visualization: Predicted vs Actual
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(y_test, y_pred_test, alpha=0.3, s=10)
axes[0].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--", lw=2)
axes[0].set_xlabel("Actual")
axes[0].set_ylabel("Predicted")
axes[0].set_title("From Scratch: Predicted vs Actual")
axes[1].scatter(y_test, sklearn_pred_test, alpha=0.3, s=10)
axes[1].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--", lw=2)
axes[1].set_xlabel("Actual")
axes[1].set_ylabel("Predicted")
axes[1].set_title("sklearn: Predicted vs Actual")
plt.tight_layout()
plt.savefig("predicted_vs_actual.png", dpi=100)
plt.show()What does this code do?
Each scatter plot shows actual house values (x-axis) vs. predicted values (y-axis) for the test set. The red dashed line represents perfect predictions. Points clustered tightly around this line indicate a well-performing model.
Both plots should show a similar spread pattern. If your from-scratch scatter looks significantly more dispersed than sklearn's, the difference in R-squared will confirm that visually.
- Save gradient_boosting.py and run the script.
- Two figure windows should appear sequentially. The first shows the feature importance bar charts. The second shows the predicted vs. actual scatter plots. Close each window to proceed to the next.
- Confirm that feature_importances.png and predicted_vs_actual.png were saved in your project folder.
Plots not appearing?
If you are running in a headless environment (remote server, WSL without display), plt.show() may not display a window. The files are still saved. Check your project folder for the .png files and open them directly.
If you see an indentation error, make sure all the visualization code is indented with 4 spaces (inside the if __name__ == "__main__": block).
Help me debug why my matplotlib plots are not showing or saving
✔️ Awesome, I've got everything!
Great. Double check you have saved your file and both PNG images appear in your project folder.
ⓧ I'd like to double check the full code
Here is the complete gradient_boosting.py at this point:
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import GradientBoostingRegressor as SklearnGBR
import matplotlib.pyplot as plt
# ============================================================
# Evaluation Metrics
# ============================================================
def mean_squared_error(y_true, y_pred):
return np.mean((y_true - y_pred) ** 2)
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
return 1 - (ss_res / ss_tot)
# ============================================================
# Train/Test Split
# ============================================================
def train_test_split(X, y, test_size=0.2, random_state=42):
rng = np.random.default_rng(random_state)
n_samples = X.shape[0]
indices = rng.permutation(n_samples)
test_count = int(n_samples * test_size)
test_idx = indices[:test_count]
train_idx = indices[test_count:]
return X[train_idx], X[test_idx], y[train_idx], y[test_idx]
# ============================================================
# Decision Tree Regressor (From Scratch)
# ============================================================
class Node:
def __init__(self, feature_idx=None, threshold=None, left=None, right=None, value=None):
self.feature_idx = feature_idx
self.threshold = threshold
self.left = left
self.right = right
self.value = value
def is_leaf(self):
return self.value is not None
class DecisionTreeRegressor:
def __init__(self, max_depth=3, min_samples_split=2):
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.root = None
self.n_features = None
self.feature_importances_ = None
def fit(self, X, y):
self.n_features = X.shape[1]
self.feature_importances_ = np.zeros(self.n_features)
self.root = self._build_tree(X, y, depth=0)
def _mse(self, y):
if len(y) == 0:
return 0.0
return np.mean((y - np.mean(y)) ** 2)
def _best_split(self, X, y):
n_samples, n_features = X.shape
best_mse_reduction = -np.inf
best_feature = None
best_threshold = None
current_mse = self._mse(y)
for feature_idx in range(n_features):
feature_values = X[:, feature_idx]
percentiles = np.percentile(feature_values, np.linspace(10, 90, 20))
thresholds = np.unique(percentiles)
for threshold in thresholds:
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
n_left = np.sum(left_mask)
n_right = np.sum(right_mask)
if n_left < 1 or n_right < 1:
continue
left_mse = self._mse(y[left_mask])
right_mse = self._mse(y[right_mask])
weighted_mse = (n_left * left_mse + n_right * right_mse) / n_samples
mse_reduction = current_mse - weighted_mse
if mse_reduction > best_mse_reduction:
best_mse_reduction = mse_reduction
best_feature = feature_idx
best_threshold = threshold
return best_feature, best_threshold, best_mse_reduction
def _build_tree(self, X, y, depth):
n_samples = X.shape[0]
if (depth >= self.max_depth or
n_samples < self.min_samples_split or
self._mse(y) == 0):
return Node(value=np.mean(y))
feature_idx, threshold, mse_reduction = self._best_split(X, y)
if feature_idx is None or mse_reduction <= 0:
return Node(value=np.mean(y))
self.feature_importances_[feature_idx] += n_samples * mse_reduction
left_mask = X[:, feature_idx] <= threshold
right_mask = ~left_mask
left_child = self._build_tree(X[left_mask], y[left_mask], depth + 1)
right_child = self._build_tree(X[right_mask], y[right_mask], depth + 1)
return Node(
feature_idx=feature_idx,
threshold=threshold,
left=left_child,
right=right_child
)
def predict(self, X):
return np.array([self._traverse(x, self.root) for x in X])
def _traverse(self, x, node):
if node.is_leaf():
return node.value
if x[node.feature_idx] <= node.threshold:
return self._traverse(x, node.left)
else:
return self._traverse(x, node.right)
# ============================================================
# Gradient Boosting Regressor (From Scratch)
# ============================================================
class GradientBoostingRegressor:
def __init__(self, n_estimators=50, learning_rate=0.1, max_depth=3,
min_samples_split=5):
self.n_estimators = n_estimators
self.learning_rate = learning_rate
self.max_depth = max_depth
self.min_samples_split = min_samples_split
self.trees = []
self.initial_prediction = None
def fit(self, X, y):
self.initial_prediction = np.mean(y)
predictions = np.full(len(y), self.initial_prediction)
self.trees = []
for i in range(self.n_estimators):
residuals = y - predictions
tree = DecisionTreeRegressor(
max_depth=self.max_depth,
min_samples_split=self.min_samples_split
)
tree.fit(X, residuals)
self.trees.append(tree)
update = tree.predict(X)
predictions += self.learning_rate * update
if (i + 1) % 10 == 0:
mse = mean_squared_error(y, predictions)
print(f" Iteration {i + 1}/{self.n_estimators} - Train MSE: {mse:.4f}")
def predict(self, X):
predictions = np.full(X.shape[0], self.initial_prediction)
for tree in self.trees:
predictions += self.learning_rate * tree.predict(X)
return predictions
def feature_importances(self):
total_importance = np.zeros(self.trees[0].n_features)
for tree in self.trees:
total_importance += tree.feature_importances_
total = total_importance.sum()
if total > 0:
total_importance /= total
return total_importance
# ============================================================
# Main
# ============================================================
if __name__ == "__main__":
# Load dataset
data = fetch_california_housing()
X, y = data.data, data.target
feature_names = data.feature_names
print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Features: {feature_names}")
print(f"Target: Median house value (in $100,000s)")
print()
# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print(f"Train set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print()
# Train from-scratch model
print("Training Gradient Boosting Regressor (from scratch)...")
gbr = GradientBoostingRegressor(
n_estimators=50,
learning_rate=0.1,
max_depth=3,
min_samples_split=5
)
gbr.fit(X_train, y_train)
print()
# Evaluate from-scratch model
y_pred_train = gbr.predict(X_train)
y_pred_test = gbr.predict(X_test)
print("From-Scratch Results:")
print(f" Train MSE: {mean_squared_error(y_train, y_pred_train):.4f}")
print(f" Test MSE: {mean_squared_error(y_test, y_pred_test):.4f}")
print(f" Train R^2: {r_squared(y_train, y_pred_train):.4f}")
print(f" Test R^2: {r_squared(y_test, y_pred_test):.4f}")
print()
# Feature importances
importances = gbr.feature_importances()
print("Feature Importances (from scratch):")
for name, imp in sorted(zip(feature_names, importances), key=lambda x: x[1], reverse=True):
print(f" {name}: {imp:.4f}")
print()
# Compare with sklearn
print("Training sklearn GradientBoostingRegressor for comparison...")
sklearn_gbr = SklearnGBR(
n_estimators=50,
learning_rate=0.1,
max_depth=3,
min_samples_split=5,
random_state=42
)
sklearn_gbr.fit(X_train, y_train)
sklearn_pred_test = sklearn_gbr.predict(X_test)
print(f" sklearn Test MSE: {mean_squared_error(y_test, sklearn_pred_test):.4f}")
print(f" sklearn Test R^2: {r_squared(y_test, sklearn_pred_test):.4f}")
print()
# sklearn feature importances
print("Feature Importances (sklearn):")
for name, imp in sorted(zip(feature_names, sklearn_gbr.feature_importances_), key=lambda x: x[1], reverse=True):
print(f" {name}: {imp:.4f}")
print()
# Visualization: Feature Importance Comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sorted_idx = np.argsort(importances)
axes[0].barh(range(len(feature_names)), importances[sorted_idx])
axes[0].set_yticks(range(len(feature_names)))
axes[0].set_yticklabels(np.array(feature_names)[sorted_idx])
axes[0].set_xlabel("Importance")
axes[0].set_title("Feature Importances (From Scratch)")
sklearn_importances = sklearn_gbr.feature_importances_
sorted_idx_sk = np.argsort(sklearn_importances)
axes[1].barh(range(len(feature_names)), sklearn_importances[sorted_idx_sk])
axes[1].set_yticks(range(len(feature_names)))
axes[1].set_yticklabels(np.array(feature_names)[sorted_idx_sk])
axes[1].set_xlabel("Importance")
axes[1].set_title("Feature Importances (sklearn)")
plt.tight_layout()
plt.savefig("feature_importances.png", dpi=100)
plt.show()
# Visualization: Predicted vs Actual
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(y_test, y_pred_test, alpha=0.3, s=10)
axes[0].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--", lw=2)
axes[0].set_xlabel("Actual")
axes[0].set_ylabel("Predicted")
axes[0].set_title("From Scratch: Predicted vs Actual")
axes[1].scatter(y_test, sklearn_pred_test, alpha=0.3, s=10)
axes[1].plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--", lw=2)
axes[1].set_xlabel("Actual")
axes[1].set_ylabel("Predicted")
axes[1].set_title("sklearn: Predicted vs Actual")
plt.tight_layout()
plt.savefig("predicted_vs_actual.png", dpi=100)
plt.show()Your from-scratch implementation is validated. Both models agree on which features matter most and produce competitive prediction accuracy. Next up, you will push this further with a Secret Mission that adds early stopping to prevent overfitting.
Secret mission
Implement Early Stopping
Your gradient boosting model trains for a fixed number of iterations, but what if it starts overfitting halfway through? In production, models need to know when to stop. Implement early stopping with a validation set and plot learning curves to visualize exactly when your model begins to overfit.
Clean Up Your Resources
Clean Up Your Resources
This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files for future reference or delete them entirely.
Resources you used:
- gradient_boosting.py (your from-scratch implementation).
- feature_importances.png (feature importance comparison chart).
- predicted_vs_actual.png (predicted vs actual scatter plots).
- learning_curve.png (early stopping learning curve from the Secret Mission).
Keep everything running
No action needed. Choose this if you plan to keep experimenting with your gradient boosting implementation or want to reference the code later.
- Your files stay where they are. No processes are running and nothing is consuming resources.
Pause - I'll come back to this later
Since this is a local project with no running processes, there is nothing to pause. Your files will remain on disk until you delete them.
- No action required. Come back anytime and run the script again with the same results.
Delete - I don't want to use this again
Remove all project files to start fresh.
- Delete the Python file and generated images by running this command from the directory where you created them:
- On macOS/Linux:
- rm gradient_boosting.py feature_importances.png predicted_vs_actual.png learning_curve.png
- On Windows:
- del gradient_boosting.py feature_importances.png predicted_vs_actual.png learning_curve.png
That's a Wrap!
That's a Wrap!
Nice work! You just built a complete Gradient Boosting Regressor from the ground up using nothing but Python and NumPy.
You've proven that you understand the full algorithm by implementing every piece yourself and validating it against a production-grade library on real housing data.
What you learned:
- Built a complete Decision Tree Regressor from scratch using NumPy with recursive MSE-based splitting and percentile-based candidate thresholds.
- Implemented Friedman's Gradient Boosting algorithm: initialization with the target mean, pseudo-residual computation, sequential tree fitting, and shrinkage-controlled prediction.
- Trained on 20,640 real housing samples and achieved competitive R-squared against scikit-learn's optimized implementation.
- Validated understanding by comparing feature importance rankings and prediction accuracy against a production-grade library.
- Secret Mission: Added early stopping with validation monitoring and visualized learning curves to demonstrate overfitting detection.
Ready to quiz yourself?