intermediate language-modeling 45 min read

Neural Language Models: From Bengio's NPLM to Recurrent Networks

Replace counts with coordinates: a softmax over learned embeddings smooths by geometry, backpropagation is one reverse sweep of the chain rule, and a recurrent state lifts the trigram's horizon — as far as products of Jacobians let gradients reach

Part of the Language Modeling track · View full curriculum →

Prerequisites: Perplexity, Cross-Entropy, and Compression

Overview & Motivation

Take one held-out sentence from the finance corpus: cedar ridge bancorp owns brackenbury lending outright. Interpolated Kneser–Ney, the best smoother of the n-gram topic, spends $45.1$ bits on the sentence. For three of them, the pair each forms with the word before it never occurred in training, so Kneser–Ney can only back off to shorter contexts and hope. A small neural network that sees exactly the same two-word context spends $28.3$ bits. It has never seen those word pairs either. What it has is a position for every context in a learned space, and contexts near this one have told it what tends to come next.

This topic is about that move: from counts to coordinates. We build two neural language models from scratch, Bengio’s neural probabilistic language model (NPLM) and an Elman recurrent network, train them on the n-gram topic’s $69{,}312$ training tokens, and prove what each piece does.

Smoothing by geometry. The NPLM’s log-probabilities are Lipschitz in the context’s coordinates, so an unseen context is predicted like the seen contexts near it.
Training is maximum likelihood. The softmax is the categorical exponential family. At every stationary point the model’s average prediction equals the training frequencies, which has a sharp consequence: a word never seen as a target has no maximum-likelihood fit at all.
Backpropagation is reverse accumulation. One backward sweep over a computational graph gives every partial derivative at about twice the cost of the forward pass.
Recurrence and its limits. A recurrent state carries the whole sentence instead of two words. Its gradients through time are products of Jacobians, which we bound; on the trained network the bound is vacuous and the gradients still fade.

The perplexity topic made two promises about this one: neural models would be trained by minimizing its cross-entropy, and they would be compared in a unit that does not depend on their tokenizer. So every comparison here is an exact arithmetic code in bits per character, with unknown words spelled out under one speller for all models. The finance case study then takes the networks to public EDGAR filings, where adapting to the text being coded matters more than anything learned in training.

From counts to coordinates

Kneser–Ney treats a context as an atom. Two contexts share information only when one is a suffix of the other, through backoff. The NPLM instead places every context at a point, and lets nearness mean what the data says it means.

Notation follows the earlier topics. A sentence is padded with two <s> symbols in front, and every word and the final </s> is a predicted token. $V$ is the n-gram topic’s vocabulary, $\lvert V \rvert = 3{,}427$ , and an out-of-vocabulary word counts as <unk>. The held-out text has $2{,}049$ sentences, $23{,}055$ predicted tokens, and $154{,}362$ characters. Logarithms are natural unless written $\log_2$ .

Definition 1 (Neural probabilistic language model).

An embedding table $E$ has one row $E_w \in \mathbb{R}^{d}$ for each word, plus one for <s>. A context $u = (u_1, u_2)$ is the point $x(u) = [E_{u_1};\, E_{u_2}] \in \mathbb{R}^{2d}$ , the two embeddings stacked. The model computes

a = \tanh(x W_1 + b_1), \qquad z = a W_2 + b_2, \qquad p(w \mid u) = \operatorname{softmax}(z)_w = \frac{e^{z_w}}{\sum_{j \in V} e^{z_j}},

with row vectors throughout, $W_1 \in \mathbb{R}^{2d \times h}$ , and $W_2 \in \mathbb{R}^{h \times \lvert V \rvert}$ . Here $d = 64$ and $h = 128$ , for $677{,}987$ parameters (Bengio et al. 2003, without their direct input-to-output connections).

Proposition 1 (Positive and proper).

Every word has positive probability in every context, and each $p(\cdot \mid u)$ sums to one.

The exponential is positive and the softmax divides by the total, so this is immediate. It means the n-gram topic’s zero-count problem cannot occur. Corollary 1 below shows why positivity is not yet the same thing as smoothing.

Theorem 1 (Smoothing by geometry).

For all context points $x, x' \in \mathbb{R}^{2d}$ and every word $w$ ,

\bigl\lvert \log p(w \mid x) - \log p(w \mid x') \bigr\rvert \;\le\; \sqrt{2}\; \lVert W_2 \rVert_2\, \lVert W_1 \rVert_2\; \lVert x - x' \rVert_2 ,

where $\lVert \cdot \rVert_2$ is the Euclidean norm of a vector and the spectral norm of a matrix.

Proof.

Let $f(x) = \log p(w \mid x) = z_w - \log \sum_j e^{z_j}$ . Its derivative with respect to $z$ is $e_w - p$ , where $e_w$ is the indicator vector of $w$ . By the chain rule, in column form,

\nabla f(x) = W_1\, \mathrm{diag}(1 - a^2)\, W_2\, (e_w - p).

The entries of $\mathrm{diag}(1 - a^2)$ lie in $(0, 1]$ , so its norm is at most one. For the last factor, $\lVert e_w - p \rVert_2^2 = (1 - p_w)^2 + \sum_{j \ne w} p_j^2 \le (1 - p_w)^2 + \bigl(\sum_{j \ne w} p_j\bigr)^2 = 2(1 - p_w)^2 \le 2$ . Hence $\lVert \nabla f \rVert_2 \le \sqrt{2}\, \lVert W_1 \rVert_2 \lVert W_2 \rVert_2$ everywhere. Apply the mean value theorem (formalCalculus, Mean Value Theorem & Taylor) to $g(t) = f(x' + t(x - x'))$ on $[0, 1]$ : $f(x) - f(x') = \nabla f(\xi) \cdot (x - x')$ for a point $\xi$ on the segment, and the Cauchy–Schwarz inequality finishes.

∎

A context never seen in training is predicted like the seen contexts near it, and “near” is learned from every context at once: an embedding row moves whenever its word occurs anywhere (Corollary 2). For the trained NPLM the constant is $207.2$ , a worst case. Over about $2{,}000$ random pairs of distinct held-out contexts, the largest change in any word’s log-probability divided by the bound has median $0.0193$ and maximum $0.0403$ . The theorem holds with room to spare, and the room is the point: the constant certifies continuity, not how smooth the learned map actually is.

Where does geometry pay? Split the held-out tokens by what training saw:

held-out tokens	count	Kneser–Ney	NPLM	Elman RNN
bigram never seen in training	$5{,}112$	$7.525$	$6.735$	$6.247$
bigram seen, trigram not	$2{,}374$	$3.563$	$3.293$	$3.066$
trigram seen	$15{,}569$	$1.214$	$1.260$	$1.142$

The entries are nats per token. With the same two-word context as Kneser–Ney, the NPLM is $0.79$ nats better where the counts are absent and $0.046$ worse where they exist. Counts are the better estimate when you have them, and coordinates when you do not. The recurrent network of a later section wins on every row, because it sees more than two words.

Training is maximum likelihood

The NPLM is trained by minimizing its cross-entropy on the training text. By the perplexity topic’s Corollary 1 that is maximum likelihood under a strictly proper score. Here we look inside the objective, starting with its last layer.

Proposition 2 (The softmax is the categorical exponential family).

Let $A(z) = \log \sum_j e^{z_j}$ and $p = \operatorname{softmax}(z)$ . Then

\nabla A(z) = p, \qquad \nabla^2 A(z) = \frac{\partial p}{\partial z} = \mathrm{diag}(p) - p\, p^{\top} = \mathrm{Cov}_p(e_Y),

the covariance of the indicator vector $e_Y$ when $Y \sim p$ . This matrix is positive semidefinite, its null space is spanned by the all-ones vector $\mathbf{1}$ , and its rank is $K - 1$ for $K$ classes. Consequently the loss $\ell(z) = -\log p_y = A(z) - z_y$ is convex in $z$ , with gradient $p - e_y$ .

Proof.

Differentiating, $\partial p_i / \partial z_j = p_i(\delta_{ij} - p_j)$ , which is the stated matrix. For any vector $v$ , $v^{\top}\bigl(\mathrm{diag}(p) - p p^{\top}\bigr) v = \sum_i p_i v_i^2 - \bigl(\sum_i p_i v_i\bigr)^2 = \mathrm{Var}_p(v_Y) \ge 0$ . The variance is zero exactly when $v_Y$ is constant with probability one, which, since every $p_i \gt 0$ , means $v$ is a multiple of $\mathbf{1}$ . The gradient of $\ell$ is $\nabla A(z) - e_y = p - e_y$ .

∎

This is the general fact that a log-partition function’s gradient is the mean and its Hessian is the covariance (formalStatistics, Exponential Families, Theorem 7.1), for the categorical family. The geometry is worth a sentence: the gradient $p - e_y$ is the arrow from the observed vertex of the probability simplex to the model’s prediction, and descent moves the prediction toward the vertex.

Theorem 2 (Score identity and moment matching).

Fix features $a_1, \dots, a_n \in \mathbb{R}^{h}$ , for instance a network’s last hidden layer, with targets $y_1, \dots, y_n$ . Let $z_i = a_i W + b$ and $L(W, b) = \frac1n \sum_i \ell(z_i; y_i)$ . Then

\nabla_W L = \frac1n \sum_{i=1}^{n} a_i^{\top} (p_i - e_{y_i}), \qquad \nabla_b L = \frac1n \sum_{i=1}^{n} (p_i - e_{y_i}),

and $L$ is convex in $(W, b)$ . At any stationary point,

\frac1n \sum_{i=1}^{n} p_i = \frac1n \sum_{i=1}^{n} e_{y_i} \qquad \text{and} \qquad \sum_{i=1}^{n} a_i^{\top} p_i = \sum_{i=1}^{n} a_i^{\top} e_{y_i} :

the model’s average predicted distribution equals the empirical distribution of the targets, feature by feature. The first identity holds at every stationary point of a whole network, since the output bias is among its parameters.

Proof.

The logits are affine in $(W, b)$ , so the chain rule and Proposition 2 give both gradients. A convex function composed with an affine map is convex. At a stationary point both gradients vanish, and the two identities are those equations rearranged.

∎

This is formalStatistics’ canonical-link score identity $X^{\top}(y - \mu)$ (Generalized Linear Models, Theorem 1) for $K$ classes, the multinomial case that page defers. Its first identity has a consequence that matters for every language model.

Corollary 1 (A never-seen target has no maximum-likelihood fit).

If a word $w$ never occurs as a target, then at every parameter value

\frac{\partial L}{\partial b_w} = \frac1n \sum_{i=1}^{n} p_i(w) \;\gt\; 0 .

So $L$ has no stationary point, gradient descent lowers $b_w$ at every step, and the infimum of $L$ is approached only as $p(w \mid \cdot) \to 0$ .

Proof.

The $w$ -th coordinate of $\nabla_b L$ in Theorem 2 has no $-e_{y_i}$ term, because no target is $w$ , and each $p_i(w) \gt 0$ by Proposition 1. A function with a nonzero gradient everywhere has no stationary point, and each descent step changes $b_w$ by a negative multiple of a positive number.

∎

The n-gram topic’s zero-count problem returns, one level up. The softmax is positive, but fitting it by maximum likelihood is not smoothing, and the training text has a word that never occurs as a target: <unk>. Trained for one epoch on the raw text (the networks below train for up to four), the NPLM’s average probability for <unk> on held-out contexts is $2.9 \times 10^{-5}$ , below the uniform $1/\lvert V \rvert = 2.9 \times 10^{-4}$ . Each held-out unknown word then costs $9.76$ nats. On the EDGAR excerpts, where $57.6\%$ of tokens are unknown, the network scores $8.2058$ nats per token, worse than a uniform guess at $8.1394$ .

The cure is the n-gram topic’s Good–Turing logic: words seen once stand in for words not seen at all. In every epoch, each training token of a word seen once becomes <unk> with probability $\tfrac12$ , in the context and as a target alike. This is word dropout (Kiperwasser and Goldberg 2016), and Theorem 2 predicts its effect. Expected over the dropout, $\tfrac12 \cdot n_1 / N = 0.0135$ of the targets are <unk>. The trained NPLM’s average probability for <unk> over the training contexts is $0.0134$ , against $0.0140$ in one draw. Each held-out unknown word now costs $3.13$ nats, and on EDGAR the network beats the uniform guess ( $6.4165$ nats per token). That comparison is not matched in training length: this network trains for four epochs and the raw one for one. The matched comparison tells the same story: after its first epoch, the dropout network scores $5.475$ nats per token on EDGAR, against the raw network’s $8.206$ . Every network below is trained this way.

Proposition 3 (The network is not convex).

