Project: Fit a Line with Gradient Descent

কী বানাবে

একটা flat ভাড়ার site বাড়িওয়ালাকে বলতে চায় "তোমার মতো flat মোটামুটি এত টাকায় ভাড়া হয়"। তোমার কাছে ঢাকার 12 টা flat এর size আর মাসিক ভাড়া আছে। তুমি একটা সোজা লাইন fit করবে, rent = w × size + b, আর করবে ঠিক সেভাবে যেভাবে সব neural network train হয়: loss লেখো, তার gradient লেখো, তারপর হাতে gradient descent চালাও। কোনো ML library ছাড়া।

পথে প্রথম যে আসল training bug টায় বেশিরভাগ মানুষ পড়ে, সেটাতেও পড়বে: এমন learning rate যেটা হয় খুব ধীর, নয়তো ফেটে যায়। সেটা ঠিক করবে feature scaling দিয়ে। শেষে exact formula র সাথে মিলিয়ে প্রমাণ করবে তোমার answer ঠিক।

যা আগে জানা লাগবে

• Calculus and Gradients: partial derivative আর gradient check।
• Optimization: update এর নিয়ম আর learning rate।
• Python 3 বা Go। কিছু install করতে হবে না।

Architecture

Training loop এর চারটা অংশ, আর এই course এ যত model train করবে সবগুলোতে এই চারটাই:

ধাপে ধাপে বানাও

Step 1: loss। সব flat এর উপর mean squared error: MSE = (w·x + b − y)² এর গড়।

Step 2: gradient। প্রতিটা term এ chain rule লাগাও। একটা flat এর ভুল e = w·x + b − y, আর:

∂(e²)/∂w = 2e × x        ∂(e²)/∂b = 2e × 1

সব flat এর উপর এগুলোর গড় নিলে পাও dw আর db।

Step 3: loop। w = b = 0 থেকে শুরু করে বারবার w -= lr * dw; b -= lr * db।

Step 4: input scale করো। Size গুলো 1,000 এর আশেপাশে, অথচ intercept গুণ হয় 1 দিয়ে, তাই w এর gradient b এর চেয়ে মোটামুটি হাজার গুণ বড়। একটা learning rate দুইটার জন্য মানায় না। সমাধান: size গুলোকে mean 0 আর spread 1 এ আনো, (x − mean) / std, train করো, তারপর weight গুলো আবার square feet এ ফিরিয়ে নাও।

Step 5: check। সোজা লাইনের জন্য একটা exact answer আছে (least squares)। তোমার gradient descent এর result এর সাথে মিলতে হবে।

"""Fit rent = w * size + b to 12 Dhaka flats with gradient descent written by hand."""
import math

SIZE = [650, 800, 900, 1000, 1100, 1200, 1250, 1400, 1500, 1650, 1800, 2000]  # square feet
RENT = [12, 15, 16, 19, 20, 22, 24, 25, 28, 30, 33, 36]                       # thousand taka a month
N = len(SIZE)

def mse(w, b, xs):
    return sum((w * x + b - y) ** 2 for x, y in zip(xs, RENT)) / N

def gradients(w, b, xs):
    dw = sum(2 * (w * x + b - y) * x for x, y in zip(xs, RENT)) / N
    db = sum(2 * (w * x + b - y) for x, y in zip(xs, RENT)) / N
    return dw, db

def train(xs, lr, steps, report=()):
    w = b = 0.0
    history = []
    for step in range(1, steps + 1):
        dw, db = gradients(w, b, xs)
        w, b = w - lr * dw, b - lr * db
        if step in report:
            history.append((step, mse(w, b, xs)))
    return w, b, history

# Step 1: raw square feet. Any learning rate is either too slow or explodes.
for lr, runs in ((1e-7, (10, 5000)), (1e-6, (10, 50))):
    for steps in runs:
        w, b, _ = train(SIZE, lr, steps)
        print(f"raw sizes, lr={lr:g}: after {steps} steps loss = {mse(w, b, SIZE):.4g}, b = {b:.4g}")

# Step 2: scale sizes to mean 0, spread 1. Now lr = 0.1 works.
mean = sum(SIZE) / N
std = math.sqrt(sum((x - mean) ** 2 for x in SIZE) / N)
scaled = [(x - mean) / std for x in SIZE]
w, b, history = train(scaled, 0.1, 100, report=(1, 5, 10, 20, 40, 60, 100))
top = history[0][1]
print("\nscaled sizes, lr=0.1")
for step, loss in history:
    print(f"step {step:>3}  loss {loss:8.3f}  {'#' * max(1, round(loss / top * 40))}")

