Autoregressive models: the chain rule, MADE, PixelCNN and transformer language models

Before this lesson: lesson 01 (maximum likelihood is a forward KL; a VAE gives only a bound on $\log p(x)$, here it is exact); cross-entropy from the information-theory track; convolutional networks and attention with masks from the deep-learning track. Time: about 70 minutes, plus the notebook.

Why this matters #

A VAE (lesson 02) gives only a lower bound on $\log p(x)$. An autoregressive (AR) model gives the exact value, and it trains with nothing more than cross-entropy. This one idea powers GPT-style language models (next-token prediction), WaveNet for raw audio, PixelCNN for images, and the token models that generate the codes of a VQ-VAE. In science, most learned weather emulators are autoregressive in time: GenCast samples the next state of the atmosphere given the previous two and feeds the result back in. The cost structure of this lesson (cheap training, sequential sampling) and its failure mode (errors that build up over a long rollout) apply to all of them.

1. Predict first #

2. The picture #

The chain rule with two bits #

Take two binary variables $x_1, x_2 \in \{0, 1\}$. Every joint distribution of them can be written as

So three numbers describe it completely: $p(x_1 = 1)$, $p(x_2 = 1 \mid x_1 = 0)$ and $p(x_2 = 1 \mid x_1 = 1)$. This is the chain rule of probability. There is no assumption yet: it holds for every distribution. For $D$ variables it reads

where $x_{<i} = (x_1, \dots, x_{i-1})$ is the prefix of $x_i$.

A tree is the picture. Draw one level per variable. Each edge carries a conditional probability. The probability of a full outcome is the product of the edges on its path from the root to its leaf.

Two directions of use.

• Evaluate $p(x)$ for an observed $x$. Every conditional is known, because its inputs $x_{<i}$ are part of the data. All $D$ factors can be computed in parallel.
• Sample a new $x$. Draw $x_1$, then $x_2$ given that $x_1$, and so on. Step $i$ needs the result of step $i - 1$, so this is sequential. It is called ancestral sampling.

Where the modeling happens. With $D$ binary variables, a full table of the joint needs $2^D - 1$ numbers: about $10^{308}$ for a tiny $32 \times 32$ binary image. An AR model replaces each conditional table by a neural network, for example $p(x_i = 1 \mid x_{<i}) = \sigma(\text{logit}_i(x_{<i}))$ with $\sigma$ the sigmoid. Three choices matter: the order of the variables, the family of each conditional (Bernoulli, a 256-way softmax over pixel values, a mixture of logistics, a Gaussian), and the network.

One network for all conditionals. To get all $D$ conditionals in one pass, a single network produces $D$ outputs, and output $i$ must never see $x_i, \dots, x_D$. Masks do this. MADE masks the weights of an MLP, PixelCNN masks the convolution kernels, and a transformer masks the attention matrix. Section 3 builds each one.

A 2D example: two moons by counting #

Discretize the plane into a $32 \times 32$ grid. Let $x_1$ be the column of a point and $x_2$ its row. The chain rule gives $p(x_1, x_2) = p(x_1)\,p(x_2 \mid x_1)$, and for discrete variables the maximum-likelihood estimate of each factor is a table of relative frequencies: count the points in each column, and in each cell of the column.

This model has $31 + 32 \cdot 31 = 1023$ free numbers: exactly as many as a full table of the 1024 cells, because with full tables the chain rule is not an approximation. That is fine for a small grid but impossible for an image. A neural AR model shares parameters across the conditionals; the counting model is the same idea without sharing.

3. Definitions and main results #

The autoregressive property #

Fix an order $x_1, \dots, x_D$. An autoregressive model is a set of conditionals $q_\theta(x_i \mid x_{<i})$, $i = 1, \dots, D$, such that:

1. for every value of $x_{<i}$, $q_\theta(\cdot \mid x_{<i})$ is a normalized distribution over $x_i$;
2. (autoregressive property) the parameters of factor $i$ are a function of $x_{<i}$ only.

The model is $q_\theta(x) = \prod_{i=1}^{D} q_\theta(x_i \mid x_{<i})$.

Theorem (an autoregressive model is normalized). For discrete $x \in \mathcal X^D$, $\sum_x q_\theta(x) = 1$. (For continuous $x$, replace the sums by integrals.)

Why this is true. Sum out the last variable first. Only the last factor contains it, and that factor sums to 1. What is left is the same model with one variable less.

