The cost function and gradient descent

w = np.array([55.0, 7.5, -28.0, 90.0])

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

df = pd.read_csv(DATA / "uy_ijara.csv")
num = ["maydon_m2"]                      # ONE feature first, easier to see
X = df[num].to_numpy(dtype=float)
y = df["narx_usd"].to_numpy(dtype=float)

def predict(X, w, b):
    return X @ w + b

def mse(X, y, w, b):
    err = predict(X, w, b) - y
    return float(np.mean(err ** 2))

# a few different guesses
for w_t, b_t in [(np.array([5.0]), 200.0), (np.array([7.5]), 180.0), (np.array([12.0]), 50.0)]:
    print(f"w={w_t[0]:>5.1f}  b={b_t:>5.0f}  ->  MSE = {mse(X, y, w_t, b_t):>12,.0f}".replace(",", " "))

show(pd.DataFrame([
    {"Reason": "Removing the sign", "Note": "A +50 and a -50 error must not cancel each other out"},
    {"Reason": "Punishing big errors more", "Note": "An error of 10 -> 100; an error of 100 -> 10 000"},
    {"Reason": "An easy derivative", "Note": "The derivative of a square is linear, so gradients are convenient"},
]))

w_range = np.linspace(0, 20, 200)
cost = [mse(X, y, np.array([wi]), 180.0) for wi in w_range]
best_i = int(np.argmin(cost))

figure, ax = plt.subplots(figsize=(9.5, 5))
ax.plot(w_range, cost, "b-", lw=2)
ax.plot(w_range[best_i], cost[best_i], "ro", ms=11,
        label=f"lowest: w={w_range[best_i]:.2f}")
ax.set_xlabel("w (the weight of maydon_m2)")
ax.set_ylabel("MSE (cost)")
ax.set_title("The cost function is shaped like a bowl", fontweight="bold")
ax.legend(); ax.grid(alpha=.3)
fig("kosa.png", caption="MSE is shaped like a bowl. The lowest point is the best weight. Gradient descent walks down to it.")

print(f"Lowest point: w = {w_range[best_i]:.3f},  MSE = {cost[best_i]:,.0f}".replace(",", " "))

def gradient(X, y, w, b):
    """The derivatives of MSE with respect to w and b."""
    n = len(y)
    err = predict(X, w, b) - y            # shape (n,)
    dw = (2 / n) * (X.T @ err)            # shape (m,)
    db = (2 / n) * np.sum(err)
    return dw, db

for wi in [3.0, 7.5, 15.0]:
    dw, db = gradient(X, y, np.array([wi]), 180.0)
    direction = "w should GO UP" if dw[0] < 0 else "w should GO DOWN"
    print(f"w={wi:>5.1f}  ->  dw={dw[0]:>12,.1f}   {direction}".replace(",", " "))

def gradient_descent(X, y, alpha=0.0001, steps=200):
    w = np.zeros(X.shape[1])
    b = 0.0
    history = []
    for i in range(steps):
        dw, db = gradient(X, y, w, b)
        w = w - alpha * dw
        b = b - alpha * db
        history.append({"step": i, "w": float(w[0]), "b": b, "mse": mse(X, y, w, b)})
    return w, b, pd.DataFrame(history)

w_fin, b_fin, history = gradient_descent(X, y, alpha=0.0002, steps=400)
show(history.iloc[[0, 1, 2, 5, 20, 100, 399]].round(3))

print(f"Start : w=0.000     MSE={history.mse.iloc[0]:>12,.0f}".replace(",", " "))
print(f"End   : w={w_fin[0]:.3f}   MSE={history.mse.iloc[-1]:>12,.0f}".replace(",", " "))
print()
from sklearn.linear_model import LinearRegression
lr = LinearRegression().fit(X, y)
print(f"sklearn's answer: w={lr.coef_[0]:.3f}  b={lr.intercept_:.2f}")
print(f"Our answer      : w={w_fin[0]:.3f}  b={b_fin:.2f}")

print("How many steps unscaled data needs:")
print(f"{'steps':>10} {'w':>9} {'b':>10} {'MSE':>12}")
print("-" * 44)
for n in [400, 5000, 50000, 300000]:
    w_t, b_t, h_t = gradient_descent(X, y, alpha=0.0002, steps=n)
    mse_text = f"{h_t.mse.iloc[-1]:,.0f}".replace(",", " ")
    print(f"{n:>10} {w_t[0]:>9.3f} {b_t:>10.2f} {mse_text:>12}")
print()
print(f"sklearn:   w={lr.coef_[0]:.3f}  b={lr.intercept_:.2f}")

figure, axes = plt.subplots(1, 2, figsize=(14, 4.8))
axes[0].plot(history.step, history.mse, "b-", lw=2)
axes[0].set_xlabel("Step"); axes[0].set_ylabel("MSE")
axes[0].set_title("The cost falls step by step", fontweight="bold")
axes[0].grid(alpha=.3)

axes[1].plot(w_range, cost, "k-", lw=1.5, alpha=.5, label="cost function")
some_steps = history.iloc[::12]
axes[1].plot(some_steps.w, [mse(X, y, np.array([wi]), 180.0) for wi in some_steps.w],
             "ro-", ms=5, lw=1, label="gradient descent path")