# Step 3: turn the scaled weights back into rent per square foot.
per_sqft = w / std
base = b - w * mean / std
print(f"\nrent = {per_sqft * 1000:.2f} per 1000 sqft x size + {base:.2f}")
print(f"1300 sqft flat -> {per_sqft * 1300 + base:.1f} thousand taka")

# Step 4: check against the exact least-squares answer.
my = sum(RENT) / N
exact_w = sum((x - mean) * (y - my) for x, y in zip(SIZE, RENT)) / sum((x - mean) ** 2 for x in SIZE)
exact_b = my - exact_w * mean
print(f"exact:  {exact_w * 1000:.2f} per 1000 sqft, intercept {exact_b:.2f}")

// Fit rent = w * size + b to 12 Dhaka flats with gradient descent written by hand.
package main

import (
	"fmt"
	"math"
	"strings"
)

var size = []float64{650, 800, 900, 1000, 1100, 1200, 1250, 1400, 1500, 1650, 1800, 2000} // square feet
var rent = []float64{12, 15, 16, 19, 20, 22, 24, 25, 28, 30, 33, 36}                       // thousand taka a month
var n = float64(len(size))

func mse(w, b float64, xs []float64) float64 {
	s := 0.0
	for i, x := range xs {
		e := w*x + b - rent[i]
		s += e * e
	}
	return s / n
}

func gradients(w, b float64, xs []float64) (dw, db float64) {
	for i, x := range xs {
		e := w*x + b - rent[i]
		dw += 2 * e * x
		db += 2 * e
	}
	return dw / n, db / n
}

type point struct {
	step int
	loss float64
}

func train(xs []float64, lr float64, steps int, report map[int]bool) (w, b float64, history []point) {
	for step := 1; step <= steps; step++ {
		dw, db := gradients(w, b, xs)
		w, b = w-lr*dw, b-lr*db
		if report[step] {
			history = append(history, point{step, mse(w, b, xs)})
		}
	}
	return
}

func main() {
	// Step 1: raw square feet. Any learning rate is either too slow or explodes.
	for _, r := range []struct {
		lr   float64
		runs []int
	}{{1e-7, []int{10, 5000}}, {1e-6, []int{10, 50}}} {
		lr := r.lr
		for _, steps := range r.runs {
			w, b, _ := train(size, lr, steps, nil)
			fmt.Printf("raw sizes, lr=%g: after %d steps loss = %.4g, b = %.4g\n", lr, steps, mse(w, b, size), b)
		}
	}

	// Step 2: scale sizes to mean 0, spread 1. Now lr = 0.1 works.
	mean, varSum := 0.0, 0.0
	for _, x := range size {
		mean += x
	}
	mean /= n
	for _, x := range size {
		varSum += (x - mean) * (x - mean)
	}
	std := math.Sqrt(varSum / n)
	scaled := make([]float64, len(size))
	for i, x := range size {
		scaled[i] = (x - mean) / std
	}
	report := map[int]bool{1: true, 5: true, 10: true, 20: true, 40: true, 60: true, 100: true}
	w, b, history := train(scaled, 0.1, 100, report)
	top := history[0].loss
	fmt.Println("\nscaled sizes, lr=0.1")
	for _, p := range history {
		bars := int(math.Max(1, math.RoundToEven(p.loss/top*40)))
		fmt.Printf("step %3d  loss %8.3f  %s\n", p.step, p.loss, strings.Repeat("#", bars))
	}

	// Step 3: turn the scaled weights back into rent per square foot.
	perSqft := w / std
	base := b - w*mean/std
	fmt.Printf("\nrent = %.2f per 1000 sqft x size + %.2f\n", perSqft*1000, base)
	fmt.Printf("1300 sqft flat -> %.1f thousand taka\n", perSqft*1300+base)

	// Step 4: check against the exact least-squares answer.
	my, num := 0.0, 0.0
	for _, y := range rent {
		my += y
	}
	my /= n
	for i, x := range size {
		num += (x - mean) * (rent[i] - my)
	}
	exactW := num / varSum
	fmt.Printf("exact:  %.2f per 1000 sqft, intercept %.2f\n", exactW*1000, my-exactW*mean)
}

Output:

raw sizes, lr=1e-07: after 10 steps loss = 0.3622, b = 1.307e-05
raw sizes, lr=1e-07: after 5000 steps loss = 0.2677, b = 6.107e-05
raw sizes, lr=1e-06: after 10 steps loss = 7.603e+10, b = -0.1486
raw sizes, lr=1e-06: after 50 steps loss = 2.03e+43, b = -2.429e+15

scaled sizes, lr=0.1
step   1  loss  380.887  ########################################
step   5  loss   64.103  #######
step  10  loss    7.099  #
step  20  loss    0.321  #
step  40  loss    0.242  #
step  60  loss    0.242  #
step 100  loss    0.242  #

