The normal equation vs gradient descent

import pandas as pd
import numpy as np
import time

df = pd.read_csv(DATA / "uy_ijara.csv")
num = ["xonalar", "maydon_m2", "metro_km", "mebel"]
X_raw = df[num].to_numpy(dtype=float)
y = df["narx_usd"].to_numpy(dtype=float)

# add a column of ones for the bias
X = np.c_[np.ones(len(X_raw)), X_raw]
print(f"X shape (with the bias column): {X.shape}")
print(X[:3].round(2))

def normal_equation(X, y):
    return np.linalg.inv(X.T @ X) @ X.T @ y

t0 = time.perf_counter()
theta = normal_equation(X, y)
ne_time = time.perf_counter() - t0

print("Normal equation result:")
print(f"  bias (b)  : {theta[0]:>8.3f}")
for name, v in zip(num, theta[1:]):
    print(f"  {name:10s}: {v:>8.3f}")
print()
print(f"Compute time: {ne_time*1000:.3f} ms  (in ONE go)")

from sklearn.linear_model import LinearRegression

lr = LinearRegression().fit(X_raw, y)
compare = pd.DataFrame({
    "parameter": ["bias"] + num,
    "normal equation": np.round(theta, 4),
    "sklearn": np.round(np.r_[lr.intercept_, lr.coef_], 4),
})
compare["diff"] = (compare["normal equation"] - compare["sklearn"]).abs().round(10)
show(compare)
print()
print("The same:", bool(np.allclose(theta, np.r_[lr.intercept_, lr.coef_])))

from sklearn.preprocessing import StandardScaler

def mse(X, y, theta):
    return float(np.mean((X @ theta - y) ** 2))

def gd(X, y, alpha, steps):
    theta = np.zeros(X.shape[1])
    for _ in range(steps):
        grad = (2 / len(y)) * X.T @ (X @ theta - y)
        theta = theta - alpha * grad
    return theta

# GD works well on scaled data (Chapter 13)
Xs_raw = StandardScaler().fit_transform(X_raw)
Xs = np.c_[np.ones(len(Xs_raw)), Xs_raw]

theta_ne_s = normal_equation(Xs, y)
print(f"{'method':<28} {'MSE':>12} {'time (ms)':>11}")
print("-" * 53)
t0 = time.perf_counter(); tn = normal_equation(Xs, y); tv = (time.perf_counter()-t0)*1000
print(f"{'Normal equation':<28} {mse(Xs, y, tn):>12,.1f} {tv:>11.3f}".replace(",", " "))
for n in [100, 1000, 10000]:
    t0 = time.perf_counter(); tg = gd(Xs, y, 0.1, n); tv = (time.perf_counter()-t0)*1000
    print(f"{'GD (' + str(n) + ' steps)':<28} {mse(Xs, y, tg):>12,.1f} {tv:>11.3f}".replace(",", " "))

print("How long the normal equation takes as the number of columns grows:")
print(f"{'columns':>10} {'X.T@X shape':>16} {'normal equation':>18} {'GD (200 steps)':>17}")
print("-" * 66)
r = np.random.default_rng(42)
for m in [10, 100, 500, 1500]:
    n = 3000
    Xb = r.normal(size=(n, m))
    yb = Xb @ r.normal(size=m) + r.normal(scale=.5, size=n)
    t0 = time.perf_counter(); normal_equation(Xb, yb); nt = (time.perf_counter()-t0)*1000
    t0 = time.perf_counter(); gd(Xb, yb, 0.01, 200); gt = (time.perf_counter()-t0)*1000
    print(f"{m:>10} {f'{m}x{m}':>16} {nt:>15.1f} ms {gt:>14.1f} ms")

# two columns are IDENTICAL (perfectly correlated); this happens in real data too
X_broken = np.c_[np.ones(100), r.normal(size=100)]
X_broken = np.c_[X_broken, X_broken[:, 1]]       # column 3 is a copy of column 2
y_broken = r.normal(size=100)

print(f"X shape: {X_broken.shape}")
print(f"determinant of X.T @ X: {np.linalg.det(X_broken.T @ X_broken):.3e}   <- very close to zero")
print()
try:
    np.linalg.inv(X_broken.T @ X_broken)
    print("inv() worked (but the result is unreliable)")
except np.linalg.LinAlgError as e:
    print(f"inv() raised an ERROR: {e}")
print()
print("THE FIX: pinv(), the pseudo-inverse (sklearn uses it too)")
theta_p = np.linalg.pinv(X_broken) @ y_broken
print(f"  pinv() worked, result: {theta_p.round(4)}")

