From Equation to Implementation: Reproducing Linear Regression

Why this project

Every idea in this course — vectors, matrices, dot products, derivatives, partial derivatives, gradients, and the mechanics of NumPy — exists so that you can build a model that learns from data. Linear regression is the simplest model where all of these pieces come together, and it is the honest foundation underneath neural networks, logistic regression, and most of classical machine learning. If you can derive it on paper and reproduce it in code from scratch, you understand the machinery, not just the API call.

In this capstone you will implement linear regression end to end with NumPy only — no scikit-learn, no PyTorch, no TensorFlow. You will derive the loss and its gradients by hand, translate those equations into vectorized code, verify your gradients numerically, train with gradient descent, and study how the learning rate and the data itself shape the outcome.

The model

Given $n$ examples, each with $d$ features, we predict a real number per example with an affine model:

Written per example, $\hat{y}_i = \sum_{j=1}^{d} X_{ij} w_j + b$.

The loss

We measure error with mean squared error (MSE):

Training means choosing $w$ and $b$ to make $L$ small.

Every variable and its shape

• $X$ — feature matrix, shape $(n, d)$. Row $i$ is example $i$; column $j$ is feature $j$.
• $y$ — target vector, shape $(n,)$. The true value for each example.
• $w$ — weight vector, shape $(d,)$. One weight per feature.
• $b$ — bias (intercept), a scalar. Shifts the whole prediction.
• $\hat{y}$ — prediction vector, shape $(n,)$, equal to $X w + b$.
• $r = \hat{y} - y$ — residual vector, shape $(n,)$.
• $L$ — the loss, a single scalar.
• $\frac{\partial L}{\partial w}$ (written $dw$) — gradient w.r.t. weights, shape $(d,)$.
• $\frac{\partial L}{\partial b}$ (written $db$) — gradient w.r.t. bias, a scalar.
• $\eta$ — the learning rate, a positive scalar step size.

Keeping shapes straight is the single most important habit here: $X w$ contracts the $d$ axis, leaving $(n,)$; broadcasting $+ b$ keeps it $(n,)$; and $X^{\top} r$ contracts the $n$ axis, giving back a $(d,)$ gradient.

The gradients you will derive

Differentiating $L = \frac{1}{n} \sum_i (\hat{y}_i - y_i)^2$ and using the chain rule with $\hat{y}_i = \sum_j X_{ij} w_j + b$ gives the two results that this whole project rests on:

Gradient descent then repeatedly steps downhill:

The plan

You will build this up one verifiable step at a time:

1. Implement prediction — code $\hat{y} = X w + b$ and confirm it returns shape $(n,)$.
2. Derive MSE — write the loss and implement it as a scalar.
3. Derive gradients — differentiate $L$ by hand to get $dw$ and $db$.
4. Implement gradient descent — one update step for $w$ and $b$.
5. Write the training loop — repeat the step for many epochs, recording the loss.
6. Gradient check — compare your analytic gradients to a central finite-difference estimate; the max absolute difference must be below $10^{-5}$.
7. Plot predictions — scatter $y$ vs $\hat{y}$ (or fit line on 1-D data).
8. Plot the loss curve — $L$ per epoch should fall smoothly.
9. Experiment with learning rates — watch too-small crawl, good converge, and too-large diverge.
10. Add multiple features — confirm everything generalizes from $d = 1$ to $d > 1$ with no code changes.
11. Analyze failure cases — high noise, bad scaling, collinear features, exploding learning rates.
12. Write the report — a short research-style write-up with concrete numbers.

Compare your gradient-descent solution against the closed-form normal equations to prove your optimizer converges to the right place.