Ordinary Least Squares (OLS) › Implementation
Machine Learning from ScratchAll chapters Ordinary Least Squares (OLS) › ImplementationOne section · Sections in this chapter5 / 8 ↓
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Let's implement OLS regression from scratch using NumPy to understand the underlying mathematics, then compare it with scikit-learn's implementation for practical applications.
Step 1: Generate Sample Data
First, we'll create a synthetic dataset with known relationships to test our implementation:
import numpy as np
# Set random seed for reproducibility
np.random.seed(42)
# Create sample data
n_samples = 100
X = np.random.randn(n_samples, 2) # 2 features
true_coef = np.array([2.5, -1.8, 3.2]) # [intercept, coef1, coef2]
y = true_coef[0] + X @ true_coef[1:] + np.random.normal(0, 0.5, n_samples)Step 2: Build the Design Matrix
The design matrix needs to include a column of ones for the intercept term:
# Add intercept column to design matrix
X_design = np.column_stack([np.ones(n_samples), X])
print(f"Design matrix shape: {X_design.shape}")
print(f"First few rows of design matrix:\n{X_design[:5]}")Design matrix shape: (100, 3) First few rows of design matrix: [[ 1. 0.49671415 -0.1382643 ] [ 1. 0.64768854 1.52302986] [ 1. -0.23415337 -0.23413696] [ 1. 1.57921282 0.76743473] [ 1. -0.46947439 0.54256004]]
The design matrix has 100 rows (observations) and 3 columns (intercept + 2 features). The first column is all ones, which allows us to estimate the intercept term in our linear model.
Step 3: Compute OLS Coefficients
Now we'll implement the OLS closed-form solution using matrix operations:
# OLS closed-form solution: β̂ = (X^T X)^(-1) X^T y
beta_hat = np.linalg.inv(X_design.T @ X_design) @ X_design.T @ y
print(f"Estimated coefficients: {beta_hat}")
print(f"True coefficients: {true_coef}")Estimated coefficients: [ 2.54639669 -1.70463983 3.11392895] True coefficients: [ 2.5 -1.8 3.2]
The estimated coefficients should be close to the true values, with small differences due to the random noise we added. This demonstrates that OLS is working correctly by recovering the underlying linear relationship despite the presence of noise.
Step 4: Make Predictions and Calculate Performance
Let's generate predictions and evaluate the model's performance:
# Make predictions
y_pred = X_design @ beta_hat
# Calculate performance metrics
mse = np.mean((y - y_pred) ** 2)
r_squared = 1 - np.sum((y - y_pred) ** 2) / np.sum((y - np.mean(y)) ** 2)
print(f"Mean Squared Error: {mse:.4f}")
print(f"R-squared: {r_squared:.4f}")
print(
f"\nFitted equation: y = {beta_hat[0]:.3f} + {beta_hat[1]:.3f}x₁ + {beta_hat[2]:.3f}x₂"
)Mean Squared Error: 0.2775 R-squared: 0.9762 Fitted equation: y = 2.546 + -1.705x₁ + 3.114x₂
The model should show good performance with a high R-squared value, meaning it explains a significant portion of the variance in the target variable. The MSE indicates the average squared prediction error, and the fitted equation shows the estimated effect of each feature on the target variable.
Step 5: Scikit-learn Implementation
For comparison, let's implement the same model using scikit-learn:
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
# Fit OLS using scikit-learn
model = LinearRegression()
model.fit(X, y)
# Make predictions
y_pred_sklearn = model.predict(X)
# Calculate metrics
mse_sklearn = mean_squared_error(y, y_pred_sklearn)
r2_sklearn = r2_score(y, y_pred_sklearn)Scikit-learn Results: Intercept: 2.546 Coefficients: [-1.70463983 3.11392895] Mean Squared Error: 0.2775 R-squared: 0.9762 Comparison with NumPy implementation: Intercept difference: 0.000000 Coefficient differences: [4.4408921e-16 4.4408921e-16]
The scikit-learn implementation produces identical results to our NumPy implementation (differences are essentially zero due to numerical precision). This confirms that both methods solve the same mathematical problem, with scikit-learn providing a reliable, feature-rich interface for practical applications.
Step 6: Visualizing the Results
Let's create a visualization to see how well our OLS model fits the data:
The scatter plots show how well the OLS model captures the relationships between each feature and the target variable. The close alignment between actual and predicted values (shown in blue and red respectively) confirms the model's effectiveness.
Key Parameters
Below are some of the main parameters that affect how the OLS model works and performs.
fit_intercept: Whether to calculate the intercept for this model (default: True). Set to False if you want to force the line through the origin.normalize: Whether to normalize the regressors before regression (default: False). Normalization can help with numerical stability and interpretation.copy_X: Whether to copy X or overwrite the original (default: True). Set to False to save memory for large datasets.n_jobs: Number of jobs to use for computation (default: None). Useful for parallel processing with large datasets.
Key Methods
The following are the most commonly used methods for interacting with the OLS model.
fit(X, y): Fits the linear model to the training data. X should be a 2D array of features, y should be a 1D array of target values.predict(X): Predicts target values for new data using the fitted model. Returns predicted values as a 1D array.score(X, y): Returns the coefficient of determination (R²) for the given test data and labels. Higher values indicate better fit.get_params(): Returns the parameters of the estimator. Useful for model inspection and hyperparameter tuning.