axes[1].set_xlabel("w"); axes[1].set_ylabel("MSE")
axes[1].set_title("Walking down the bowl", fontweight="bold")
axes[1].legend(); axes[1].grid(alpha=.3)
fig("gd_yoli.png", caption="Left: the cost goes down at every step. Right: what those steps look like on the bowl.")

print(f"{'alpha':>10} {'w after 400 steps':>22} {'MSE':>16} {'status':>14}")
print("-" * 68)
for a in [0.000001, 0.00002, 0.0002, 0.001, 0.01]:
    w_t, b_t, h = gradient_descent(X, y, alpha=a, steps=400)
    m = h.mse.iloc[-1]
    if not np.isfinite(m) or m > 1e12:
        status, m_text = "EXPLODED", "inf"
    elif abs(w_t[0] - lr.coef_[0]) < 0.5:
        status, m_text = "converged", f"{m:,.0f}".replace(",", " ")
    else:
        status, m_text = "too slow", f"{m:,.0f}".replace(",", " ")
    print(f"{a:>10} {w_t[0]:>22.3f} {m_text:>16} {status:>14}")
print()
print(f"sklearn's answer: w = {lr.coef_[0]:.3f}")

from sklearn.preprocessing import StandardScaler
Xs = StandardScaler().fit_transform(X)
print("On scaled data a MUCH bigger alpha works:")
print(f"{'alpha':>8} {'w':>10} {'MSE':>14}")
print("-" * 34)
for a in [0.001, 0.01, 0.1, 0.5]:
    w_t, b_t, h = gradient_descent(Xs, y, alpha=a, steps=400)
    print(f"{a:>8} {w_t[0]:>10.2f} {h.mse.iloc[-1]:>13,.0f}".replace(",", " "))

num4 = ["xonalar", "maydon_m2", "metro_km", "mebel"]
X4 = StandardScaler().fit_transform(df[num4].to_numpy(dtype=float))
w4, b4, h4 = gradient_descent(X4, y, alpha=0.1, steps=1500)

lr4 = LinearRegression().fit(X4, y)
compare = pd.DataFrame({
    "feature": num4,
    "we found": w4.round(2),
    "sklearn": lr4.coef_.round(2),
})
show(compare)
print(f"bias: ours {b4:.2f}  |  sklearn {lr4.intercept_:.2f}")
print(f"MSE : ours {h4.mse.iloc[-1]:,.0f}  |  sklearn {mse(X4, y, lr4.coef_, lr4.intercept_):,.0f}".replace(",", " "))

show(pd.DataFrame([
    {"Kind": "Batch GD", "Per step": "ALL rows", "Speed": "Slow",
     "Path": "Smooth", "When": "Small data"},
    {"Kind": "Stochastic GD (SGD)", "Per step": "ONE row", "Speed": "Fast",
     "Path": "Very noisy", "When": "Huge data, online training"},
    {"Kind": "Mini-batch GD", "Per step": "32-256 rows", "Speed": "Fast",
     "Path": "Moderately noisy", "When": "The MOST USED in practice"},
]))

def mini_batch_gd(X, y, alpha=0.1, epochs=30, batch=64, seed=42):
    r = np.random.default_rng(seed)
    w, b = np.zeros(X.shape[1]), 0.0
    history = []
    n = len(y)
    for e in range(epochs):
        order = r.permutation(n)
        for start in range(0, n, batch):
            idx = order[start:start + batch]
            dw, db = gradient(X[idx], y[idx], w, b)
            w, b = w - alpha * dw, b - alpha * db
        history.append(mse(X, y, w, b))
    return w, b, history

w_mb, b_mb, hist_mb = mini_batch_gd(X4, y)
print(f"Batch GD      (1500 steps): MSE = {h4.mse.iloc[-1]:>10,.0f}".replace(",", " "))
print(f"Mini-batch GD (30 epochs) : MSE = {hist_mb[-1]:>10,.0f}".replace(",", " "))
print(f"sklearn                   : MSE = {mse(X4, y, lr4.coef_, lr4.intercept_):>10,.0f}".replace(",", " "))

show(pd.DataFrame([
    {"Where": "Linear regression (Chapter 15)", "Role of gradient descent": "Finds the weights"},
    {"Where": "Logistic regression (Chapter 16)", "Role of gradient descent": "The same, with another cost function"},
    {"Where": "Ridge / Lasso (Chapter 17)", "Role of gradient descent": "A penalty is added to the cost"},
    {"Where": "Neural network (Chapter 31)", "Role of gradient descent": "Backpropagation is the gradient along a chain"},
    {"Where": "XGBoost (Chapter 20)", "Role of gradient descent": "It's in the name: Gradient Boosting"},
]))

show(pd.DataFrame([
    {"Concept": "Cost function", "In one sentence": "Measures the error as one number (MSE)"},
    {"Concept": "Why squared", "In one sentence": "Removes the sign, punishes big errors, easy derivative"},
    {"Concept": "The cost is a bowl", "In one sentence": "The task: find the bottom of the bowl"},
    {"Concept": "Gradient", "In one sentence": "The slope: its sign is the direction, its size the steepness"},
    {"Concept": "The GD rule", "In one sentence": "w := w - alpha * gradient"},
    {"Concept": "Learning rate", "In one sentence": "Small is slow, big explodes; the window is narrow"},
    {"Concept": "Why scaling is needed", "In one sentence": "It makes the bowl round, so one alpha works everywhere"},
    {"Concept": "Mini-batch", "In one sentence": "The variant used most in practice"},
]))