Intro to ML Lab

# Learning as optimization: choose f from a model class to minimize the sample loss (curve fitting)
import numpy as np

rng = np.random.default_rng(0)
truth = lambda x: np.sin(2 * np.pi * x)
x_test = np.linspace(0, 1, 200); y_test = truth(x_test)


def sample(n):
    x = rng.random(n)
    return x, truth(x) + 0.25 * rng.standard_normal(n)


def fit(x, y, degree):
    """Least squares: the f in the class of degree-d polynomials that minimizes the mean squared error."""
    V = np.vander(x, degree + 1)
    w, *_ = np.linalg.lstsq(V, y, rcond=None)
    return lambda t: np.vander(t, degree + 1) @ w


mse = lambda f, x, y: np.mean((f(x) - y) ** 2)
x, y = sample(15)
print("15 training points")
print(" degree  train MSE  test MSE")
for d in (1, 3, 5, 9, 12):
    f = fit(x, y, d)
    print(f"  {d:2d}     {mse(f, x, y):.3f}     {mse(f, x_test, y_test):.3f}")      # high degrees overfit

print("degree 9 with more data (the unreasonable effectiveness of data)")
for n in (15, 50, 200, 1000):
    xs, ys = sample(n)
    print(f"  n = {n:4d}: test MSE {mse(fit(xs, ys, 9), x_test, y_test):.3f}")

# The same idea with gradient descent: a line y = a x + b fitted by minimizing the sample score S(f)
x, y = sample(50)
a = b = 0.0
for step in range(2000):
    err = a * x + b - y
    a -= 0.5 * np.mean(err * x); b -= 0.5 * np.mean(err)
print(f"gradient descent line: y = {a:.2f} x + {b:.2f}, MSE {np.mean((a * x + b - y) ** 2):.3f}")

# One dataset, three learning paradigms: supervised, unsupervised and semi-supervised (Figure 11.2)
import numpy as np

rng = np.random.default_rng(1)
X = np.vstack([rng.normal([-1.5, 0], 0.7, (100, 2)), rng.normal([1.5, 0.5], 0.7, (100, 2))])
y = np.repeat([0, 1], 100)

# Supervised: logistic regression learns from all the labels
w, b = np.zeros(2), 0.0
for _ in range(2000):
    p = 1 / (1 + np.exp(-(X @ w + b)))
    w -= 0.1 * X.T @ (p - y) / len(y); b -= 0.1 * np.mean(p - y)
sup = (X @ w + b > 0).astype(int)
print(f"supervised (all 200 labels):      accuracy {np.mean(sup == y):.1%}")

# Unsupervised: k-means sees no labels at all; it can only group the points
C = X[rng.choice(len(X), 2, replace=False)]
for _ in range(50):
    lab = ((X[:, None] - C) ** 2).sum(-1).argmin(1)
    C = np.array([X[lab == j].mean(0) for j in range(2)])
acc = max(np.mean(lab == y), np.mean(lab != y))           # cluster ids are arbitrary
print(f"unsupervised (no labels):         agreement with the classes {acc:.1%}")

# Semi-supervised: 4 labels, then self-training with pseudo-labels for the rest
known = np.array([0, 1, 100, 101])
labels = np.full(len(y), -1); labels[known] = y[known]
for _ in range(10):
    mask = labels >= 0
    cent = np.array([X[mask & (labels == k)].mean(0) for k in (0, 1)])
    d = ((X[:, None] - cent) ** 2).sum(-1)
    conf = np.abs(d[:, 0] - d[:, 1])
    unl = np.flatnonzero(~mask)
    if not len(unl):
        break
    pick = unl[np.argsort(-conf[unl])[:20]]             # pseudo-label the 20 most confident points
    labels[pick] = d[pick].argmin(1)
semi = ((X[:, None] - np.array([X[labels == k].mean(0) for k in (0, 1)])) ** 2).sum(-1).argmin(1)
print(f"semi-supervised (4 labels):       accuracy {np.mean(semi == y):.1%}")

