Linear Regression From Scratch

import numpy as np

def fit_closed_form(X, y):
    X1 = np.hstack([X, np.ones((len(X), 1))])          # absorb the bias term
    return np.linalg.lstsq(X1, y, rcond=None)[0]        # lstsq, not inv — see below

def fit_gradient_descent(X, y, lr=0.05, steps=2000):
    X1 = np.hstack([X, np.ones((len(X), 1))])
    w = np.zeros(X1.shape[1])
    history = []
    for _ in range(steps):
        residual = X1 @ w - y                            # (n,)
        loss = float((residual ** 2).mean())
        grad = 2.0 / len(y) * (X1.T @ residual)          # (d+1,)
        w -= lr * grad
        history.append(loss)
    return w, history

Hold on to

• Closed-form solutions give you a ground truth for checking iterative ones
• Prefer a least-squares solve over an explicit matrix inverse
• The shape of a loss curve names the failure

Try each one before opening the solution. Getting it wrong first is most of where the learning happens.

1. Fit both ways on the same data and confirm the weights agree to four decimals.
   Hint

   Use enough gradient steps and a learning rate small enough not to oscillate.

   Solution

   With scaled features, a learning rate around 0.05 and a couple of thousand steps brings gradient descent within 1e-4 of the closed-form solution. If they disagree, the usual causes are unscaled features (which make the problem ill-conditioned and slow one direction to a crawl) or too few steps.

   import numpy as np
   
   rng = np.random.default_rng(0)
   X = rng.normal(size=(200, 3))
   true_w = np.array([1.5, -2.0, 0.5])
   y = X @ true_w + 4.0 + rng.normal(scale=0.1, size=200)
   
   w_closed = fit_closed_form(X, y)
   w_gd, history = fit_gradient_descent(X, y, lr=0.05, steps=5000)
   assert np.abs(w_closed - w_gd).max() < 1e-4

   Check your work

   Paste this after your own code. If it runs without raising, you have it.

   import numpy as np
   assert np.abs(w_closed - w_gd).max() < 1e-4
   assert history[-1] < history[0]
   print("ok")

2. Duplicate a feature column and observe what the closed form does. Explain using rank.
   Hint

   A duplicated column makes the design matrix rank-deficient.

   Solution

   With an exact duplicate the system is underdetermined: infinitely many weight vectors give identical predictions, because any amount can be shifted between the two identical columns. `np.linalg.lstsq` returns the minimum-norm solution, which splits the coefficient evenly between them rather than failing. An explicit `inv(X.T @ X)` would raise a singular-matrix error or return garbage. This is why you never read a coefficient in isolation when features are correlated — and it is the same rank deficiency that defines a robot singularity.

   import numpy as np
   
   X2 = np.hstack([X, X[:, :1]])            # column 0 duplicated
   print(np.linalg.matrix_rank(X2), X2.shape[1])   # 3 vs 4 — rank deficient
   
   w = np.linalg.lstsq(np.hstack([X2, np.ones((len(X2), 1))]), y, rcond=None)[0]
   print(w[0], w[3])                        # the coefficient split between the twins

3. Sweep the learning rate across four orders of magnitude and collect the five curve shapes above.
   Hint

   Try learning rates spanning something like 1e-5 to 10.

   Solution

   Too small: the curve is still steeply descending when the steps run out. Just right: smooth decay to a plateau. Slightly too large: oscillation, a sawtooth that descends unevenly. Too large: divergence to enormous values. Far too large: NaN within a few steps, because the overshoot compounds until the floats overflow. The boundary sits near 2 divided by the largest eigenvalue of the curvature.

   for lr in (1e-5, 1e-3, 0.05, 0.5, 2.0):
       _, history = fit_gradient_descent(X, y, lr=lr, steps=200)
       print(f"lr={lr:<8} first={history[0]:.3e} last={history[-1]:.3e}")

Sign in to track your progress through the lab.