Ju Lin's AI Weblog

micromlp: A From-Scratch Neural Net That Predicts Housing Prices

I built micromlp: a single file of Python, no dependencies, no PyTorch. It downloads a real dataset, builds a 2-layer MLP, implements automatic differentiation from scratch, trains with gradient descent, and makes predictions. The task: the California housing dataset from chapter 2 of Hands-On Machine Learning — predict a district’s median house value from its census stats.

This is inspired by Karpathy’s microgpt. I’ve used PyTorch for years. I’d never actually written backprop by hand. There’s a difference between knowing .backward() exists and knowing what it does when you call it — I found it worth closing that gap once, on a toy, instead of taking it on faith forever.

The four pieces

Same shape as any neural net project, just shrunk down to where you can hold the whole thing in your head:

  1. Autodiff. A Value class that wraps a number, remembers what operation produced it and from which inputs, and can walk that graph backward to compute every gradient in one pass.
  2. The model. A Neuron, a Layer, an MLP — Value objects wired together. Two layers, sixteen hidden units, one linear output.
  3. The data. 20,640 California census districts, one row per district, downloaded once and cached locally.
  4. The task. Sixteen features per district, one target: median_house_value.

I already wrote about the autodiff part in more depth in Partial Derivatives, Reverse-Mode Autodiff — same idea, same six-line backward(). Here it’s wrapped in a Value class instead of loose functions, so it composes into a network:

class Value:
    __slots__ = ('data', 'grad', '_children', '_local_grads')

    def __init__(self, data, children=(), local_grads=()):
        self.data = data                # scalar value of this node, from the forward pass
        self.grad = 0.0                 # d(loss)/d(this node), filled in by the backward pass
        self._children = children       # nodes that feed into this one
        self._local_grads = local_grads # d(this node)/d(each child)

    def __add__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        return Value(self.data + other.data, (self, other), (1, 1))

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        return Value(self.data * other.data, (self, other), (other.data, self.data))

    def __pow__(self, other):
        return Value(self.data**other, (self,), (other * self.data**(other - 1),))

    def exp(self):
        e = math.exp(self.data)
        return Value(e, (self,), (e,))

    def log(self):
        return Value(math.log(self.data), (self,), (1 / self.data,))

    def tanh(self):
        t = math.tanh(self.data)
        return Value(t, (self,), (1 - t * t,))

    def relu(self):
        return Value(max(0, self.data), (self,), (float(self.data > 0),))

    def __neg__(self): return self * -1
    def __radd__(self, other): return self + other
    def __sub__(self, other): return self + (-other)
    def __rsub__(self, other): return other + (-self)
    def __rmul__(self, other): return self * other
    def __truediv__(self, other): return self * other**-1
    def __rtruediv__(self, other): return other * self**-1

    def backward(self):
        topo, visited = [], set()

        def build(v):
            if v not in visited:
                visited.add(v)
                for child in v._children:
                    build(child)
                topo.append(v)

        build(self)
        self.grad = 1.0
        for v in reversed(topo):
            for child, local_grad in zip(v._children, v._local_grads):
                child.grad += local_grad * v.grad

The model is Value objects wired into neurons:

class Neuron:
    def __init__(self, nin, nonlin=True):
        self.w = [Value(random.uniform(-1, 1) / nin**0.5) for _ in range(nin)]
        self.b = Value(0.0)
        self.nonlin = nonlin

    def __call__(self, x):
        act = sum((wi * xi for wi, xi in zip(self.w, x)), self.b)
        return act.tanh() if self.nonlin else act

    def parameters(self):
        return self.w + [self.b]


class Layer:
    def __init__(self, nin, nout, **kwargs):
        self.neurons = [Neuron(nin, **kwargs) for _ in range(nout)]

    def __call__(self, x):
        out = [n(x) for n in self.neurons]
        return out[0] if len(out) == 1 else out

    def parameters(self):
        return [p for n in self.neurons for p in n.parameters()]


class MLP:
    def __init__(self, nin, nouts):
        sizes = [nin] + nouts
        self.layers = [
            Layer(sizes[i], sizes[i + 1], nonlin=(i != len(nouts) - 1))
            for i in range(len(nouts))
        ]

    def __call__(self, x):
        for layer in self.layers:
            x = layer(x)
        return x

    def parameters(self):
        return [p for layer in self.layers for p in layer.parameters()]

