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. I use the California housing dataset from chapter 2 of Hands-On Machine Learning. The task: predict a district’s median house value from its census stats.

This is inspired by Karpathy’s microgpt. 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.

The four pieces

It breaks into four pieces:

Autodiff. A Value class wraps a number and remembers what operation produced it, and from which inputs. It walks that graph backward to compute every gradient in one pass.
The model. A Neuron, a Layer, an MLP — Value objects wired together. Two layers, sixteen hidden units, one linear output.
The data. 20,640 California census districts, one row per district, downloaded once and cached locally.
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. Here it’s wrapped in a Value class:

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:

ACTIVATIONS = {
    "tanh": lambda act: act.tanh(),
    "relu": lambda act: act.relu(),
    "linear": lambda act: act,
}


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

    def __call__(self, x):
        act = sum((wi * xi for wi, xi in zip(self.w, x)), self.b)
        return ACTIVATIONS[self.g](act)

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

A neuron takes a weighted sum of its inputs, adds a bias, and passes the result through an activation function g. g names one of three options: tanh, relu, or linear, which just returns the sum untouched.

Then a layer:

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

    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()]

A layer is a list of neurons. Each one sees the same input x and produces its own output.

And finally, an MLP:

class MLP:
    def __init__(self, n_in, n_outs, g="tanh"):
        sizes = [n_in] + n_outs
        self.layers = [
            Layer(sizes[i], sizes[i + 1], g=(g if i != len(n_outs) - 1 else "linear"))
            for i in range(len(n_outs))
        ]

    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()]

An MLP chains layers together: each layer’s output becomes the next layer’s input.

MLP always forces the last layer’s g to "linear", no matter what g the hidden layers use. The output is a raw number, not squashed into [-1, 1] by tanh. That’s the difference between this and a classifier: house prices aren’t a yes/no.

The data

Download the CSV once, then read from the cached file on every run after that:

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 ML workflow in Hands-On ML. Each row is one California census block group, with:

eight numeric columns (location, age, room counts, income)
one categorical column, ocean_proximity
a target column, 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 instead, test-row information leaks into a statistic the model trains 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. Hands-On ML calls the ratios 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 go in, min-max normalized against training-set statistics. That’s the same normalize() function for any dataset. Nothing here is specific to housing.

Loss and evaluation

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

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 work well. Otherwise the network spends its first hundred iterations just learning that outputs should be six digits.

Training

Mini-batch gradient descent, four hundred iterations:

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 schedule. median_income alone is well-correlated with median_house_value here. This isn’t a weak-signal problem, so it doesn’t need aggressive tuning to converge.

Running it

End to end, from a clean cache:

$ 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. That’s well under the $95.6K you’d get by just guessing the training mean for every district. It’s a clean, monotonic curve, because the signal here is real and dense: every one of the 15,740 training examples pulls the gradient toward a useful answer.

Full code: gist.github.com/soasme/micromlp.py.