S085 · Linear Regression from Scratch (numpy)

# linreg_scratch.py — linear regression, three ways, ~120 lines.
import numpy as np
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
 
np.random.seed(42)
 
# --- Data: 1000 rows, 5 features, some noise ---
X, y = make_regression(n_samples=1000, n_features=5, noise=10.0, random_state=42)
 
def add_bias(X):
    """Prepend a column of 1s so bias is folded into θ[0]."""
    return np.hstack([np.ones((X.shape[0], 1)), X])
 
Xb = add_bias(X)                   # shape (1000, 6)
print(f"X shape: {X.shape}, Xb (with bias): {Xb.shape}")
 
# =====================================================
# 1. NORMAL EQUATION — θ = (XᵀX)⁻¹ Xᵀy
# =====================================================
def fit_normal(X, y):
    """Closed-form OLS. One line if you don't count the docstring."""
    return np.linalg.inv(X.T @ X) @ X.T @ y
 
theta_normal = fit_normal(Xb, y)
yhat_normal = Xb @ theta_normal
print(f"[normal eq]   RMSE = {mean_squared_error(y, yhat_normal, squared=False):.4f}")
 
# =====================================================
# 2. BATCH GRADIENT DESCENT
# =====================================================
def fit_gd(X, y, lr=0.05, epochs=1000, verbose=False):
    """Batch gradient descent for MSE."""
    n, d = X.shape
    theta = np.zeros(d)                          # initialise at origin
    history = []
    for epoch in range(epochs):
        yhat = X @ theta                         # (n,)
        residual = yhat - y                      # (n,)
        grad = (1 / n) * (X.T @ residual)        # (d,)
        theta = theta - lr * grad                # update
        loss = 0.5 * np.mean(residual ** 2)      # MSE
        history.append(loss)
        if verbose and epoch % 100 == 0:
            print(f"  epoch {epoch:4d}: loss={loss:.4f}, ||grad||={np.linalg.norm(grad):.4f}")
    return theta, history
 
theta_gd, history = fit_gd(Xb, y, lr=0.05, epochs=1000, verbose=True)
yhat_gd = Xb @ theta_gd
print(f"[batch GD]    RMSE = {mean_squared_error(y, yhat_gd, squared=False):.4f}")
 
# =====================================================
# 3. MINI-BATCH STOCHASTIC GD
# =====================================================
def fit_sgd(X, y, lr=0.05, epochs=50, batch_size=32):
    """SGD — noisier gradient, faster convergence for huge datasets."""
    n, d = X.shape
    theta = np.zeros(d)
    for epoch in range(epochs):
        idx = np.random.permutation(n)
        for start in range(0, n, batch_size):
            batch = idx[start:start + batch_size]
            Xb_, yb_ = X[batch], y[batch]
            grad = (1 / len(batch)) * (Xb_.T @ (Xb_ @ theta - yb_))
            theta = theta - lr * grad
    return theta
 
theta_sgd = fit_sgd(Xb, y)
yhat_sgd = Xb @ theta_sgd
print(f"[mini-batch]  RMSE = {mean_squared_error(y, yhat_sgd, squared=False):.4f}")
 
# =====================================================
# 4. sklearn — sanity check
# =====================================================
sk = LinearRegression().fit(X, y)                # note: raw X, no bias col
theta_sk = np.concatenate([[sk.intercept_], sk.coef_])
yhat_sk = sk.predict(X)
print(f"[sklearn]     RMSE = {mean_squared_error(y, yhat_sk, squared=False):.4f}")
 
# =====================================================
# 5. Compare coefficients — should be near-identical
# =====================================================
print("\ncoefficient comparison (bias + 5 features):")
print(f"  normal:  {np.round(theta_normal, 3)}")
print(f"  GD:      {np.round(theta_gd, 3)}")
print(f"  SGD:     {np.round(theta_sgd, 3)}")
print(f"  sklearn: {np.round(theta_sk, 3)}")
 
# =====================================================
# 6. Learning-rate sanity check
# =====================================================
print("\n=== Learning rate demo ===")
for lr in [1e-4, 1e-2, 5e-2, 3e-1, 1.0]:
    _, hist = fit_gd(Xb, y, lr=lr, epochs=200)
    final = hist[-1]
    tag = "✅" if final < 100 else ("❌ diverged" if np.isnan(final) or final > 1e6 else "❌ slow")
    print(f"  lr={lr:>6}: final loss = {final:>10.2f}  {tag}")