Proof. Induction on $D$. For $D = 1$, $\sum_{x_1} q_\theta(x_1) = 1$ by condition 1.

1. In $\sum_{x_1} \cdots \sum_{x_D} \prod_{i=1}^{D} q_\theta(x_i \mid x_{<i})$, factor $i$ depends only on $x_1, \dots, x_i$ (condition 2). So only the last factor contains $x_D$.
2. Move $\sum_{x_D}$ next to the last factor: $\sum_{x_D} q_\theta(x_D \mid x_{<D}) = 1$ by condition 1.
3. What remains is the same sum for $D - 1$ variables, which is 1 by induction. $\blacksquare$

Why condition 2 matters. Take $D = 1$ and let the "conditional" see $x_1$ itself: $q(x_1 = 1 \mid x_1) = x_1$. Then $q(1) = 1$ and $q(0) = 1 - 0 = 1$. The training log-likelihood is perfect (0 nats), but $q(0) + q(1) = 2$: this is not a distribution. In the proof, step 2 fails: the last factor now depends on $x_D$ also through its parameters, so its sum over $x_D$ is no longer 1. A leak of $x_i$ into factor $i$ shows up as a training loss that is too good to be true.

Training: maximum likelihood with teacher forcing #

The negative log-likelihood of one data point is a sum of $D$ terms:

Each term is an ordinary cross-entropy loss: "predict $x_i$ from the true prefix $x_{<i}$". Feeding the true prefix from the data, not the model's own samples, is called teacher forcing. With masks, one forward pass gives all $D$ terms. Averaged over the data, the loss is the cross-entropy $H(p_{\text{data}}, q_\theta) = H(p_{\text{data}}) + \mathrm{KL}(p_{\text{data}} \,\|\, q_\theta)$, so training minimizes the forward KL, as in lesson 01. For language models, the mean loss per token $\bar\ell$ in nats gives the perplexity $e^{\bar\ell}$.

Cost. Evaluating $\log q_\theta(x)$ takes one parallel pass. Sampling takes $D$ sequential passes, because $x_i$ needs $x_{<i}$ first.

Exposure bias #

At training time, the prefix is always real data. At sampling time, it is the model's own output. A small error at step $i$ can move the prefix to a region the model never saw in training, where the later conditionals are worse. So errors can build up over a long sequence. This mismatch is called exposure bias.