rent = 17.93 per 1000 sqft x size + 0.54
1300 sqft flat -> 23.9 thousand taka
exact:  17.93 per 1000 sqft, intercept 0.54

কী দেখলাম:

• Raw size, lr = 1e-7: loss দেখতে ঠিকই লাগে (0.27), কিন্তু 5,000 step পরেও b মাত্র 0.00006। w শিখেছে, b প্রায় নড়েইনি, কারণ তার gradient হাজার গুণ ছোট। Model ভুল intercept নিয়ে আটকে আছে, আর শুধু loss দেখে তুমি টেরও পেতে না।
• Raw size, lr = 1e-6: দশ গুণ বড় করতেই ফেটে গেল, 10 step এ 7.6 × 10¹⁰, 50 step এ 2 × 10⁴³। দুই knob এর জন্য একসাথে কাজ করে এমন কোনো learning rate নেই।
• Scaled size, lr = 0.1: loss 20 step এ 381 থেকে 0.32 এ নামে, তারপর 0.242 এ থিতু হয়। Bar গুলোই loss curve: শুরুতে খাড়া, তারপর সমতল।
• Answer: ভাড়া ≈ প্রতি 1,000 sqft এ 17.93 হাজার টাকা, সাথে 0.54। 1,300 sqft এর একটা flat এর ভাড়া হওয়া উচিত প্রায় 23.9 হাজার টাকা।
• Check: exact least-squares answer ও একই, 17.93 আর 0.54। তোমার loop আসল সেরা লাইনটাই পেয়েছে।

Test করো

import contextlib
import io

with contextlib.redirect_stdout(io.StringIO()):  # main.py prints its demo on import
    import main
from main import SIZE, gradients, mse, train

def test_gradient_matches_nudging():
    for w, b in [(0.0, 0.0), (0.01, 3.0), (-0.02, 10.0)]:
        dw, db = gradients(w, b, SIZE)
        h = 1e-6
        nw = (mse(w + h, b, SIZE) - mse(w - h, b, SIZE)) / (2 * h)
        nb = (mse(w, b + h, SIZE) - mse(w, b - h, SIZE)) / (2 * h)
        assert abs(dw - nw) / abs(nw) < 1e-4 and abs(db - nb) / abs(nb) < 1e-4

def test_loss_never_goes_up_with_a_good_learning_rate():
    _, _, hist = train(main.scaled, 0.1, 100, report=tuple(range(1, 101)))
    losses = [l for _, l in hist]
    assert all(b <= a + 1e-12 for a, b in zip(losses, losses[1:]))

def test_matches_exact_least_squares():
    assert abs(main.per_sqft - main.exact_w) < 1e-6
    assert abs(main.base - main.exact_b) < 1e-3

if __name__ == "__main__":
    for name, fn in list(globals().items()):
        if name.startswith("test_"):
            fn()
            print("pass", name)

pass test_gradient_matches_nudging
pass test_loss_never_goes_up_with_a_good_learning_rate
pass test_matches_exact_least_squares

• Gradient check: হাতে লেখা gradient কে কয়েকটা জায়গায় নাড়িয়ে মাপা মানের সাথে মিলতে হবে। একটা sign বা 2 এর factor ভুল করলে এটা সাথে সাথে fail করে, learning rate tune করে সময় নষ্ট করার আগেই।
• Loss কখনো বাড়ে না, এই convex problem এ ভালো learning rate দিলে। বাড়লে হয় learning rate বেশি, নয়তো gradient ভুল।
• Exact solution এর সাথে মিলে। শুধু সহজ model এ এটা সম্ভব, তাই যতদিন পারো কাজে লাগাও: পুরো loop টা শুরু থেকে শেষ পর্যন্ত ঠিক আছে, এটা প্রমাণ করে।

Trade-offs

আরো এগিয়ে নাও

1. দ্বিতীয় একটা feature যোগ করো, কত তলায় flat, তাহলে rent = w1·size + w2·floor + b। বাড়তি partial derivative টা লেখো আর gradient check দিয়ে মিলিয়ে নাও।
2. একটা outlier যোগ করো (1,000 sqft এর flat 90 হাজারে) আর দেখো লাইন কীভাবে হেলে যায়। তারপর Huber loss দিয়ে চেষ্টা করো।
3. প্রতি step এ random 4 টা flat এর mini-batch দিয়ে train করো, আর full batch এর loss curve এর সাথে তুলনা করো।
4. Matplotlib দিয়ে ঠিকঠাক loss curve আঁকো, loss এর axis log scale এ।