What are we trying to do?

Examples in, a model out, new examples from the model

Our toy data

200 numbers from a recipe we keep hidden

The model sees only these 200 numbers, not the recipe. Goal: learn from them, then produce a 201st that looks like it belongs.

A tempting answer

The average is a terrible sample

The mean is a number that almost never occurs. Generating needs the whole spread of the data, not one summary of it.

Our first model: one Gaussian · how it scores a number

How probable is $x = 1.5$ under a bump at 2?

$$\N(x;\mu,\sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}\,\exp\!\Big(-\frac{(x-\mu)^2}{2\sigma^2}\Big)$$

$$\frac{x-\mu}{\sigma} = \frac{1.5 - 2}{0.5} = -1$$

$$\exp\!\big(-\tfrac12\cdot(-1)^2\big) = e^{-0.5} = 0.607$$

$$\frac{1}{0.5\sqrt{2\pi}} = 0.798 \qquad\Rightarrow\qquad 0.798\times0.607 = \mathbf{ 0.484 }$$

The model's two parameters, $\mu$ and $\sigma$, decide how high the curve is at every $x$. Higher means more probable.

Step 1 · Learn

Score a choice of $\mu, \sigma$ on all 200 numbers

$$\log p_\theta(\text{data}) = \sum_{i=1}^{200} \log \N(x_i;\mu,\sigma^2), \qquad \theta = (\mu, \sigma)$$

$\log p = −3746.0$

$\log p = −1790.0$

$\log p = −430.0$

Learning = keep the $\theta$ with the highest score. That is maximum likelihood, the same rule that trained the name model.

Step 1 · Learn, then Step 2 · Generate

The best Gaussian has a formula; then we sample from it

$$\hat\mu = \frac1N\sum_{i=1}^{N} x_i = −0.22 \qquad \hat\sigma^2 = \frac1N \sum_{i=1}^{N} (x_i-\hat\mu)^2 = 4.31$$

$$x_{\text{new}} = \hat\mu + \hat\sigma\,\eps,\quad \eps \sim \N(0,1)$$

Setting the derivative of the score to zero gives these two numbers: the same answer as the third candidate.

Step 2: 11 of 30 samples land in the empty gap. Learning worked; the model was too simple.

The same fit, in PyTorch

Learning = gradient descent on the negative log-likelihood

x = torch.tensor(data)                  # the 200 numbers
mu        = torch.tensor(1.0, requires_grad=True)
log_sigma = torch.tensor(0.0, requires_grad=True)
opt = torch.optim.Adam([mu, log_sigma], lr=0.05)

for step in range(300):
    dist = Normal(mu, log_sigma.exp())
    loss = -dist.log_prob(x).mean()
    opt.zero_grad(); loss.backward(); opt.step()

x_new = dist.sample((30,))              # step 2: generate

Watch the fit happen

Four moments of that loop

The curve starts in the wrong place, widens to cover both clumps, and stops at the best single bump. The loss falls from 3.82 to 2.15.

Checkpoint

So far: learn and generate work; one bump is not enough

A richer model: tell a story

Hidden choice first, then the number

This is how we will generate: run the two steps. The coin is never written down in the data, so it is a hidden (latent) variable $z$.

How the new model scores a number

To score $x$, add up every way it could have been made

$$p(x) = \underbrace{\tfrac12\,\N(x;-2,0.5^2)}_{\text{made by bump A}} + \underbrace{\tfrac12\,\N(x;\,2,0.5^2)}_{\text{made by bump B}} = \sum_{z} \cz{p(z)}\,\cdec{p(x\mid z)}$$

$$\text{at } x=1.5:\quad \tfrac12\cdot 0.000 + \tfrac12\cdot 0.484 = 0.242$$

One hidden coin flip made the model rich enough: the samples now fall in the bumps, not in the gap.

Fitting $K$ bumps in PyTorch

Same loop, a richer log_prob

logits = torch.zeros(K, requires_grad=True)   # bump weights
means  = torch.randn(K, requires_grad=True)
log_sd = torch.zeros(K, requires_grad=True)

def log_p(x):                                # log p(x), shape [N]
    comp = Normal(means, log_sd.exp()).log_prob(x[:, None])   # [N, K]
    return torch.logsumexp(log_softmax(logits) + comp, dim=1)

loss = -log_p(x).mean()        # then the same 3-line update

Step 2 · Generate, from the two-bump fit

