Linear Regression from Scratch in Python

Linear regression is the “hello world” of machine learning. It fits a straight line (or hyperplane) to your data by minimizing the distance between predictions and actual values. Building it from scratch teaches you the foundations that every ML model builds on.

The hypothesis

For simple linear regression with one feature:

$$y = wx + b$$

• w (weight/slope): how much y changes per unit of x
• b (bias/intercept): the value of y when x = 0

For multiple features: y = w₁x₁ + w₂x₂ + ... + wₙxₙ + b

The cost function

Mean Squared Error (MSE) measures how wrong our predictions are:

def cost_function(X, y, w, b):
    m = len(y)
    predictions = X @ w + b
    cost = (1 / (2 * m)) * np.sum((predictions - y) ** 2)
    return cost

The goal: find w and b that minimize this cost.

Gradient descent

Gradient descent iteratively adjusts w and b in the direction that reduces the cost.

def gradient_descent(X, y, w, b, learning_rate, iterations):
    m = len(y)
    history = []

    for i in range(iterations):
        predictions = X @ w + b
        errors = predictions - y

        dw = (1 / m) * (X.T @ errors)
        db = (1 / m) * np.sum(errors)

        w = w - learning_rate * dw
        b = b - learning_rate * db

        if i % 100 == 0:
            cost = cost_function(X, y, w, b)
            history.append(cost)

    return w, b, history

Full implementation

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X.squeeze() + np.random.randn(100) * 0.5

w = np.zeros(X.shape[1])
b = 0.0
learning_rate = 0.1
iterations = 1000

w, b, history = gradient_descent(X, y, w, b, learning_rate, iterations)
print(f"w = {w[0]:.4f}, b = {b:.4f}")
# Expected: w ≈ 3.0, b ≈ 4.0

plt.scatter(X, y, alpha=0.5)
x_line = np.linspace(0, 2, 100)
plt.plot(x_line, w[0] * x_line + b, color="red", linewidth=2)
plt.xlabel("X")
plt.ylabel("y")
plt.title("Linear Regression from Scratch")
plt.show()

Comparing to scikit-learn

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X, y)

print(f"sklearn w = {model.coef_[0]:.4f}, b = {model.intercept_:.4f}")
print(f"ours   w = {w[0]:.4f}, b = {b:.4f}")

The values should be nearly identical. scikit-learn uses the normal equation (closed-form solution) which is exact, while gradient descent converges to the same result iteratively.

Feature scaling

Gradient descent converges faster when features are on similar scales.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

Multiple features

X_multi = np.column_stack([
    np.random.rand(100),  # square footage
    np.random.rand(100),  # bedrooms
    np.random.rand(100),  # age
])
y_multi = 100 + 50 * X_multi[:, 0] + 30 * X_multi[:, 1] - 10 * X_multi[:, 2] + np.random.randn(100) * 5

w = np.zeros(3)
b = 0.0
w, b, history = gradient_descent(X_multi, y_multi, w, b, 0.1, 2000)
print(f"Weights: {w}")
print(f"Bias: {b:.2f}")

R² score

def r2_score(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)

predictions = X @ w + b
print(f"R² = {r2_score(y, predictions):.4f}")

Summary

Linear regression minimizes MSE using gradient descent (or the normal equation). The weight tells you the slope, the bias tells you the intercept, and R² tells you how well your line fits. Every neural network, logistic regression, and complex ML model uses these same primitives — cost functions, gradients, and iterative optimization.