The last layer is built with nonlin=False, so it skips the tanh — the output is a raw number, not squashed into [-1, 1]. That’s the difference between this and a classifier: house prices aren’t a yes/no.

That’s the whole engine. No tensors, no vectorization, one Python object per scalar. Slow, and that’s fine — the point was to see the machinery, not to compete with PyTorch.

The data

HOUSING_URL = "https://raw.githubusercontent.com/ageron/handson-ml2/master/datasets/housing/housing.csv"
CACHE_PATH = os.path.join(os.path.dirname(__file__), "housing.csv")


def fetch_housing(url=HOUSING_URL, cache_path=CACHE_PATH):
    if not os.path.exists(cache_path):
        with urllib.request.urlopen(url, timeout=20) as resp:
            raw = resp.read()
        with open(cache_path, "wb") as f:
            f.write(raw)
        print(f"downloaded housing.csv to {cache_path}")
    else:
        print(f"loaded housing.csv from cache: {cache_path}")

    with open(cache_path, newline="") as f:
        rows = list(csv.DictReader(f))
    print(f"{len(rows)} rows")
    return rows

This is the dataset Géron uses to teach the whole ML workflow in Hands-On ML: one row per California census block group, eight numeric columns (location, age, room counts, income) plus one categorical column (ocean_proximity), and a target — median_house_value. It has two well-known warts, both worth handling on purpose instead of by accident:

total_bedrooms is missing on about 1% of rows. I impute it with the median, computed from the training split only — if you compute it from the whole dataset, information about the test rows leaks into a statistic the model ends up trained against.

