A Deep Dive into Linear Regression: Math, Code, and Intuition

Linear regression is the foundational building block of predictive modeling. It is the perfect entry point to understand how machine learning models map inputs to continuous numeric outputs.

In this deep dive, we will unpack linear regression from three perspectives: the mathematical formulations that govern it, a clean Python/NumPy implementation built from scratch, and the visual intuition of how it learns.

Linear Regression Deep Dive: 3D scatter plot with a neon-cyan linear regression fit line

---

1. The Mathematical Foundation

Linear regression assumes a linear relationship between a set of input features x and a continuous target variable y.

The Hypothesis Function

For a dataset with single-dimensional features, the hypothesis function is the equation of a straight line:

h_\theta(x) = \theta_0 + \theta_1 x

Where:

• \theta_0 is the intercept (bias).
• \theta_1 is the slope (weight).

For multi-dimensional features (where x is a vector), we represent the hypothesis in vectorized form by prepending a dummy feature x_0 = 1 to the input vector:

h_\theta(x) = \theta^T x = \theta_0 x_0 + \theta_1 x_1 + \dots + \theta_n x_n

The Cost Function (Mean Squared Error)

To measure how "wrong" our model predictions are, we define a cost function. The most common metric for regression is the Mean Squared Error (MSE). We scale it by \frac{1}{2m} to simplify derivatives during optimization:

J(\theta) = \frac{1}{2m} \sum_{i=1}^{m} \left( h_\theta(x^{(i)}) - y^{(i)} \right)^2

Where m is the total number of training examples. The goal of training is to find the parameter vector \theta that minimizes J(\theta).

Optimization: Gradient Descent

To minimize the cost function, we use Gradient Descent. Think of the cost function as a bowl-shaped terrain (a convex surface). We start with random parameters and take small steps in the direction of the steepest descent.

The parameters are updated iteratively:

\theta_j := \theta_j - \alpha \frac{\partial}{\partial \theta_j} J(\theta)

Where \alpha is the learning rate (step size). By calculating the partial derivative of the cost function, the update rule for each parameter \theta_j is:

\theta_j := \theta_j - \alpha \frac{1}{m} \sum_{i=1}^{m} \left( h_\theta(x^{(i)}) - y^{(i)} \right) x_j^{(i)}

---

2. Implementing Linear Regression from Scratch

Let's write a pure Python implementation of Linear Regression using NumPy to perform the gradient descent updates.

import numpy as np

class LinearRegressionGD:
    def __init__(self, learning_rate=0.01, epochs=1000):
        self.lr = learning_rate
        self.epochs = epochs
        self.weights = None
        self.bias = None
        self.loss_history = []

    def fit(self, X, y):
        n_samples, n_features = X.shape
        # Initialize weights to zeros
        self.weights = np.zeros(n_features)
        self.bias = 0.0

        for epoch in range(self.epochs):
            # Compute predictions (h_theta = Xw + b)
            y_predicted = np.dot(X, self.weights) + self.bias
            
            # Compute cost (Mean Squared Error)
            loss = (1 / (2 * n_samples)) * np.sum((y_predicted - y) ** 2)
            self.loss_history.append(loss)

            # Compute gradients
            dw = (1 / n_samples) * np.dot(X.T, (y_predicted - y))
            db = (1 / n_samples) * np.sum(y_predicted - y)

            # Update parameters
            self.weights -= self.lr * dw
            self.bias -= self.lr * db

            if epoch % (self.epochs // 10) == 0:
                print(f"Epoch {epoch:4d}: Cost = {loss:.6f}")

    def predict(self, X):
        return np.dot(X, self.weights) + self.bias

Testing the Implementation

Let's create dummy housing data and fit our scratch model:

# Generate synthetic dataset: House Size (sq ft / 1000) vs House Price ($100k)
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
y = y.squeeze()  # Flatten target into 1D array

# Instantiate and fit the model
model = LinearRegressionGD(learning_rate=0.1, epochs=500)
model.fit(X, y)

print(f"\nLearned Weights (Slope): {model.weights[0]:.4f}")
print(f"Learned Bias (Intercept): {model.bias:.4f}")
# Expected outputs: Slope close to 3, Intercept close to 4

---

3. Developing Intuition

The Learning Rate (\alpha)

Choosing the right learning rate is critical:

• Too Small: Gradient descent will eventually find the minimum, but it will take a very long time, requiring high computational resources.
• Too Large: The algorithm can overshoot the minimum, bounce back and forth, and actually diverge (the cost increases with every step).

         Large Learning Rate (Divergence)         Small Learning Rate (Slow)
               \             /                             \
                \   *       /                               *
                 \ / \     /                                 \
                  *   \   /                                   *
                       \ /                                     \
                        * (Minimum)                             * (Minimum)

The Geometry of Loss

For a single feature regression, the loss function J(\theta_0, \theta_1) is a 3D paraboloid. Gradient descent calculates the vector direction of steepest slope at our current coordinates, then moves our parameters down that slope. Because the MSE cost function for linear regression is convex, we are mathematically guaranteed to find the global minimum if our learning rate is chosen correctly.

Conclusion

Linear regression exemplifies the core mechanics of machine learning: defining a model hypothesis, setting up a mathematical cost function, and optimization through gradient descent. Whether you run a simple linear model in Scikit-Learn or build multi-million-parameter neural networks, the underlying optimization loop of evaluating loss and applying gradients remains exactly the same.