Implement Linear Regression from Scratch

Math Foundations

Linear Regression predicts continuous output vector $y$ as a linear combination of input feature matrix $X \in \mathbb{R}^{N \times D}$ and weight vector $w \in \mathbb{R}^D$ plus bias scalar $b$:

$$\hat{y} = X w + b$$

The objective is to minimize Mean Squared Error (MSE) Loss:

$$\mathcal{L}(w, b) = \frac{1}{N} \sum_{i=1}^N (\hat{y}_i - y_i)^2 = \frac{1}{N} | X w + b - y |^2$$

Calculating partial derivatives with respect to $w$ and $b$:

$$\frac{\partial \mathcal{L}}{\partial w} = \frac{2}{N} X^T (\hat{y} - y), \quad \frac{\partial \mathcal{L}}{\partial b} = \frac{2}{N} \sum_{i=1}^N (\hat{y}_i - y_i)$$

Gradient Descent updates parameters iteratively with learning rate $\alpha$:

$$w \leftarrow w - \alpha \frac{\partial \mathcal{L}}{\partial w}, \quad b \leftarrow b - \alpha \frac{\partial \mathcal{L}}{\partial b}$$

NumPy Implementation from Scratch

import numpy as np


class LinearRegressionFromScratch:

  def __init__(self, lr=0.01, n_iters=1000):
    self.lr = lr
    self.n_iters = n_iters
    self.w = None
    self.b = None

  def fit(self, X, y):
    n_samples, n_features = X.shape
    # Initialize parameters
    self.w = np.zeros(n_features)
    self.b = 0.0

    # Gradient descent optimization loop
    for _ in range(self.n_iters):
      # 1. Forward Pass
      y_pred = np.dot(X, self.w) + self.b

      # 2. Compute Gradients
      error = y_pred - y
      dw = (2 / n_samples) * np.dot(X.T, error)
      db = (2 / n_samples) * np.sum(error)

      # 3. Update Parameters
      self.w -= self.lr * dw
      self.b -= self.lr * db

  def predict(self, X):
    return np.dot(X, self.w) + self.b

Closed Form Solution (Normal Equation)

Instead of iterative Gradient Descent, optimal weights can be solved directly using the Normal Equation:

$$w^* = (X^T X)^{-1} X^T y$$

• Pros: Exact analytical solution in a single step with zero hyperparameter tuning.
• Cons: Matrix inversion $(X^T X)^{-1}$ takes $\mathcal{O}(D^3)$ time, becoming slow when number of features $D > 10,000$.

Say this out loud

Implementing Linear Regression from scratch requires initializing weights to zero, computing linear forward predictions using matrix dot products, calculating Mean Squared Error loss gradients, and updating parameters iteratively via Gradient Descent. The Normal Equation offers a closed form alternative for smaller feature spaces.

Followups to expect

1. How do you add L2 Regularization (Ridge Regression) to the code? Add penalty term $2 \lambda w$ to the weight gradient calculation: dw = (2 / n_samples) * np.dot(X.T, error) + 2 * lambda_param * self.w.
2. Why vectorization with np.dot is faster than Python loops? Vectorized matrix operations execute low-level C and BLAS linear algebra instructions, utilizing CPU SIMD vector instructions for parallel speedup.