# ML for optimization, "learning to configure algorithms" (Figure 11.14): a bandit picks the step size of a search
import math
import random

def rastrigin(x):
    return 10 * len(x) + sum(v * v - 10 * math.cos(2 * math.pi * v) for v in x)

STEPS = [0.02, 0.1, 0.5, 2.0]


def search(choose, rng, iters=3000, dim=5):
    """A (1+1) random search: perturb the current point with Gaussian noise of a chosen size, keep it if better."""
    x = [rng.uniform(-5.12, 5.12) for _ in range(dim)]; fx = rastrigin(x)
    n, gain = [0] * len(STEPS), [0.0] * len(STEPS)
    for t in range(iters):
        k = choose(t, n, gain, rng)
        y = [min(5.12, max(-5.12, v + rng.gauss(0, STEPS[k]))) for v in x]
        fy = rastrigin(y)
        reward = max(0.0, fx - fy)                       # reward = improvement
        n[k] += 1; gain[k] += reward
        if fy < fx:
            x, fx = y, fy
    return fx, n


def ucb(t, n, gain, rng, c=0.5):
    for k in range(len(STEPS)):
        if n[k] == 0:
            return k
    return max(range(len(STEPS)), key=lambda k: gain[k] / n[k] + c * math.sqrt(math.log(t + 1) / n[k]))


rng = random.Random(0)
fixed = []
for k, s in enumerate(STEPS):
    res = [search(lambda t, n, g, r, k=k: k, rng)[0] for _ in range(20)]
    fixed.append(sum(res) / 20)
    print(f"fixed step {s:4.2f}: mean best f = {fixed[-1]:6.2f}")
print(f"a fixed step picked blindly averages {sum(fixed) / len(fixed):6.2f}")
runs = [search(ucb, rng) for _ in range(20)]
print(f"UCB-chosen step: mean best f = {sum(r[0] for r in runs) / 20:6.2f}; pulls per step size in one run {runs[0][1]}")

# Combinatorial optimization: the TSP search space, then construction (L2C) vs. improvement (L2I) heuristics
import math
import random

for n in (5, 10, 14, 20):
    print(f"{n:2d} cities: (n - 1)! = {math.factorial(n - 1):.3e} tours, (n - 1)!/2 = {math.factorial(n - 1) // 2:.3e} undirected")
# 14 cities: (14 - 1)! = 13! = 6,227,020,800, not 622,702,080

rng = random.Random(5)
P = [(rng.random(), rng.random()) for _ in range(60)]
L = lambda t: sum(math.dist(P[t[i]], P[t[i - 1]]) for i in range(len(t)))


def nearest_neighbor(start=0):
    """Construction: build a tour by always visiting the closest unvisited city."""
    tour, left = [start], set(range(len(P))) - {start}
    while left:
        nxt = min(left, key=lambda j: math.dist(P[tour[-1]], P[j]))
        tour.append(nxt); left.remove(nxt)
    return tour


def two_opt(tour):
    """Improvement: reverse a segment whenever that shortens the tour, until no reversal helps."""
    improved = True
    while improved:
        improved = False
        for i in range(1, len(tour) - 1):
            for j in range(i + 1, len(tour)):
                a, b, c, d = P[tour[i - 1]], P[tour[i]], P[tour[j]], P[tour[(j + 1) % len(tour)]]
                if math.dist(a, c) + math.dist(b, d) < math.dist(a, b) + math.dist(c, d) - 1e-12:
                    tour[i:j + 1] = reversed(tour[i:j + 1]); improved = True
    return tour


t = nearest_neighbor()
print(f"nearest neighbor (construct): {L(t):.3f}")
print(f"after 2-opt (improve):        {L(two_opt(t)):.3f}")
rand = list(range(60)); rng.shuffle(rand)
print(f"random tour {L(rand):.3f} -> 2-opt {L(two_opt(rand)):.3f}")