Sampling runs the story: choose a bump, then a number

weights = softmax(logits)
k = torch.multinomial(weights, 1000,
                      replacement=True)    # the coin
x_new = means[k] + log_sd.exp()[k] * torch.randn(1000)

Step 1: pick a bump for each sample, with the learned weights.

Step 2: draw a number from that bump: its centre plus its width times standard noise.

1,000 new numbers that follow the data, and none of them is a copy of a training number.

But how many bumps should we fit?

More bumps always fit the training data better

1, 2, 3, 4 bumps, each fitted with the same loop. The training score rises every time. So which one do we keep?

Score on data the fit has never seen

New data decides: two bumps

Draw 200 new numbers from the same hidden recipe and score each fit on them.

1 → 2 bumps: $−2.153 \to −1.376$ per number. A real improvement.

3 or 4 bumps: the training score creeps up, the new-data score does not. The extra bumps fit noise in the 200 training numbers.

Choose the size of a model by how well it scores data it has not seen. That holds for bumps, and for networks.

Step 1 · Learn, looked at closely

Hidden in the fit: which bump made each point?

If we knew the coin for every point

fitting is easy: average the A-points for bump A's centre, the B-points for bump B's.

We don't, so ask for each point

given this $x$, how probable is it that bump B made it? That is $p(\mathrm{B}\mid x)$.

Running the story backwards, from a number to its likely hidden cause, is called inference. The gradient descent we just ran was doing it all along: four slides on, we will see where.

The backwards step, worked by hand

You see $x = 0.1$. Which bump made it?

Bayes' rule is just this table: multiply, then divide by the total. Note the total is exactly the $p(x)$ of the two-bump formula.

Do it for every $x$

The backwards step, as a curve

Repeat the table for every $x$ and plot $p(\mathrm{B}\mid x)$.

Far right: certainly B. Far left: certainly A. Near 0: unsure.

For these two bumps the curve is exactly $\mathrm{sigmoid}(16x)$: a logistic regression in disguise.

These are the colours on the 200 points two slides ago.

The backwards step, in one line

Posterior = softmax of (log prior + log likelihood)

log_joint = log_softmax(logits) + Normal(means, log_sd.exp()).log_prob(x)   # [K]
posterior = softmax(log_joint)                 # divide by the total: Bayes

Softmax is "exponentiate, then divide by the total": exactly the last column of the Bayes table.

Why learning needs the backwards step

The gradient of $\log p(x)$ is weighted by the posterior

$$\frac{\partial \log p(x)}{\partial \mu_k} \;=\; \hl{p(k\mid x)}\;\cdot\;\frac{x - \mu_k}{\sigma_k^2}$$

At $x = 0.1$: bump B's centre is pulled by $0.832 \times (0.1 - 2)/0.25 = −6.32$, bump A's by $0.168 \times (0.1 + 2)/0.25 = 1.41$.

Each point pulls each bump in proportion to how likely that bump made it. Learning needs the posterior for every point: a 2-row table here, an intractable patch for images.

Four words for the rest of the lecture

Prior, likelihood, evidence, posterior

Generating uses the prior and the likelihood. Learning needs the evidence and the posterior.

The plan, as a picture

Squeeze, fit bumps, sample, un-squeeze

Two questions before any code: each bump is a Gaussian over what? And how many bumps?

Which dimension is each Gaussian?

One bump, three possible sizes

each bump is a Gaussian over	centre	covariance	numbers per bump, with its weight
1 number (the toy)	1	1 (a width)	3
784 pixels, every pixel on its own	784	784 widths	1,569
784 pixels, with all pixel pairs	784	307,720	308,505
50 PCA numbers (what we fit)	50	1,275	1,326

Samples from one 784-pixel Gaussian with independent pixels: each pixel jitters on its own, so you get static around the average digit. Strokes need pixels that move together.

Keep the covariance, shrink the dimension: a full Gaussian over 50 numbers costs about 230 times less than one over 784 pixels.

The squeeze, and the way back

PCA and inverse PCA: generation needs both

PCA, image to $d$ numbers: coordinates along the $d$ directions in which digits vary most.

$$\underbrace{z}_{d} = \underbrace{V^{\top}}_{d\times 784}\,(\underbrace{x - \bar x}_{784})$$

Inverse PCA, $d$ numbers back to an image: add the directions back onto the mean digit.

