Linear Regression

Method of Least Squares

This is the closed-form solution — instead of searching for the best weights step by step, we solve for them directly. The weights below are the exact ones that minimize the sum of squared errors, derived by setting the derivative of the cost function to zero and solving algebraically.

import numpy as np

def fit(X, y):
    X = np.insert(X, 0, 1, axis=1)
    XT_X = np.dot(X.T, X)
    XT_X_inv = np.linalg.inv(XT_X)
    XT_y = np.dot(X.T, y)
    betas = np.dot(XT_X_inv, XT_y)
    return betas[0], betas[1:]

It's elegant, and for small datasets it's the fastest path to an exact answer. The catch is the matrix inversion $(\mathbf{X}^T\mathbf{X})^{-1}$, which costs roughly $O(n^3)$ as the number of features $n$ grows, and simply fails when $\mathbf{X}^T\mathbf{X}$ isn't invertible. That's the gap gradient descent fills: instead of solving for the answer in one shot, it approaches it gradually, which scales far better to large datasets and large feature counts.

Gradient Descent Method

Imagine standing on a hillside in thick fog, trying to reach the lowest point in the valley. You can't see the bottom, but you can feel which direction the ground slopes beneath your feet — so you take a small step downhill, feel the slope again, and repeat. That's gradient descent. The "hill" is the cost function, the "slope" is its gradient with respect to the weights, and each step nudges the weights a little closer to the values that minimize prediction error.

Cost Function (Mean Squared Error):

To know which way is "downhill," we need the slope of this cost function with respect to each weight — that's the gradient. Here's how it's derived, one step at a time.

Deriving the Gradient:

Weight and Bias Updates:

The learning rate $\alpha$ controls how big each step is. Set it too small and training crawls toward the minimum, needing far more iterations than necessary. Set it too large and the weights overshoot the minimum entirely, sometimes bouncing back and forth or diverging instead of converging — choosing $\alpha$ well is one of the most practical decisions in training any gradient-based model.

There's still an open question, though: when we compute that gradient, how much data do we use to compute it? That choice splits gradient descent into three variants.

Batch Gradient Descent: computes the gradient using the entire training set before making a single update. The path it takes toward the minimum is smooth and stable, since every step is based on the full picture — but for large datasets, recomputing the gradient over millions of rows for every single update is slow and memory-hungry.

Mini-Batch Gradient Descent: splits the training data into small batches — commonly 32, 64, or 128 samples — and updates the weights once per batch. This is the workhorse of modern machine learning: it's far faster than batch gradient descent since each update only needs a slice of the data, and the small amount of noise from batch-to-batch variation actually helps the optimizer avoid getting stuck in shallow local minima, while still being far more stable than updating on a single point at a time.

Stochastic Gradient Descent (SGD): takes mini-batching to its extreme — the batch size is just one. The weights update after every single training example, making each individual step fast and well-suited to streaming or online learning where data arrives continuously. The trade-off is a noisier path to the minimum: the cost doesn't decrease smoothly but jitters as it descends, since every update is based on the (sometimes misleading) gradient of just one data point.

In short: batch trades speed for stability, stochastic trades stability for speed, and mini-batch sits in between — which is why it's the default choice in practice. The NumPy implementation below is stochastic gradient descent: notice that the weight update happens inside the inner loop, once per training example, rather than once per full pass over the dataset.

import numpy as np

def fit(X, y, epochs, learning_rate):
    n_samples, n_features = X.shape
    X_aug = np.insert(X, 0, 1, axis=1)
    weights = np.random.randn(n_features + 1) * 0.01

    for epoch in range(epochs):
        for i in range(n_samples):
            y_pred = np.dot(X_aug[i], weights)
            error = y_pred - y[i]
            gradient = error * X_aug[i]
            weights -= learning_rate * gradient
    return weights