The NPLM’s training loss is not a convex function of its parameters.

Proof.

Relabeling the hidden units, by permuting the columns of $W_1$ and the entries of $b_1$ and the rows of $W_2$ the same way, leaves the network’s function unchanged, so the loss $L$ is unchanged: $L(\pi\theta) = L(\theta)$ . If $L$ were convex, its value at the midpoint would satisfy $L\bigl(\tfrac12(\theta + \pi\theta)\bigr) \le \tfrac12 L(\theta) + \tfrac12 L(\pi\theta) = L(\theta)$ . For the trained NPLM, with $\pi$ reversing the order of the hidden units, on the first $4{,}096$ training events (one draw of the word dropout below), $L(\theta) = L(\pi\theta) = 2.1004$ nats while the midpoint scores $3.1837$ .

∎

So the convergence rates of formalML’s Gradient Descent, which assume convexity, do not apply. Every network on this page is trained by Adam (Kingma and Ba 2015) from one random start, and nothing certifies that the result is a global optimum. Training stops at the epoch with the lowest loss on a separate dev text: a fresh $844$ -sentence run of the corpus generator, which also chooses the other tuned quantities below. The held-out text tunes nothing.

Backpropagation is reverse accumulation

Training needs the gradient of the loss with respect to $677{,}987$ parameters, at every step. Computing each partial derivative separately would cost a forward pass per parameter. Reverse accumulation gets all of them at once.

Definition 2 (Computational graph and adjoint).

A computational graph lists nodes $v_1, \dots, v_N$ in topological order. Each non-input node is $v_i = \varphi_i(v_{\mathrm{pa}(i)})$ for a differentiable primitive $\varphi_i$ of its parents, and $v_N$ is a scalar. Let $F_i(v_1, \dots, v_i)$ be $v_N$ computed by running the graph from node $i + 1$ on. The adjoint of $v_i$ is $\bar v_i = \partial F_i / \partial v_i$ : the derivative of the output with the earlier nodes held fixed and the later ones recomputed.

Theorem 3 (Reverse accumulation).

The adjoints satisfy $\bar v_N = 1$ and, for $i \lt N$ ,

\bar v_i = \sum_{c \,:\, i \in \mathrm{pa}(c)} \bar v_c \, \frac{\partial \varphi_c}{\partial v_i},

so one sweep in reverse topological order computes every adjoint. 2. If each primitive’s vector–Jacobian products, one per parent, cost at most $\omega$ times evaluating the primitive, the sweep costs at most $\omega$ times the forward pass, plus one addition per edge and coordinate. For a matrix product $Y = XW$ , the products $\bar Y W^{\top}$ and $X^{\top} \bar Y$ each cost exactly as many multiply–adds as $XW$ , so $\omega = 2$ .

Proof.

(1) By construction $F_{i-1}(v_1, \dots, v_{i-1}) = F_i\bigl(v_1, \dots, v_{i-1}, \varphi_i(v_{\mathrm{pa}(i)})\bigr)$ . The multivariate chain rule (formalCalculus, The Jacobian & Multivariate Chain Rule, Theorem 2) gives, for every $j \lt i$ ,

\frac{\partial F_{i-1}}{\partial v_j} = \frac{\partial F_i}{\partial v_j} + \frac{\partial F_i}{\partial v_i}\,\frac{\partial \varphi_i}{\partial v_j},

where the last term is zero unless $j$ is a parent of $i$ . The sweep keeps one accumulator per node, starts with $\bar v_N = 1$ and zeros elsewhere, and on reaching node $i$ adds $\bar v_i\, \partial \varphi_i / \partial v_j$ to each parent’s accumulator, which is this equation. By induction on the number of nodes processed, once nodes $N, \dots, i + 1$ are done, the accumulator of every remaining node $v_j$ , $j \le i$ , holds $\partial F_i / \partial v_j$ . This holds at the start, since $F_N = v_N$ , and processing node $i$ turns each $\partial F_i / \partial v_j$ into $\partial F_{i-1} / \partial v_j$ by the equation above (non-parents gain zero). In particular the accumulator of $v_i$ holds $\partial F_i / \partial v_i = \bar v_i$ , and it is final because no unprocessed node depends on $v_i$ . (2) Each node’s products are formed once and each edge’s contribution is added once.

∎

Reverse mode goes back to Linnainmaa (1976), the cost bound to Baur and Strassen (1983), and its use in networks to Rumelhart, Hinton, and Williams (1986); Griewank and Walther (2008) treat it in general. A node with several children sums their contributions, which is where the graph’s fan-out shows up: a word’s embedding that appears in two positions of a context gets the sum of both positions’ gradients.

Remark (Checking a gradient exactly: the complex step).

A hand-written gradient needs a reference. Finite differences subtract two nearly equal numbers and lose half their digits. For a real-analytic function, $f(\theta + i h e_j) = f(\theta) + i h\, \partial_j f(\theta) + O(h^2)$ , so $\operatorname{Im} f(\theta + i h e_j)/h$ is the derivative with no subtraction at all (Squire and Trapp 1998). With $h = 10^{-30}$ it is exact to rounding. Every hand-written gradient on this page, the NPLM’s, the recurrent network’s, and the long-lag Elman and LSTM networks’ below, is checked against it on random instances, with central differences as a second check. The companion notebook runs Theorem 3 itself on $200$ random graphs, $198$ of them with fan-out (matrix products, $\tanh$ , products and sums of vectors, log-sum-exp), and every input’s adjoint matches its complex-step derivative to $10^{-11}$ .

