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

One practical note that the from-scratch implementation above quietly hides: it only converges nicely because the synthetic feature is already on a [0, 2] scale. Feed it a real housing dataset with square footage in the thousands and bedroom counts in single digits, and the same learning rate that worked here will explode the cost to NaN within twenty epochs. Gradient descent is extremely sensitive to feature scaling — the loss surface becomes a long narrow valley, and the step that is reasonable along one axis is catastrophic along another. Standardize your features before training, or watch your loss history and cut the learning rate when it goes up instead of down.

Also worth being honest about: for ordinary least squares you usually should not use gradient descent at all. The normal equation \theta = (X^T X)^{-1} X^T y gives the exact optimum in one step, and it is the faster choice for anything under roughly 10,000 features. Gradient descent earns its keep only when the matrix inversion becomes intractable or when the model is no longer convex — which is to say, when you have moved on to neural networks. Implementing it by hand here is worth doing for the intuition, not because it is the right tool for a straight line.

The deeper limitation is in the word "linear" itself. If the true relationship curves, no amount of tuning fixes it; you will get a well-optimized model that is systematically wrong. Plot your residuals against the fitted values before you trust any regression. If they show structure rather than random scatter, the model is telling you it needs a different hypothesis, not a better learning rate.