Linear regression is a statistical method that models the relationship between a dependent variable and one or more independent variables by fitting a linear equation to the observed data. In its simplest form, it helps us understand how one variable changes when another is modified.

While there are many Python packages like Scikit-Learn that offer functions and methods to perform linear regression, here we will implement it from scratch using NumPy. Let’s explore how to do this, breaking down the process into clear, manageable steps.

Setting Up the Environment

First, let’s import the necessary libraries and create some sample data:

import numpy as np
import matplotlib.pyplot as plt

# Set random seed for reproducibility
np.random.seed(42)

# Generate sample data
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1) * 0.1

In our setup, we are creating synthetic data where we know the true relationship: y = 4 + 3x with some added noise. The intercept is 4, and the slope is 3. Our goal is to see how well our linear regression methods can recover these values.

Implementing Linear Regression Using Normal Equation

The normal equation provides a direct mathematical solution using linear algebra. It’s derived from minimizing the sum of squared residuals and gives us θ = (X^T X)^(-1) X^T y, where θ contains our regression coefficients. While this might look intimidating, NumPy makes it straightforward to implement:

# Add bias term to X
X_b = np.c_[np.ones((100, 1)), X]

# Calculate parameters using normal equation
theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

print(f"Intercept: {theta[0][0]:.2f}")
print(f"Slope: {theta[1][0]:.2f}")

Output:

Intercept: 4.02
Slope: 2.98

These results are remarkably close to our true values (4 and 3)! The normal equation method found these values in one step, with no iterations needed. Let’s visualize the results:

# Plotting the results
plt.scatter(X, y, label='Data Points')
plt.plot(X, X_b.dot(theta), color='red', label='Regression Line')
plt.xlabel('X')
plt.ylabel('y')
plt.legend(fontsize=14)
plt.title('Linear Regression using Normal Equation', fontsize=16)
plt.show()

The plot shows how well our regression line (in red) fits the data points (in blue). The tight clustering of points around the line indicates a strong linear relationship.

Implementing Gradient Descent

While the normal equation is elegant, it can be computationally expensive for large datasets because it involves matrix inversion. Gradient descent offers an iterative alternative that can be more efficient for big data. It works by taking small steps in the direction that reduces the error the most:

# Initialize parameters
theta = np.random.randn(2, 1)
learning_rate = 0.01
n_iterations = 1000

for iteration in range(n_iterations):
    # Compute predictions
    y_pred = X_b.dot(theta)
    
    # Compute gradients
    gradients = 2/100 * X_b.T.dot(y_pred - y)
    
    # Update parameters
    theta = theta - learning_rate * gradients

print(f"Intercept: {theta[0][0]:.2f}")
print(f"Slope: {theta[1][0]:.2f}")

Output:

Intercept: 3.97
Slope: 3.02

Notice how gradient descent arrived at very similar values to the normal equation method! The small differences are due to the iterative nature of the algorithm and our chosen learning rate.

Making Predictions with the Model

Once we have our trained model, we can use it to make predictions:

# Create test data
X_test = np.array([[1.5], [3.0]])
X_test_b = np.c_[np.ones((2, 1)), X_test]

# Make predictions
predictions = X_test_b.dot(theta)

print("\nPredictions for X values:")
for x, pred in zip(X_test, predictions):
    print(f"X = {x[0]:.1f}, Predicted y = {pred[0]:.2f}")

# Evaluate model performance
y_pred = X_b.dot(theta)
mse = np.mean((y - y_pred) ** 2)
r2 = 1 - (np.sum((y - y_pred) ** 2) / np.sum((y - np.mean(y)) ** 2))

print(f"\nModel Performance:")
print(f"Mean Squared Error: {mse:.4f}")
print(f"R-squared Score: {r2:.4f}")

Output:

Predictions for X values:
X = 1.5, Predicted y = 8.50
X = 3.0, Predicted y = 13.03

Model Performance:
Mean Squared Error: 0.0088
R-squared Score: 0.9972

Our model performs exceptionally well! The very low Mean Squared Error (0.0088) and high R-squared score (0.9972) indicate that our model explains about 99.72% of the variance in the data. This is expected since we generated clean data with only a small amount of noise.

Conclusions

NumPy provides foundational tools for implementing linear regression from scratch. We’ve seen two different approaches:

1. The normal equation method: Best for smaller datasets, providing an exact solution through linear algebra. It’s fast and precise but can become computationally expensive with large datasets due to the matrix inversion operation.
2. Gradient descent: More suitable for large datasets, iteratively improving the parameters. While it requires choosing a learning rate and number of iterations, it’s more memory-efficient and often faster for big data.

Understanding these implementations helps build intuition about how linear regression works under the hood. While libraries like scikit-learn offer optimized implementations, building models from scratch with NumPy provides valuable insights into the mathematical foundations of machine learning. For more tutorials on linear algebra topics, please visit our collection here. 

For beginners looking to expand their knowledge, try experimenting with different learning rates in gradient descent, adding more features to the model, or introducing regularization to prevent overfitting. The beauty of implementing these algorithms from scratch is that it demystifies what might otherwise seem like “black box” machine learning models.