print("Do the three methods give the same answer on our data:")
t1 = np.linalg.inv(X.T @ X) @ X.T @ y
t2 = np.linalg.pinv(X) @ y
t3 = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"  inv  : {t1.round(4)}")
print(f"  pinv : {t2.round(4)}")
print(f"  lstsq: {t3.round(4)}")
print(f"  All agree: {bool(np.allclose(t1, t2) and np.allclose(t2, t3))}")

show(pd.DataFrame([
    {"Property": "Iterations", "Normal equation": "NONE, one go", "Gradient descent": "Thousands of steps"},
    {"Property": "Learning rate", "Normal equation": "NOT NEEDED", "Gradient descent": "Must be chosen, delicate"},
    {"Property": "Scaling", "Normal equation": "Not needed", "Gradient descent": "REQUIRED (Chapter 13)"},
    {"Property": "Many columns (m > 10 000)", "Normal equation": "Very slow, O(m^3)", "Gradient descent": "Works"},
    {"Property": "Many rows (n > 1 million)", "Normal equation": "Memory problems", "Gradient descent": "Works with mini-batches"},
    {"Property": "Data arrives as a stream", "Normal equation": "Doesn't work", "Gradient descent": "Works with SGD"},
    {"Property": "Multicollinearity", "Normal equation": "inv() falls over (needs pinv)", "Gradient descent": "Works"},
    {"Property": "Fits other models", "Normal equation": "ONLY linear + MSE", "Gradient descent": "Logistic, NN, boosting: all of them"},
]))

show(pd.DataFrame([
    {"Model": "Linear Regression", "Closed-form formula?": "YES", "What sklearn uses": "lstsq (closed form)"},
    {"Model": "Ridge", "Closed-form formula?": "YES", "What sklearn uses": "Closed form or iterative (solver=)"},
    {"Model": "Lasso", "Closed-form formula?": "NO", "What sklearn uses": "Coordinate descent"},
    {"Model": "Logistic Regression", "Closed-form formula?": "NO", "What sklearn uses": "lbfgs / newton-cg (iterative)"},
    {"Model": "Neural network", "Closed-form formula?": "NO", "What sklearn uses": "Mini-batch GD (Adam)"},
    {"Model": "XGBoost", "Closed-form formula?": "NO", "What sklearn uses": "Gradient boosting"},
]))

show(pd.DataFrame([
    {"Situation": "Small data, few columns", "Choice": "Normal equation (LinearRegression)"},
    {"Situation": "More than 10 000 columns", "Choice": "SGDRegressor or Ridge(solver='sag')"},
    {"Situation": "Millions of rows", "Choice": "SGDRegressor (mini-batch)"},
    {"Situation": "Data arrives as a stream", "Choice": "SGDRegressor + partial_fit"},
    {"Situation": "A non-linear model", "Choice": "No other way: gradient descent"},
]))

from sklearn.linear_model import SGDRegressor
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

print("Two sklearn models, the same data:")
lr_scores = cross_val_score(LinearRegression(), Xs_raw, y, cv=5, scoring="r2")
sgd_scores = cross_val_score(SGDRegressor(max_iter=5000, tol=1e-4, random_state=42),
                             Xs_raw, y, cv=5, scoring="r2")
print(f"  LinearRegression (closed form)  : R2 = {lr_scores.mean():.4f} ±{lr_scores.std():.4f}")
print(f"  SGDRegressor     (gradient desc): R2 = {sgd_scores.mean():.4f} ±{sgd_scores.std():.4f}")
print()
print(f"Gap: {abs(lr_scores.mean() - sgd_scores.mean()):.4f}  |  spread: ±{max(lr_scores.std(), sgd_scores.std()):.4f}")
print("Verdict:", "the gap is within the spread, so both reach the same answer"
      if abs(lr_scores.mean()-sgd_scores.mean()) < max(lr_scores.std(), sgd_scores.std())
      else "the gap is noticeable")

show(pd.DataFrame([
    {"Concept": "Normal equation", "In one sentence": "w = (X^T X)^-1 X^T y: the exact answer in one go"},
    {"Concept": "Inside LinearRegression", "In one sentence": "NOT gradient descent: a closed-form solution (lstsq)"},
    {"Concept": "The danger of inv()", "In one sentence": "Fails under multicollinearity; use pinv() or lstsq()"},
    {"Concept": "Cost of the normal equation", "In one sentence": "O(m^3): expensive with many columns"},
    {"Concept": "Cost of GD", "In one sentence": "Linear in the number of columns: cheaper for large m"},
    {"Concept": "Why GD is universal", "In one sentence": "A closed form exists only for linear + MSE"},
]))