def clean_rows(rows):
    """total_bedrooms has ~1% missing values in this dataset. Hands-On ML's
    fix is to impute with the median, computed from training data only so
    the test set can't leak into a statistic the model is fit against."""
    bedrooms = [float(r["total_bedrooms"]) for r in rows if r["total_bedrooms"]]
    bedrooms.sort()
    median_bedrooms = bedrooms[len(bedrooms) // 2]
    for r in rows:
        if not r["total_bedrooms"]:
            r["total_bedrooms"] = median_bedrooms
    return median_bedrooms

And median_house_value is capped: every district actually worth more than $500,000 just got recorded as 500001.0. Training on that teaches the model a lie about the top of the market, so I drop those rows instead of pretending they’re real data points.

The features

Eight raw columns, five one-hot columns for ocean_proximity, and three ratios that Hands-On ML calls out specifically as more informative than the raw counts they’re built from — total_rooms alone doesn’t mean much without knowing how many households it’s spread across:

NUMERIC_COLUMNS = [
    "longitude", "latitude", "housing_median_age", "total_rooms",
    "total_bedrooms", "population", "households", "median_income",
]
OCEAN_CATEGORIES = ["<1H OCEAN", "INLAND", "ISLAND", "NEAR BAY", "NEAR OCEAN"]
ENGINEERED_COLUMNS = ["rooms_per_household", "bedrooms_per_room", "population_per_household"]
FEATURE_NAMES = NUMERIC_COLUMNS + ENGINEERED_COLUMNS + [f"ocean_{c}" for c in OCEAN_CATEGORIES]


def extract_features(row):
    total_rooms = float(row["total_rooms"])
    total_bedrooms = float(row["total_bedrooms"])
    households = float(row["households"])
    population = float(row["population"])

    features = {name: float(row[name]) for name in NUMERIC_COLUMNS}
    features["total_bedrooms"] = total_bedrooms
    features["rooms_per_household"] = total_rooms / households
    features["bedrooms_per_room"] = total_bedrooms / total_rooms
    features["population_per_household"] = population / households
    for cat in OCEAN_CATEGORIES:
        features[f"ocean_{cat}"] = float(row["ocean_proximity"] == cat)
    return features

Sixteen features in, min-max normalized against training-set statistics — same normalize() I’d use for any of these, nothing dataset-specific about it.

Loss and evaluation

This is regression, not classification, so the loss is mean squared error and the output is a raw number, not a probability:

PRICE_SCALE = 100_000.0  # train on median_house_value / 100k, so targets sit near 1-5


def mse_loss(pred, target):
    return (pred - target) ** 2


def rmse(model, X, y):
    se = sum((model(x).data - target) ** 2 for x, target in zip(X, y))
    return (se / len(y)) ** 0.5 * PRICE_SCALE

Scaling the target down by 100,000 before training isn’t cosmetic — it keeps the loss and the gradients in a range where tanh hidden units and a normal learning rate behave, instead of the network spending its first hundred iterations just learning that outputs should be six digits.

Training

def train(Xtrain, ytrain, Xtest, ytest, hidden=16, iters=400, batch_size=256, lr=0.05):
    random.seed(1)
    model = MLP(len(FEATURE_NAMES), [hidden, 1])
    for it in range(iters):
        batch_idx = [random.randrange(len(Xtrain)) for _ in range(batch_size)]
        loss = Value(0.0)
        for i in batch_idx:
            pred = model(Xtrain[i])
            loss = loss + mse_loss(pred, ytrain[i])
        loss = loss * (1.0 / batch_size)

        for p in model.parameters():
            p.grad = 0.0
        loss.backward()

        for p in model.parameters():
            p.data -= lr * p.grad

        if it % 50 == 0 or it == iters - 1:
            print(f"iter {it:4d}  train loss {loss.data:.4f}  test RMSE ${rmse(model, Xtest, ytest):,.0f}")
    return model

Plain mini-batch SGD, no momentum, no learning-rate schedule. median_income alone is famously well-correlated with median_house_value in this dataset, so unlike a weak-signal problem, this one doesn’t need aggressive learning rates or careful balancing to move — it just needs to run.

Running it

$ python3 micromlp.py
downloaded housing.csv to housing.csv
20640 rows

train examples: 15740, test examples: 3935
baseline (always predict the train-set mean): RMSE $95,605

iter    0  train loss 5.1530  test RMSE $173,384
iter   50  train loss 0.6881  test RMSE $79,907
iter  100  train loss 0.6123  test RMSE $75,151
iter  150  train loss 0.5459  test RMSE $71,276
iter  200  train loss 0.4978  test RMSE $68,384
iter  250  train loss 0.5186  test RMSE $65,821
iter  300  train loss 0.4303  test RMSE $63,916
iter  350  train loss 0.3888  test RMSE $62,943
iter  399  train loss 0.4760  test RMSE $62,591

sample predictions on held-out districts:
  predicted $169,687  actual $113,300
  predicted $215,223  actual $223,800
  predicted $195,557  actual $172,100
  predicted $125,612  actual $87,900
  predicted $193,666  actual $240,300
  predicted $204,689  actual $81,300
  predicted $154,692  actual $67,500
  predicted $220,351  actual $141,900

Test RMSE drops from $173K at initialization to $62.6K after 400 iterations, well under the $95.6K you’d get by just guessing the training mean for every district. That’s a clean, monotonic curve — a much easier story than a rare-event classifier, because the signal here is real and dense: every one of the 15,740 training examples pulls the gradient somewhere useful, instead of 96% of them agreeing on “predict no.”

The sample predictions still miss by a lot on individual districts — $204,689 predicted against $81,300 actual is a bad guess by any measure. A 2-layer, 16-hidden-unit network trained for 400 iterations on 16 features was never going to nail every outlier; it’s finding the broad shape of the relationship (income and location predict price), not the local noise.

What I’d do differently

Sixteen hidden units and 400 iterations is enough to beat the baseline meaningfully, not enough to compete with a tuned gradient-boosted tree, which is what actually wins this dataset in practice. Géron’s own book gets RMSE down near $50K with random forests. The gap between micromlp’s $62.6K and that isn’t a bug — it’s the cost of a from-scratch autodiff engine that processes one scalar at a time. A next version could add more capacity, more iterations, or feature crosses (latitude × longitude captures neighborhood effects that neither coordinate does alone) — I haven’t done any of that. The point of this project was the pipeline, not the leaderboard.

Full code: gist.github.com/soasme/micromlp.py. Delete housing.csv to re-download; the source data doesn’t change, so the numbers above should reproduce exactly.