Two precise remarks. First, exposure bias does not bias the likelihood: $\log q_\theta(x)$ is still exact. Second, a perfect model ($q_\theta = p_{\text{data}}$) has no mismatch, because its own prefixes follow the data distribution. It is a problem of imperfect models at generation time. Fixes: scheduled sampling (sometimes feed the model's own predictions during training), sequence-level training (for example RL fine-tuning), and, for time-stepping emulators, training on multi-step rollouts. For large language models, how much it matters in practice is debated.

MADE: masks from degrees #

MADE (masked autoencoder for distribution estimation; Germain et al. 2015) takes an MLP with $D$ inputs and $D$ outputs ($KD$ outputs for $K$ parameters per variable) and multiplies each weight matrix elementwise by a fixed binary mask. Every unit gets a degree:

• input $j$ has degree $j$;
• hidden unit $k$ of layer $l$ gets a degree $m_l(k) \in \{1, \dots, D - 1\}$. It means "this unit may depend on $x_1, \dots, x_{m_l(k)}$". A simple choice is $m_l(k) = 1 + (k \bmod (D - 1))$ for $k = 0, 1, \dots$; the paper samples the degrees at random;
• output $i$ (the parameters of factor $i$) has degree $i$.

The masks have the shape of the weight matrix (rows = units of this layer, columns = units of the layer below). With $m_0(j) = j$:

In words: a hidden unit may read units of the same or lower degree; an output may read only units of strictly lower degree.

Theorem (MADE has the autoregressive property). Output $i$ depends only on inputs $j < i$.

Proof.

1. Entry $(i, j)$ of $M^{\text{out}} M^L \cdots M^1$ counts the paths from input $j$ to output $i$ through allowed connections. If it is 0, output $i$ cannot depend on $x_j$, whatever the weights are.
2. Along any path $j \to k_1 \to \dots \to k_L \to i$, the hidden rule gives $j \le m_1(k_1) \le \dots \le m_L(k_L)$, and the output rule gives $m_L(k_L) < i$.
3. So every path has $j < i$. For $j \ge i$ there is no path: the matrix of path counts is strictly lower triangular. $\blacksquare$

Why the degrees are $1, \dots, D - 1$. A unit of degree $D$ could only feed outputs with $i > D$, which do not exist, so it is wasted. A unit of degree 0 has no inputs (it is a constant). Output 1 gets no connections at all: its parameters are just the biases, which is right for the unconditional factor $p(x_1)$.

Order-agnostic training (short). The MADE paper also trains one network for many orders at once: it samples a new order and new degrees for each minibatch. At test time, an ensemble over several orders gives a better likelihood. A model trained for one order works only for that order.

PixelCNN: masked convolutions #

For images, PixelCNN (van den Oord et al. 2016) uses the raster order: row by row, left to right. Pixel $(r, c)$ may depend on all rows above and on the pixels to its left in row $r$. A convolution is made causal by setting some kernel weights to zero. For $3 \times 3$ kernels (1 = weight kept; the center is the middle entry):

• Layer 1 uses mask A (center excluded). The feature at $(r, c)$ will predict input pixel $(r, c)$, so it must not see it.
• Later layers use mask B (center included). The layer-$(l-1)$ feature at $(r, c)$ already depends only on pixels before $(r, c)$, so reading it again is safe. Mask B in layer 1 would leak $x(r, c)$.

For RGB images the papers also order the channels R, G, B inside a pixel, and masks A and B also say whether a channel may see itself. Here we use one channel.

The blind spot. In a masked $3 \times 3$ layer, the feature at $(r, c)$ reads positions $(r + \Delta r, c + \Delta c)$ with either $\Delta r = -1$ and $\Delta c \in \{-1, 0, 1\}$ (the row above), or $\Delta r = 0$ and $\Delta c \le 0$ (same row, not to the right). After $L$ layers, input $(r', c')$ reaches output $(r, c)$ only as a sum of $L$ such offsets. Exactly $r - r'$ of them go up, and only those can move right, by at most 1 each. So

Pixel $(r - 1, c + 2)$ is never seen, although it comes before $(r, c)$ in raster order. The unseen region is a wedge to the upper right, and more layers do not remove it.

Transformer language models: causal masks and the KV cache #

A decoder-only transformer models tokens $x_1, \dots, x_T$ with the same chain rule. The input at positions $1, \dots, T$ is the sequence shifted right by one, $(\text{BOS}, x_1, \dots, x_{T-1})$, where BOS is a start token, and position $i$ outputs the distribution of $x_i$. The causal mask sets the attention weight from position $i$ to every later position $j > i$ to zero (the score is set to $-\infty$ before the softmax). So position $i$ sees $\text{BOS}, x_1, \dots, x_{i-1}$: exactly the prefix of $x_i$. The three architectures put the strict inequality in different places:

The KV cache. Under a causal mask, the hidden state at position $j$, in every layer, depends only on the inputs up to position $j$. So when a new token is appended, the keys and values of all earlier positions do not change. Decoding step $t$ computes the query, key and value of the new position only, appends the key and value to a cache, and attends over the cache. This is exact, not an approximation. Queries are not cached, because each query is used only once, at its own step.

Without a cache, generating $T = 1024$ tokens processes $\sum_{t=1}^{T} t = 524{,}800$ positions through the projections and MLPs, instead of 1024, and the attention work grows like $T^3$ instead of $T^2$. The price is memory: $2LTd$ numbers per sequence for standard multi-head attention ($L$ layers, model width $d$, one key and one value vector per layer and position); grouped-query attention stores fewer. A plain MADE has no such saving: it recomputes the full pass at each of its $D$ sampling steps.

Answers to section 1. (1) The network can copy $x_i$ into its prediction, so the training loss goes to 0 nats. But the result is not a distribution: in the $D = 1$ example above, the "probabilities" sum to 2. (2) $\log p(x)$ takes one pass, 1 ms. Sampling takes $32 \cdot 32 = 1024$ sequential passes, about 1 s.

4. Worked example #

(a) Two bits. $p(x_1 = 1) = 0.3$, $p(x_2 = 1 \mid x_1 = 0) = 0.8$, $p(x_2 = 1 \mid x_1 = 1) = 0.4$.

The sum is $0.14 + 0.56 + 0.18 + 0.12 = 1$, as the theorem says. The log-likelihood of $(0, 1)$ is $\log 0.7 + \log 0.8 = -0.3567 - 0.2231 = -0.5798$ nats.

The reverse order describes the same joint with three other numbers: $p(x_2 = 1) = 0.56 + 0.12 = 0.68$, $p(x_1 = 1 \mid x_2 = 1) = 0.12/0.68 = 0.1765$ and $p(x_1 = 1 \mid x_2 = 0) = 0.18/0.32 = 0.5625$. With full tables the order does not matter. With a restricted network for each conditional, it does.

(b) MADE masks, $D = 4$, one hidden layer of 5 units. Degrees $m(k) = 1 + (k \bmod 3)$ for $k = 0, \dots, 4$, so $m = (1, 2, 3, 1, 2)$. Row $k$ of $M^1$ has ones in the columns $j \le m(k)$. Row $i$ of $M^2$ has ones where $m(k) < i$:

Row $i$ of $M^2 M^1$ counts the paths from each input to output $i$. Check one entry: output 3 reads the hidden units with degree below 3, which are units 0, 1, 3, 4 (degrees 1, 2, 1, 2). All four read $x_1$, so there are 4 paths from $x_1$; units 1 and 4 also read $x_2$, so there are 2 paths from $x_2$. Output 1 sees nothing (only its bias), output 2 sees $x_1$, output 3 sees $x_1, x_2$, and output 4 sees $x_1, x_2, x_3$. The hidden unit of degree 3 feeds only output 4.

(c) The blind spot, by counting offsets. Stack masks A, B, B ($L = 3$) and look at output pixel $(4, 4)$, with rows and columns numbered from 0.

• Row 3 (one step up): $c' - 4 \le 1$, so columns up to 5.
• Row 2 (two steps up): columns up to 6. Row 1 (three steps up): columns up to 7.
• Row 4 itself: layer 1 (mask A) must move left, so only columns 1 to 3.

Pixels such as $(3, 6)$, $(2, 7)$ and $(1, 8)$ come before $(4, 4)$ in raster order, but no path reaches $(4, 4)$ from them. (Column 0 and row 0 are unseen for a different reason: three layers reach at most three steps away.)

5. In code #

NumPy checks the claims of sections 3 and 4: the joint table of the two bits, the MADE mask product, the normalization theorem on a tiny causal transformer (summed over all $3^3 = 27$ sequences), the causality that makes a KV cache exact, and the blind spot.

import numpy as np
from scipy.signal import convolve2d

# (a) Two bits: the chain rule gives the joint table
p1, p2 = 0.3, {0: 0.8, 1: 0.4}                 # p(x1=1) and p(x2=1 | x1)
joint = np.array([[(1 - p1) * (1 - p2[0]), (1 - p1) * p2[0]],
                  [p1 * (1 - p2[1]),       p1 * p2[1]]])
print(joint.round(2).tolist(), joint.sum().round(12), np.log(joint[0, 1]).round(4))
# [[0.14, 0.56], [0.18, 0.12]] 1.0 -0.5798

# (b) MADE masks from degrees (D = 4, one hidden layer of 5 units)
D, H = 4, 5
m_in = np.arange(1, D + 1)                     # input j has degree j
m_h = 1 + np.arange(H) % (D - 1)               # hidden degrees [1 2 3 1 2]
M1 = (m_h[:, None] >= m_in[None, :]).astype(int)   # hidden rule: >=
M2 = (m_in[:, None] > m_h[None, :]).astype(int)    # output rule: strict >
print(M2 @ M1)                                 # path counts: strictly lower triangular
# [[0 0 0 0]
#  [2 0 0 0]
#  [4 2 0 0]
#  [5 3 1 0]]

# (c) A tiny causal self-attention language model sums to 1 over all V^T sequences
rng = np.random.default_rng(0)
V, T, BOS, d = 3, 3, 3, 8                      # vocabulary {0,1,2}, length 3, BOS id 3
E, P = rng.normal(size=(V + 1, d)), rng.normal(size=(T, d))
Wq, Wk, Wv, Wo = (rng.normal(size=(d, d)) for _ in range(4))
Wout = rng.normal(size=(d, V))

def log_prob(seq, causal=True):
    inp = [BOS] + list(seq[:-1])               # shift right: position i predicts token i
    h = E[inp] + P
    S = (h @ Wq) @ (h @ Wk).T / np.sqrt(d)
    if causal:
        S = np.where(np.tril(np.ones((T, T))) == 1, S, -np.inf)   # no attention to j > i
    A = np.exp(S - S.max(1, keepdims=True)); A /= A.sum(1, keepdims=True)
    out = h + (A @ (h @ Wv)) @ Wo
    logits = out @ Wout
    logp = logits - np.log(np.exp(logits).sum(1, keepdims=True))
    return logp[np.arange(T), seq].sum(), out

seqs = [np.array(s) for s in np.ndindex(V, V, V)]                     # all 27 sequences
print(round(sum(np.exp(log_prob(s)[0]) for s in seqs), 12))              # 1.0
print(round(sum(np.exp(log_prob(s, causal=False)[0]) for s in seqs), 4)) # 0.9954: not 1
# change the second token: only the last position, which reads it, changes
_, a = log_prob(np.array([0, 1, 2]))
_, b = log_prob(np.array([0, 2, 2]))
print(np.abs(a - b).max(1).round(4))           # [0.     0.     4.3719]

# (d) PixelCNN receptive field: masks A, B, B on a 9 x 9 image, output pixel (4, 4)
A_ = np.array([[1, 1, 1], [1, 0, 0], [0, 0, 0]])
B_ = np.array([[1, 1, 1], [1, 1, 0], [0, 0, 0]])
g = np.zeros((9, 9)); g[4, 4] = 1               # backward pass from the output pixel
for K in (B_, B_, A_):                          # the last layer first
    g = convolve2d(g, K, mode="same")           # adjoint of a cross-correlation
for r in range(5):
    print("".join("o" if (r, c) == (4, 4) else "#" if g[r, c] > 0 else "." for c in range(9)))
# .........
# .#######.
# .######..
# .#####...
# .###o....

o is the output pixel and # marks the pixels it depends on. The empty wedge to the upper right of o is the blind spot. The leaky model in (c) sums to 0.9954 with this seed: a random leaky model can land on either side of 1, and a trained one goes far above it, because training rewards copying.

6. Where it shows up #

• Language models. GPT-style models are causal transformers trained with next-token cross-entropy, which is exactly the AR maximum-likelihood loss. Perplexity is $\exp$ of the mean negative log-likelihood per token. A later track on large language models will cover pretraining and inference.
• Audio and images. WaveNet uses dilated causal 1D convolutions over raw audio samples. PixelCNN++ (a mixture of logistics per pixel) reaches 2.92 bits per dimension on CIFAR-10.
• Two-stage generators. A VQ-VAE (lesson 02) turns an image into a grid of discrete codes. An AR prior over the codes (a PixelCNN in VQ-VAE, a transformer in VQGAN; DALL-E uses a transformer over the codes of a discrete VAE) generates new grids, and the decoder turns them into images.
• Normalizing flows. An AR model with Gaussian conditionals, $x_i = \mu_i(x_{<i}) + \sigma_i(x_{<i})\,u_i$ with $u_i \sim \mathcal N(0, 1)$, is a one-layer affine autoregressive flow. MAF stacks such layers with MADE inside (lesson 04).
• ML for science. Weather emulators step forward in time. GenCast samples each 12-hour step from $p(X^{t+1} \mid X^t, X^{t-1})$ with a diffusion model and is trained on single steps. GraphCast (deterministic) is fine-tuned on rollouts of 2 to 12 steps of 6 hours to limit error growth: the time-series form of exposure bias. For a long signal, such as a seismogram with $10^4$ samples, the $D$ sequential passes are the main reason to prefer flows or diffusion models, which update all samples in each network pass.

7. Common mistakes #

• Leaky masks. A non-strict MADE output mask, mask B in the first PixelCNN layer, or a transformer without the one-token shift. The training loss looks excellent, samples are garbage, and the probabilities over all states do not sum to 1. Check the mask product, or sum the model over all states of a tiny example.
• Thinking sampling is parallel because training is. Teacher forcing parallelizes only evaluation. Sampling stays sequential; the KV cache makes each step cheaper, not the number of steps smaller.
• Confusing exposure bias with a biased likelihood. The likelihood is exact. Exposure bias is about the distribution of the prefix at sampling time.
• Caching queries, or caching in a bidirectional model. Only past keys and values are reused. Without a causal mask, past states change when a token is added, so a cache would be wrong.
• Thinking more PixelCNN layers remove the blind spot. Every layer moves at most one column to the right per row up, so the wedge stays. The fix is a different architecture (two stacks).

8. Check yourself #

1. How many free parameters does a full probability table for $D = 10$ binary variables have? How many does a logistic AR model, $p(x_i = 1 \mid x_{<i}) = \sigma(w_i^\top x_{<i} + b_i)$, have?
   Answer

   $2^{10} - 1 = 1023$. Factor $i$ has $i - 1$ weights and one bias, so the AR model has $\sum_{i=1}^{10} i = 55$. For $D = 784$ (MNIST): $307{,}720$ against $2^{784} - 1$.
2. In MADE with $D = 5$, which hidden degrees are allowed? Why are 0 and 5 excluded?
   Answer

   $1, 2, 3, 4$. A degree-5 unit could only feed outputs with $i > 5$, which do not exist. A degree-0 unit sees no input, so it is a constant.
3. In worked example (b), which inputs does output 3 depend on, and through how many paths?
   Answer

   $x_1$ through 4 paths and $x_2$ through 2 paths: row 3 of $M^2 M^1$ is $(4, 2, 0, 0)$.
4. A stack of masked $3 \times 3$ layers (A, then B, B, ...) of any depth. Can output $(r, c)$ see pixel $(r - 2, c + 3)$? Pixel $(r - 3, c + 3)$?
   Answer

   $(r - 2, c + 3)$: never, because $c' - c = 3 > 2 = r - r'$ (blind spot). $(r - 3, c + 3)$: yes, with at least 3 layers (three up-right offsets; the row above is allowed in both masks).
5. A transformer LM generates $T$ tokens. How many sequential steps does it need with and without a KV cache, and what changes?
   Answer

   $T$ steps either way, because sampling is sequential. Without a cache, each step reprocesses the whole prefix (about $T^2/2$ positions in total). With a cache, each step processes one new position and attends over the $t$ cached keys and values.
6. Why does a model trained with teacher forcing behave differently at sampling time? Why does a perfect model not have this problem?
   Answer

   Training conditions on true prefixes, sampling on generated ones (exposure bias). If $q_\theta = p_{\text{data}}$, generated prefixes follow the data distribution, so there is no mismatch. The problem comes from errors of imperfect models that compound along the sequence.

9. Interview angle #

• "Why can a GPT train on all positions in parallel but must generate one token at a time?" — The causal mask gives every conditional from the true prefix in one pass (teacher forcing). Generation needs $x_{<i}$ before $x_i$.
• "How does MADE enforce the autoregressive property?" — Degrees, $\ge$ in the hidden masks, a strict output mask; the mask product is strictly lower triangular.
• "What is the PixelCNN blind spot, and how is it fixed?" — Each row up allows at most one column to the right, so a wedge to the upper right is never seen. The gated PixelCNN uses vertical and horizontal stacks.
• "Why is the KV cache exact? Why are queries not cached?" — Causality: past keys and values never change. Each query is used once.
• "What is exposure bias? Does it matter for LLMs?" — The train/sample prefix mismatch. Fixes: scheduled sampling, sequence-level training, multi-step rollout training. Its size for large LMs is debated.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Fit a grid AR model of two moons by counting; check that the probabilities sum to 1; sample from it.
2. Build MADE masks from degrees; check the dependency pattern of a random masked MLP with finite differences.
3. Exposure bias on a Markov chain: compare the error of an imperfect model under true prefixes and under its own prefixes.
4. Compute the receptive field of stacked PixelCNN masks and find the blind-spot pixels.
5. Write causal self-attention in NumPy and check that output $i$ does not change when later inputs change.

11. Watch and read #

• Murphy, Probabilistic Machine Learning: Advanced Topics, free PDF: chapter 22 (§22.1 the chain rule, §22.2 NADE and MADE, §22.3 PixelCNN, §22.4 transformers).
• Germain, Gregor, Murray and Larochelle, MADE: Masked Autoencoder for Distribution Estimation (2015): §4 (masks and the proof), §4.2 (order-agnostic training).
• van den Oord et al., Conditional Image Generation with PixelCNN Decoders (2016): §2.1 (gated units), §2.2 (the blind spot, vertical and horizontal stacks). Masks A and B with color channels are defined in the earlier Pixel Recurrent Neural Networks paper, §3.4.
• Stanford CS236, Deep Generative Models: course notes on autoregressive models.
• Karpathy, pytorch-made: a MADE in about 150 lines of PyTorch.