Corollary 2 (The NPLM's gradient).

For one event with target $y$ , write $g_z = p - e_y$ . Reverse accumulation through Definition 1 gives

\frac{\partial \ell}{\partial W_2} = a^{\top} g_z, \quad \frac{\partial \ell}{\partial b_2} = g_z, \quad g_s = (g_z W_2^{\top}) \odot (1 - a^2), \quad \frac{\partial \ell}{\partial W_1} = x^{\top} g_s, \quad \frac{\partial \ell}{\partial b_1} = g_s,

and the embedding row of each context word receives its slice of $g_s W_1^{\top}$ . Rows of words absent from a batch receive nothing.

So a word’s embedding moves whenever the word appears, in whatever context, by an amount that depends on what followed. That is the mechanism by which Theorem 1’s geometry is learned. Whether words used alike end up near each other is the subject of the later topic on word embeddings. We do not prove it here.

Recurrent networks and gradients through time

The NPLM sees two words, like the trigram it replaces. Take pellucid software based in dublin expanded its headcount clinicbase footprint. The last word is predictable from expanded its, four words back, and from nothing in its two-word context. The NPLM spends $57.8$ bits on the sentence. A recurrent network, whose state summarizes everything read so far, spends $41.3$ .

Definition 3 (Elman recurrent language model).

With $h_0 = 0$ and $x_t$ the embedding of the previous word (<s> first),

h_t = \tanh(x_t W_x + h_{t-1} W_h + b), \qquad p(\cdot \mid \text{prefix}) = \operatorname{softmax}(h_t W_o + b_o).

Here the embeddings have $64$ dimensions and the state $128$ , for $686{,}179$ parameters (Elman 1990; Mikolov et al. 2010). The state is reset at every sentence, so held-out sentences stay independent as the perplexity topic’s standard errors assume.

The state has no fixed horizon, and the model is not Markov of any order in the words. Does it still define a probability distribution over sentences, or could probability leak to sentences that never end?

Theorem 4 (A softmax recurrent model is tight).

Write $W_o^{(j)}$ for the column of $W_o$ belonging to word $j$ , and “eos” for the index of </s>. Let

\varepsilon = \Bigl( \sum_{j} \exp\bigl( \lVert W_o^{(j)} - W_o^{(\mathrm{eos})} \rVert_1 + b_{o,j} - b_{o,\mathrm{eos}} \bigr) \Bigr)^{-1} .

Then the probability of ending, $p(\mathrm{eos} \mid h)$ , is at least $\varepsilon$ for every state with $\lVert h \rVert_\infty \le 1$ , and so for every state the network can reach. Consequently, for the length $N$ of a sentence sampled from the model, $P(N \gt n) \le (1 - \varepsilon)^n$ . The probabilities of all finite sentences sum to one, and every moment of $N$ is finite.

Proof.

By Hölder’s inequality, $z_j - z_{\mathrm{eos}} = h \cdot (W_o^{(j)} - W_o^{(\mathrm{eos})}) + b_{o,j} - b_{o,\mathrm{eos}} \le \lVert h \rVert_\infty \lVert W_o^{(j)} - W_o^{(\mathrm{eos})} \rVert_1 + b_{o,j} - b_{o,\mathrm{eos}}$ , and $p(\mathrm{eos} \mid h) = 1 / \sum_j e^{z_j - z_{\mathrm{eos}}}$ , which gives the first claim since $\tanh$ keeps every coordinate in $(-1, 1)$ . The probability that none of the first $n$ tokens ends the sentence is $\mathbb{E}\bigl[\prod_{t \le n} (1 - p(\mathrm{eos} \mid h_t))\bigr] \le (1 - \varepsilon)^n$ . So $P(N = \infty) = \lim_n P(N \gt n) = 0$ , and $\mathbb{E}[N^k] \le \sum_n (n + 1)^k (1 - \varepsilon)^n \lt \infty$ .

∎

A model with this property is called tight (Du et al. 2023). Unbounded states break the argument, and for recurrent networks with unbounded activations even deciding consistency is hard (Chen et al. 2018). Tightness is what lets the perplexity topic treat a model’s sentences as a sentence source with finite moments. For the trained network $\varepsilon = 7.1 \times 10^{-17}$ : positive, so the model is tight, but as a numerical tail bound useless. Sampled sentences end quickly anyway, and the smallest probability of ending along $500$ held-out sentences is far above $\varepsilon$ .

Corollary 3 (Backpropagation through time).

Unroll the recurrence over a sentence and apply Theorem 3. With $\delta_t$ the gradient at the pre-activation of step $t$ and $g_t$ the gradient arriving at $h_t$ from the loss at step $t$ itself,

\delta_t = \bigl(g_t + \delta_{t+1} W_h^{\top}\bigr) \odot (1 - h_t^2), \qquad \frac{\partial L}{\partial W_h} = \sum_{t} h_{t-1}^{\top} \delta_t ,

and likewise for $W_x$ , $b$ , and the embeddings: weights shared across steps receive the sum of their per-step gradients (Werbos 1990).

The gradient that reaches an early state from a late loss passes through every step in between, and each step multiplies it by a Jacobian.

Theorem 5 (Products of Jacobians).

For the recurrence of Definition 3,

\frac{\partial h_t}{\partial h_k} = \prod_{j=t}^{k+1} \mathrm{diag}(1 - h_j^2)\, W_h^{\top},

with the factor for $j = t$ on the left, and therefore

\Bigl\lVert \frac{\partial h_t}{\partial h_k} \Bigr\rVert_2 \;\le\; \prod_{j=k+1}^{t} \bigl\lVert \mathrm{diag}(1 - h_j^2)\, W_h^{\top} \bigr\rVert_2 \;\le\; \sigma_{\max}(W_h)^{\,t-k} .

If $\sigma_{\max}(W_h) \lt 1$ , every contribution of a loss at step $t$ to the state at step $k$ shrinks at least geometrically in $t - k$ .
Growth with $t - k$ is possible only if $\sigma_{\max}(W_h) \gt 1$ .
Without the $\tanh$ , $\partial h_t / \partial h_k = (W_h^{\top})^{t-k}$ , and $\rho(W_h)^m \le \lVert W_h^m \rVert_2$ with $\lVert W_h^m \rVert_2^{1/m} \to \rho(W_h)$ , the spectral radius. So linear gradients explode when $\rho(W_h) \gt 1$ and vanish when $\rho(W_h) \lt 1$ .

Proof.

The product is the chain rule applied to $h_t = \tanh(x_t W_x + h_{t-1} W_h + b)$ , whose Jacobian in $h_{t-1}$ is $\mathrm{diag}(1 - h_t^2)\, W_h^{\top}$ . The first inequality is submultiplicativity of the operator norm, and the second uses $\lVert \mathrm{diag}(1 - h_j^2) \rVert_2 \le 1$ and $\lVert W_h^{\top} \rVert_2 = \sigma_{\max}(W_h)$ . Statements 1 and 2 read the bound in each direction. For 3, an eigenvector $W v = \lambda v$ gives $\lVert W^m v \rVert = \lvert \lambda \rvert^m \lVert v \rVert$ , so $\lVert W^m \rVert_2 \ge \rho(W)^m$ . The limit is Gelfand’s formula (Horn and Johnson 2012, Corollary 5.6.14).

∎

This is the vanishing- and exploding-gradient analysis of Bengio, Simard, and Frasconi (1994) and Pascanu, Mikolov, and Bengio (2013), whose gradient clipping we use during training. The theorem gives conditions on $W_h$ alone. The trained network satisfies none of the clean ones: $\sigma_{\max}(W_h) = 3.350$ and $\rho(W_h) = 2.131$ , so neither statement 1 nor the linear picture predicts anything. Worse for the bound, every one of its per-step factors along $200$ held-out sentences exceeds one, with a minimum of $1.62$ and a median of $2.06$ . Measured back from each sentence’s last state:

lag $m$	$1$	$2$	$4$	$8$	$10$	$12$	$14$
median $\lVert \partial h_T / \partial h_{T-m} \rVert_2$	$2.024$	$2.223$	$1.899$	$1.026$	$0.704$	$0.449$	$0.219$
median bound from per-step norms	$2.0$	$4.1$	$17.4$	$296$	$1{,}241$	$4{,}786$	$17{,}656$

The bound grows by a factor of about two per step and the measured norm falls. Saturation explains part of the gap ( $37.9\%$ of the units sit beyond $\lvert h \rvert = 0.9$ , where $1 - h^2$ is small), but not all of it, since every factor still has norm above one. The rest is alignment: the direction one factor stretches is not the direction the next one stretches, so the product shrinks while every factor expands. The measured norm peaks at lag $2$ , falls at every lag after that, and drops below one between lags $8$ and $9$ ( $1.026$ , then $0.871$ ). No bound built from per-step norms can see this. In effect the recurrent state lifts the trigram’s two-word horizon, softly, to about a dozen words.

The long short-term memory network was designed to keep a gradient alive across many steps (Hochreiter and Schmidhuber 1997).

Definition 4 (LSTM).

Gates $i_t, f_t, o_t$ (sigmoids) and $g_t$ ( $\tanh$ ) are functions of $(x_t, h_{t-1})$ , and

c_t = f_t \odot c_{t-1} + i_t \odot g_t, \qquad h_t = o_t \odot \tanh(c_t),

where $c_t$ is the cell and $f_t$ the forget gate (Gers, Schmidhuber, and Cummins 2000).

Proposition 4 (The cell path).

With $h_k = o_k \odot \tanh(c_k)$ and $o_k$ held fixed,

\frac{\partial c_t}{\partial c_k} = \prod_{j=k+1}^{t} \mathrm{diag}(f_j) + R_{t,k},

where every term of $R_{t,k}$ contains at least one factor $\partial c_j / \partial h_{j-1}$ , which is linear in the recurrent weights of the gates $i$ , $f$ , and $g$ . In particular, if those weights are zero, $\partial c_t / \partial c_k = \prod_j \mathrm{diag}(f_j)$ exactly: a product with no weight matrix in it.

Proof.

Treat the state as the pair $(c_j, h_j)$ and expand $\partial c_t / \partial c_k$ , by the chain rule, as a sum over paths through the one-step block Jacobians. No gate reads the cell, so $\partial c_j / \partial c_{j-1} = \mathrm{diag}(f_j)$ , and the path that never leaves the cell contributes $\prod_j \mathrm{diag}(f_j)$ . Every other path enters the hidden state and must come back to the cell, since the output is $c_t$ , and the only block leading from the hidden state back to the cell is

\frac{\partial c_j}{\partial h_{j-1}} = \mathrm{diag}\bigl(c_{j-1} \odot f_j (1 - f_j)\bigr) W_{h,f}^{\top} + \mathrm{diag}\bigl(g_j \odot i_j (1 - i_j)\bigr) W_{h,i}^{\top} + \mathrm{diag}\bigl(i_j \odot (1 - g_j^2)\bigr) W_{h,g}^{\top},

which is linear in the three gates’ recurrent weights.

∎

With no forget gate, as in the original LSTM, the cell path is the identity, Hochreiter and Schmidhuber’s “constant error carousel”. With one, the forget gate’s bias $b_f$ sets the cell path at initialization to about $\sigma(b_f)^{t-k}$ . A bias of one (Jozefowicz, Zaremba, and Sutskever 2015) gives $\sigma(1)^{60} = 6.9 \times 10^{-9}$ across sixty steps, and a bias of three gives $0.054$ . To see what that decides, train three networks with $16$ hidden units to read a bit, read $60$ random noise symbols, and then report the bit:

model	median gradient reaching the bit at initialization (12 starts)	accuracy after 600 training steps (one start)
Elman network	$2.45 \times 10^{-10}$	$0.486$
LSTM, forget bias $1$	$2.42 \times 10^{-7}$	$0.531$
LSTM, forget bias $3$	$\mathbf{1.42 \times 10^{-3}}$	$\mathbf{1.000}$

In this run, only the network whose gradient reaches the bit learns to report it. The direction follows Proposition 4, and the scale does not: the measured ratio between the two LSTMs, about $5{,}900$ , is far smaller than the cell paths’ ratio of about $7.9 \times 10^{6}$ , because the gradient also reaches the bit through the input gate and through paths that leave the cell. Each accuracy in the table comes from one training start. In an exploratory rerun from eight further starts per model, which the harness does not run, the bias-three network learned the task from seven, the bias-one network from none, and the Elman network from two. A gradient that reaches the bit at initialization makes learning likely, not certain, and a vanishing one makes it unlikely, not impossible. Choosing $b_f$ from the lag a task needs is the idea behind chrono initialization (Tallec and Ollivier 2018).

The softmax bottleneck

Arrange a model’s predictions as a matrix: one row per context, one column per word, entries $\log p(w \mid u)$ . For interpolated Kneser–Ney over the $1{,}000$ most frequent training contexts, that matrix has rank $611$ . For the NPLM it has rank $130$ .

Theorem 6 (The softmax bottleneck).

Suppose a model’s logits are $z(u) = s(u) W + b$ for a state $s(u) \in \mathbb{R}^{d}$ . On any $N$ contexts, let $S$ stack the states, and let $\lambda_i$ be the log-sum-exp of row $i$ of the logits. The log-probability matrix is

L = S W + \mathbf{1} b^{\top} - \lambda \mathbf{1}^{\top},

so $\operatorname{rank} L \le d + 2$ . Consequently, a target matrix $T$ of log-probabilities can be represented exactly by such a model only if adding a suitable constant to each row of $T$ brings its rank down to at most $d + 1$ , however the model is trained.

Proof.

$S W$ has rank at most $d$ , and each of the two outer products has rank at most one. Adding $\lambda_i$ back to row $i$ of $L$ removes the term $\lambda \mathbf{1}^{\top}$ and leaves $S W + \mathbf{1} b^{\top}$ , of rank at most $d + 1$ . So if $L = T$ , the row-shifted matrix $T + \lambda \mathbf{1}^{\top}$ has rank at most $d + 1$ .

∎

For the NPLM $d = h = 128$ , and its matrix has rank exactly $130 = d + 2$ , as does the recurrent network’s over $1{,}000$ held-out prefixes. A mixture of several softmaxes is not of this form, and Yang et al. (2018) proposed it for exactly this reason. Whether the true log-probability matrix of a language has high rank is their hypothesis, not a theorem. On our corpus the bottleneck is real in rank and small in energy: the best rank- $129$ approximation of Kneser–Ney’s row-centered matrix leaves only $1.48\%$ of its squared Frobenius norm (Eckart–Young, formalML’s Singular Value Decomposition). Frobenius distance between log-probabilities is not cross-entropy, so this does not bound what the bottleneck costs. A low-rank matrix of log-probabilities, or of association scores like pointwise mutual information, is also what word embeddings factorize, which is where the later topic on word2vec begins.

One sentence, four models; gradients through time

The laboratory puts the last four sections side by side. In panel A, pick a sentence: each token gets one bar per model, in bits, and the chip beside the word says what training saw. In the first sentence the networks win on the tokens whose word pair is new. In the second, the recurrent network wins on footprint, three words after expanded its. In the EDGAR sentence, an unknown word’s hatched spelling segment is the same for every model, and dynamic evaluation, which keeps learning on the text it is coding, takes over. Click a token, and panel B shows how strongly its loss depends on every earlier state of the trained network, against Theorem 5’s bound. Switch panel B to the synthetic view to vary the spectral norm of $W_h$ , with and without the $\tanh$ .

A · One sentence, four models — bits per token, and what training saw

bigram never seenbigram seen, trigram nottrigram seen in traininggold words are unknown and spelled letter by letter (hatched)

Interpolated KN, whole sentence

45.1 bits

NPLM, whole sentence

28.3 bits

Elman RNN, whole sentence

28.3 bits

RNN, dynamic evaluation, whole sentence

29.6 bits

picked token

bancorp

Interpolated KN

0.27 bits

NPLM

0.21 bits

Elman RNN

0.13 bits

RNN, dynamic evaluation

0.68 bits

Click a token to see, in panel B, how far back its loss reaches. Totals are closed-vocabulary bits; the second number is the spelling, which every model pays. On the whole held-out text, in bits per character: Interpolated KN 0.6371, NPLM 0.6001, Elman RNN 0.5545, RNN, dynamic evaluation 0.5158.

B · Gradients through time — products of Jacobians, against their bound

at the token's own state

0.178

at ⟨s⟩, the first state

0.176

Theorem 5's bound there

0.724

σ_max(W_h), ρ(W_h)

3.35, 2.13

Every per-step factor exceeds 1, so the bound only grows going back. The measured gradient does not follow it: the factors' expanding directions do not line up, and the product shrinks or holds steady while the bound explodes.

Medians over 12 initializations; gold ticks mark the LSTM's cell path σ(b_f)^60 (Proposition 4). Only the network whose gradient reaches the bit learns to report it.

The honest comparison

Each model is scored the perplexity topic’s way: a ratio estimate over held-out sentences, with a sentence-level standard error, per token and in bits per character of an open-vocabulary code. An unknown word is coded as <unk> followed by its letters under the perplexity topic’s speller. All models share one vocabulary and one speller, so the spelling bits are identical across rows ( $0.0219$ bits per character held out). The mixtures average a network’s prediction with Kneser–Ney’s, $\lambda\, p_{\text{neural}} + (1 - \lambda)\, p_{\text{KN}}$ , with $\lambda = 0.8$ for the NPLM and $0.9$ for the recurrent network, chosen on dev.

held out	nats per token (se)	perplexity	bits per character (se)	minus Kneser–Ney (paired se)
interpolated Kneser–Ney	$2.8551$ ( $0.0320$ )	$17.38$	$0.6371$ ( $0.0065$ )	—
NPLM	$2.6831$ ( $0.0301$ )	$14.63$	$0.6001$ ( $0.0059$ )	$-0.0371$ ( $0.0014$ )
NPLM with Kneser–Ney	$2.6688$ ( $0.0300$ )	$14.42$	$0.5970$ ( $0.0059$ )	$-0.0401$ ( $0.0013$ )
Elman RNN	$2.4717$ ( $0.0296$ )	$11.84$	$0.5545$ ( $0.0059$ )	$-0.0826$ ( $0.0017$ )
Elman RNN with Kneser–Ney	$2.4593$ ( $0.0294$ )	$11.70$	$0.5518$ ( $0.0059$ )	$-0.0853$ ( $0.0016$ )
Elman RNN, dynamic evaluation	$2.2919$ ( $0.0265$ )	$9.89$	$\mathbf{0.5158}$ ( $0.0054$ )	$-0.1213$ ( $0.0020$ )

Every network beats Kneser–Ney by more than twenty paired standard errors, and mixing in Kneser–Ney helps each a little. All models here share one tokenizer, vocabulary, and speller, so per-token perplexity and bits per character rank them identically; bits per character is the number that stays comparable when a later model tokenizes differently (the perplexity topic’s Theorem 5). The networks do it with more numbers: $677{,}987$ and $686{,}179$ parameters, against $42{,}458$ stored n-gram counts.

The last row is dynamic evaluation (Mikolov et al. 2010; Krause et al. 2018): after scoring each held-out sentence, take one gradient step on it, with a step size chosen on dev. It is an honest score because it is a code.

Proposition 5 (Dynamic evaluation is a code).

The arithmetic code of the perplexity topic’s Theorem 4 needs only a sequential model $Q(x_i \mid x_1 \dots x_{i-1})$ . If the model’s parameters are a deterministic function of the text already coded, it is such a model, and its codeword is within two bits of $-\log_2 Q(\text{text})$ .

Proof.

The decoder, having decoded $x_1 \dots x_{i-1}$ , can recompute the same parameters, so it can compute the same conditional and decode $x_i$ . The coding argument never used anything else about $Q$ .

∎

This is prequential coding (formalML, Minimum Description Length). Run exactly, the whole held-out text becomes one codeword of $79{,}614$ bits, within two bits of the model’s own $-\log_2 Q$ of $79{,}612.86$ . That codeword is encoded and its length checked; the decoder’s replay of every update is checked on two others, the EDGAR text’s codeword and a separate codeword of the first $150$ held-out sentences, and decoding each reproduces its text exactly. The perplexity topic’s Kneser–Ney code of the same text was $98{,}345$ bits, so the recurrent network with dynamic evaluation writes a file $19.0\%$ shorter.

Finance case study: what adapts, and what cannot

§ Finance case study: what adapts, and what cannot

The perplexity topic found that on public EDGAR MD&A excerpts, where $57.6\%$ of tokens are words the corpus never contained, Kneser–Ney predicts worse than a uniform guess, and lzma, given the training text, beats every model with $2.908$ bits per character. Here are the networks on the same text:

EDGAR	nats per token	perplexity	bits per character
uniform over $V$	$8.1394$	$3{,}427$	—
interpolated Kneser–Ney	$8.7963$	$6{,}610$	$5.6605$
NPLM	$6.4165$	$612$	$5.1176$
NPLM with Kneser–Ney	$6.1047$	$448$	$5.0465$
Elman RNN	$7.2049$	$1{,}346$	$5.2975$
Elman RNN with Kneser–Ney	$6.8996$	$992$	$5.2278$
Elman RNN, dynamic evaluation	$\mathbf{3.8986}$	$49.3$	$4.5432$
lzma, given the training text	—	—	$\mathbf{2.9080}$

The networks beat the uniform guess where Kneser–Ney cannot. On a text that is more than half unknown words, the <unk> events carry about half or more of each static model’s loss: $70.1\%$ of Kneser–Ney’s nats, $55.7\%$ of the recurrent network’s, and $49.8\%$ of the NPLM’s. Kneser–Ney’s is its uniform floor. The networks learned a rate, through word dropout and Theorem 2. Without dropout the NPLM is worse than uniform after one epoch of training ( $8.206$ nats per token), as Corollary 1 predicts, while after the same one epoch with dropout it scores $5.475$ .

In-domain training hurts out of domain, and the ranking reverses. By the selected epoch, both networks score worse on EDGAR than after their first epoch, while dev keeps improving: the NPLM goes from $5.475$ to $6.416$ nats per token, the recurrent network from $5.868$ to $7.205$ . Choosing the epoch on in-domain text picks a worse model for filings. The NPLM also beats the recurrent network on EDGAR, the reverse of the in-domain order: the longer context is longer memory of the training distribution.

Adaptation pays most. Dynamic evaluation cuts the recurrent network’s EDGAR loss from $7.20$ to $3.90$ nats per token, more than any change to the static model. It learns EDGAR’s own rate of unknown words, and its own contexts, as it codes. This is the perplexity topic’s “model that learns from the text it is coding”, now a neural one.

Proposition 6 (The spelling floor).

In any code that sends each unknown word as <unk> followed by its spelling under a speller $S$ , the bits per character are at least $\sum_{\text{unknown } w} -\log_2 S(w) \,/\, n_{\text{char}}$ .

The word-level part of the code costs at least zero bits, which is the whole proof. On EDGAR the floor is $3.6538$ bits per character, above lzma’s $2.9080$ . No model with this vocabulary and this speller can beat lzma there, however well it predicts words. The gap belongs to the closed vocabulary, which is a tokenizer’s problem (the subword topic), not a modeling one. Adaptivity helps here too: a cache that sends a repeated unknown word as a pointer to its earlier spelling, with the PPM-C escape probability of Moffat (1990), brings the spelling down to $2.7608$ bits per character, below lzma. The recurrent network with dynamic evaluation and the cache codes EDGAR at $3.6501$ bits per character, and lzma still wins. lzma adapts at every byte. Our models adapt above the word, and below it they still spell from a fixed table.

Where this goes next

Word embeddings. The NPLM’s embedding table is a by-product here. The word2vec topic makes it the product, and shows that its training objective factorizes a matrix of association scores, a low-rank matrix of the kind Theorem 6 is about.
Attention. The recurrent state compresses the whole prefix into $128$ numbers, and Theorem 5 says how fast it forgets. Attention replaces the compression with a weighted lookup over every earlier position, and the attention topic reads that lookup as kernel smoothing.
Scaling laws. Loss-versus-size curves are fit to exactly this cross-entropy. The unit matters there as it did here: bits per character survive a change of tokenizer, and nats per token do not.
Decoding. A tight model can be sampled or searched for its most probable sentence; the beam-search topic asks how well that search does.
Calibration. Theorem 2’s moment matching makes a trained model calibrated on average over its training data; the calibration topic asks whether it is calibrated event by event, and on text unlike its training data.

Honest caveats

⚠ Rigor flag

One corpus, one seed. The networks are small, trained once each from one random start, on a synthetic corpus. Their wins over Kneser–Ney are evidence about this generator, not about English, and variation across starts is not measured.
Tuning. Training is non-convex (Proposition 3), so nothing certifies the trained weights. The network sizes, learning rates, epoch counts, and the word-dropout treatment of <unk> were set in exploratory runs, some of which were scored on held-out text before the dev split existed. Everything the published harness tunes is tuned on dev.
Word dropout is a heuristic. Replacing once-seen words by <unk> half the time follows Good–Turing’s logic, but it is not an estimator with a guarantee, and the rate one half is a choice.
Three bounds are vacuous on the trained models. The Lipschitz constant of Theorem 1 is about twenty-five times the largest change observed; the tightness constant of Theorem 4 is $7.1 \times 10^{-17}$ ; the bound of Theorem 5 passes $10^{4}$ where the measured norm is $0.2$ . Each is true, and none describes what the networks actually do.
The LSTM result is partial. Proposition 4 is proved for the cell path only. The long-lag task is asserted for one start per model. In an exploratory rerun from eight further starts, not part of the harness, the bias-three network succeeded from seven, and the Elman network, whose gradient barely reaches the bit, from two.
The bottleneck is a rank statement. That language needs a high-rank model is a hypothesis, and on this corpus the energy beyond the bottleneck’s rank is small.
Exact codes need exact floats. The arithmetic codes are lossless because encoder and decoder compute identical floating-point numbers. That holds within one process and is not guaranteed across machines.
Small and narrow tests. EDGAR is $33$ sentences, and no character-level neural model was trained, so the comparison with lzma does not include the model best placed to win it.

Implementation

The companion notebook (notebookPath) imports neural_language_models.py, which owns every number on this page. It imports the perplexity topic’s speller and arithmetic coder, the n-gram topic’s smoothers and corpus, and the shared finance corpus, and trains both networks from scratch in NumPy on a laptop CPU. Running it as a script executes a harness of $23$ assertion-based tests in under a minute. The tests cover:

Proposition 1 on random and trained networks, Theorem 1 on random networks and on the trained NPLM, and the event-class table.
Proposition 2 and Theorem 2 on random softmax regressions, with complex-step Jacobians and L-BFGS stationary points; Corollary 1 by gradient descent and on the corpus, with and without word dropout; Proposition 3 on the trained network.
Theorem 3 on $200$ random computational graphs, and the hand-written gradients of both networks and of the long-lag models, against complex-step derivatives and central differences.
Theorem 4 on random networks with large $\varepsilon$ and on the trained one; Theorem 5 on random and trained networks, including the linear case and Gelfand’s formula; Proposition 4 exactly and to first order; the long-lag task from one start per model.
Theorem 6 on random instances, a mixture of softmaxes, and both trained networks; the comparison table; dynamic evaluation as an exact code, decoded back; and the finance case study.

The laboratory reads its data from a JSON file that one of those tests requires to equal the module’s viz_payload() exactly, so the figure cannot drift from the code.

References & Further Reading

paper A Neural Probabilistic Language Model — Yoshua Bengio, Réjean Ducharme, Pascal Vincent, and Christian Jauvin (2003) Journal of Machine Learning Research 3: 1137–55. Definition 1
paper Taylor Expansion of the Accumulated Rounding Error — Seppo Linnainmaa (1976) BIT Numerical Mathematics 16 (2): 146–60. Reverse-mode differentiation
paper The Complexity of Partial Derivatives — Walter Baur and Volker Strassen (1983) Theoretical Computer Science 22 (3): 317–30. The cost of a gradient is a constant times the cost of the function
paper Learning Representations by Back-Propagating Errors — David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams (1986) Nature 323 (6088): 533–36. Backpropagation in networks
paper Finding Structure in Time — Jeffrey L. Elman (1990) Cognitive Science 14 (2): 179–211. The simple recurrent network, Definition 3
paper Backpropagation Through Time: What It Does and How to Do It — Paul J. Werbos (1990) Proceedings of the IEEE 78 (10): 1550–60. Corollary 3
paper A Cache-Based Natural Language Model for Speech Recognition — Roland Kuhn and Renato De Mori (1990) IEEE Transactions on Pattern Analysis and Machine Intelligence 12 (6): 570–83. Cache models
paper Implementing the PPM Data Compression Scheme — Alistair Moffat (1990) IEEE Transactions on Communications 38 (11): 1917–21. The PPM-C escape used by the spelling cache
paper Learning Long-Term Dependencies with Gradient Descent Is Difficult — Yoshua Bengio, Patrice Simard, and Paolo Frasconi (1994) IEEE Transactions on Neural Networks 5 (2): 157–66. Vanishing gradients, Theorem 5
paper Long Short-Term Memory — Sepp Hochreiter and Jürgen Schmidhuber (1997) Neural Computation 9 (8): 1735–80. The LSTM and the constant error carousel
paper Using Complex Variables to Estimate Derivatives of Real Functions — William Squire and George Trapp (1998) SIAM Review 40 (1): 110–12. The complex-step derivative
paper Learning to Forget: Continual Prediction with LSTM — Felix A. Gers, Jürgen Schmidhuber, and Fred Cummins (2000) Neural Computation 12 (10): 2451–71. The forget gate
book Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation — Andreas Griewank and Andrea Walther (2008) 2nd ed. SIAM. Reverse accumulation and its cost
paper Recurrent Neural Network Based Language Model — Tomáš Mikolov, Martin Karafiát, Lukáš Burget, Jan Černocký, and Sanjeev Khudanpur (2010) Proceedings of Interspeech 2010, 1045–48. Recurrent language models and dynamic evaluation
book Matrix Analysis — Roger A. Horn and Charles R. Johnson (2012) 2nd ed. Cambridge University Press. Gelfand's formula (Corollary 5.6.14)
paper On the Difficulty of Training Recurrent Neural Networks — Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio (2013) Proceedings of ICML 2013, PMLR 28 (3): 1310–18. Spectral conditions and gradient clipping
paper Adam: A Method for Stochastic Optimization — Diederik P. Kingma and Jimmy Ba (2015) ICLR 2015. The optimizer used for every network here
paper An Empirical Exploration of Recurrent Network Architectures — Rafal Jozefowicz, Wojciech Zaremba, and Ilya Sutskever (2015) Proceedings of ICML 2015, PMLR 37: 2342–50. A forget-gate bias of one
book Deep Learning — Ian Goodfellow, Yoshua Bengio, and Aaron Courville (2016) MIT Press. Background on backpropagation and recurrent networks
paper Simple and Accurate Dependency Parsing Using Bidirectional LSTM Feature Representations — Eliyahu Kiperwasser and Yoav Goldberg (2016) Transactions of the ACL 4: 313–27. Word dropout to the unknown-word symbol
paper Recurrent Neural Networks as Weighted Language Recognizers — Yining Chen, Sorcha Gilroy, Andreas Maletti, Jonathan May, and Kevin Knight (2018) Proceedings of NAACL-HLT 2018, 2261–71. Consistency of recurrent models
paper Dynamic Evaluation of Neural Sequence Models — Ben Krause, Emmanuel Kahembwe, Iain Murray, and Steve Renals (2018) Proceedings of ICML 2018, PMLR 80: 2766–75. Proposition 5
paper Can Recurrent Neural Networks Warp Time? — Corentin Tallec and Yann Ollivier (2018) ICLR 2018. Chrono initialization of the forget bias
paper Breaking the Softmax Bottleneck: A High-Rank RNN Language Model — Zhilin Yang, Zihang Dai, Ruslan Salakhutdinov, and William W. Cohen (2018) ICLR 2018. Theorem 6 and the mixture of softmaxes
paper A Measure-Theoretic Characterization of Tight Language Models — Li Du, Lucas Torroba Hennigen, Tiago Pimentel, Clara Meister, Jason Eisner, and Ryan Cotterell (2023) Proceedings of ACL 2023, 9744–70. Tightness, Theorem 4