$$\underbrace{\hat x}_{784} = \bar x + \underbrace{V}_{784\times d}\,\underbrace{z}_{d}$$

top: test digits, our 3 among them · bottom: PCA to 50 numbers, then inverse PCA

The bumps only ever see $z$. Every generated sample is $d$ numbers until the inverse PCA turns it into pixels. With $d = 50$, 82.5% of the variation survives the round trip.

Squeeze all the way to 2 numbers

In 2-D you can see every bump

4,000 test digits by their 2 PCA numbers, coloured by label (the model never saw labels). Black: 10 fitted bumps, 1 and 2 widths. Red ×: 24 samples.

Inverse PCA of a grid over the plane: every point is a blend of a few blurry shapes.

the 24 red × after inverse PCA

2 numbers keep only 16.8% of the variation. Every bump is visible, every sample is a ghost: interpretable is not the same as good.

The two-bump code, scaled up

PCA, fit $K$ bumps, sample, inverse PCA

x = mnist_train.flatten(1)                 # [60000, 784], grey levels
mean = x.mean(0)
_, _, Vt = torch.linalg.svd(x - mean, full_matrices=False)
V = Vt[:d].T                              # [784, d]
z = (x - mean) @ V                         # PCA: [60000, d]

gmm = MixtureSameFamily(Categorical(logits=logits),
                        MultivariateNormal(means, scale_tril=L))
loss = -gmm.log_prob(z_batch).mean()       # the same 3-line update

z_new = gmm.sample((24,))                  # step 2: generate
x_new = (mean + z_new @ V.T).clamp(0, 1)    # inverse PCA: [24, 784]

Try it: every combination below was fitted

Slide the dimension and the number of bumps

PCA numbers $d$50 bumps $K$100

variation kept by $d$ numbers: 82.5%
parameters: 132,599
score on new images: −16.5
new images prefer $K$ = 100 at this $d$

test digits

after PCA to $d$ numbers and back

24 samples from the $K$-bump fit, after inverse PCA

How many bumps for MNIST?

The toy's rule: score on images the model has not seen

$d = 50$. Fit $K$ = 1, 2, 5, 10, 20, 50, 100, 200, 400, 800 bumps to the 60,000 training images.

Training images: more bumps fit them better, as with the toy.

10,000 new test images: the score stops improving at $K = 100$. Doubling to 200 bumps changes the score by +0.2 nats. At 800 even the training score drops: that many bumps are more than 2,000 steps can fit.

Pick the smallest $K$ that does as well as any on new images: $K = 100$, 132,599 parameters, each bump a Gaussian over 50 numbers.

Choosing $K$ without held-out images

AIC and BIC: charge for every parameter

New images and AIC agree on $K = 100$; BIC, which charges more per parameter, keeps only 50. With plenty of data, measure on new data; with little, AIC or BIC.

What the bumps learned

Each bump's centre is a style of digit

the centres of the 100-bump fit, after inverse PCA, heaviest bump first (the heaviest has weight 3.9%)

The hidden coin now has 100 sides, one prototype each, found without the labels. The backwards step $p(z\mid x)$ is still just a 100-row table.

Step 2 · Generate

Samples, as the number of bumps grows

1 bump: a single Gaussian over the 50 numbers of every digit

10 bumps

100 bumps, the size chosen on new images

Or skip PCA: model the 784 pixels directly

A mixture of Bernoullis: each bump is 784 ink probabilities

$$\log p(x\mid k) = \sum_{i=1}^{784}\Big[x_i\log\pi_{ki} + (1 - x_i)\log(1 - \pi_{ki})\Big]$$

the 20 heaviest of 200 bumps: $\pi_k$, ink probability per pixel

samples: pick a bump, flip each pixel's coin

Inside a bump every pixel is on its own, so samples are speckled; only the choice of bump ties pixels together. Same lesson as the 784-pixel Gaussian.

An honest look

New digits, not copies, and not yet sharp

top: 8 samples · bottom: the closest of the 60,000 training images

What worked

The Part 1 loop, unchanged, learns a distribution over images and samples digit-like shapes, close to real digits but not copies of them.

What did not

Blur from the straight-line squeeze, and each bump can only stretch its template along fixed directions, not bend a stroke.

Why we cannot stop here

This recipe does not scale

We need a squeeze that can bend, learned by a network, and parts that are shared. Part 2: the autoencoder.

End of Part 1

What we did, and what comes next