Solutions to Exercises

The density at the mean is $1 / (0.1\sqrt{2\pi}) \approx 3.99$. The probability of landing within 0.05 of 0 is $\Phi(0.5) - \Phi(-0.5) \approx 0.383$, where $\Phi$ is the standard Gaussian CDF. A density measures probability per unit length. Squeezing the distribution into a narrow range must raise its height so that the total area stays 1. A probability is an area, bounded by the total area of 1. For a narrow interval, $P(|X| < \delta) \approx 2\delta\, p(0)$, which is small however tall $p(0)$ is.

With a 20% base rate, flagged essays are $0.2 \times 0.95 = 0.19$ generated and $0.8 \times 0.05 = 0.04$ human, so

The detector’s error rates are properties of the detector, but the posterior also depends on the prior. When generated essays are rare, even a small false-positive rate applied to the large human majority produces more false alarms than true detections. The same arithmetic governs any classifier used to find a rare class.

Expand the variance of the sum into variances and covariances:

For i.i.d. draws every covariance is zero, which leaves $n\sigma^2 / n^2 = \sigma^2 / n$. With correlation $\rho$, each of the $n(n-1)$ covariance terms is $\rho\sigma^2$, so

However large $n$ becomes, the variance never falls below $\rho\sigma^2$. With $\rho = 0.1$ and $n = 32$, it is 0.128 instead of the i.i.d. 0.031: four times larger. A minibatch of near-duplicates behaves like a much smaller batch. The simulation below, checked in the tests, builds correlated draws from a shared component:

Up to a constant, the average negative log-likelihood is

Setting $\partial \ell / \partial \mu = -\frac{1}{\sigma^2 N}\sum_i (x_i - \mu) = 0$ gives $\hat\mu = \bar{x}$. Setting $\partial \ell / \partial \sigma = \frac{1}{\sigma} - \frac{1}{\sigma^3 N}\sum_i (x_i - \mu)^2 = 0$ gives $\hat\sigma^2 = \frac1N \sum_i (x_i - \bar{x})^2$. The tests confirm that the finite-difference gradient vanishes at these values.

For the bias, write $x_i - \bar{x} = (x_i - \mu) - (\bar{x} - \mu)$. Summing the squares, the cross terms combine into $\sum_i (x_i - \bar{x})^2 = \sum_i (x_i - \mu)^2 - N(\bar{x} - \mu)^2$. Take expectations, using $\Var[\bar{x}] = \sigma^2 / N$:

The sample mean sits closer to the data than the true mean does, so deviations from it are too small on average. Dividing by $N - 1$ instead removes the bias (NumPy’s ddof=1). With $N = 5$, the maximum-likelihood variance averages 0.8 of the truth:

With $\sigma$ fixed, the Gaussian negative log-likelihood is $\frac{1}{2\sigma^2}(y - \hat{y})^2$ plus terms that do not depend on the prediction. Minimizing its average is minimizing mean squared error. For Laplace noise,

so the loss is mean absolute error:

For a constant prediction, squared error is minimized by the mean of the targets and absolute error by the median. A few large outliers drag the mean but barely move the median, so MSE chases outliers and MAE resists them. Each loss is the right choice exactly when its noise model matches the data.

For a discrete distribution, the gradient passes through the finite sum. Then use $\nabla p = p\, \nabla \log p$:

With $f = 1$ the left side is the gradient of 1, so $\E[\nabla_\theta \log p_\theta] = 0$. Hence $\E[(f - b)\nabla_\theta \log p_\theta] = \E[f \nabla_\theta \log p_\theta]$ for any constant $b$.

For $x = 1 + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0, 1)$, the score is $\varepsilon$. Use the moments $\E[\varepsilon^2] = 1$, $\E[\varepsilon^4] = 3$, $\E[\varepsilon^6] = 15$, with odd moments zero:

The tests check all three means and standard deviations against these values.

Let $F(g) = e^{-e^{-g}}$ and $f(g) = e^{-g} e^{-e^{-g}}$ be the Gumbel CDF and density. Index $k$ wins when $z_j + g_j < z_k + g_k$ for every $j \ne k$. Condition on $g_k = g$; the other noises are independent:

The $j = k$ term in the sum comes from $f(g)$ itself. Write $S = \sum_j e^{z_j - z_k}$ and substitute $t = e^{-g}$, so that $\dd t = -e^{-g}\, \dd g$:

Scaling after the noise, $\argmax_k (z_k + g_k)/T$, picks the same index for any $T > 0$, so temperature must act on the logits first: $\argmax_k (z_k / T + g_k)$ samples $\softmax(\vz / T)$. The tests compare empirical frequencies from 200,000 draws with the softmax probabilities at $T = 0.5$ and $T = 2$.

The CDF is $F(x) = \int_0^x 2t\, \dd t = x^2$ on $[0, 1]$, so $F^{-1}(u) = \sqrt{u}$:

The mean is $\int_0^1 x \cdot 2x\, \dd x = 2/3$, and $P(X < 1/2) = F(1/2) = 1/4$. With 200,000 samples, both agree with these values to within 1%.

Surprisal is $-\log_2 p$ bits, or $-\ln p$ nats:

A language model with a cross-entropy of 2 nats per token is therefore doing far better than uniform guessing, which would cost 11.76 nats per token.

With $u(x) = 1/K$,

Gibbs' inequality makes the left side nonnegative, so $H(p) \le \log K$, with equality exactly when $p = u$. The tests check the identity for 200 random distributions.

Restrict the sum to outcomes with $p(x) > 0$, and let $Z = q(x)/p(x)$ with $x \sim p$. Jensen’s inequality for the concave logarithm gives

Equality in Jensen’s step requires $q/p$ to be constant where $p > 0$. Equality in the last step requires $q$ to put all its mass there. Together they force $q = p$.

Group the sum over samples by value. The value $x$ appears $N \hat{p}(x)$ times, so

and $H(\hat{p}, q) = H(\hat{p}) + \KL(\hat{p} \Vert q)$ splits it into a constant and a divergence. The tests confirm the equality numerically:

Over all distributions, Gibbs' inequality says the minimum is at $q = \hat{p}$: a model that only memorizes training frequencies, assigning zero probability to anything unseen. Real models are kept from that point by their parameterization, which cannot represent arbitrary tables and must share structure across inputs. They are also held back by regularization and early stopping, and by the softmax, which never outputs an exact zero. The gap between memorizing $\hat{p}$ and generalizing toward $p$ is the subject of the next part of the book.

The log-ratio of the densities is

Under $p$, $\E[(x - \mu_1)^2] = \sigma_1^2$ and $\E[(x - \mu_2)^2] = \sigma_1^2 + (\mu_1 - \mu_2)^2$. Taking expectations gives (7.6). With equal parameters the result is $0 + \tfrac12 - \tfrac12 = 0$. With equal variances it reduces to $(\mu_1 - \mu_2)^2 / (2\sigma^2)$: quadratic in the distance between the means. The tests compare the formula with numerical integration for three pairs of Gaussians.

Summing over the support of $q$, $\E_q[r] = \sum_x q(x) \frac{p(x)}{q(x)} = \sum_x p(x) = 1$. So $\E_q[r - 1] = 0$ and

For nonnegativity, the logarithm is concave, so it lies below its tangent line at 1: $\log r \le r - 1$ for every $r > 0$. Hence $k_3 = (r - 1) - \log r \ge 0$, with equality only at $r = 1$. Adding the zero-mean term $r - 1$ cancels much of $k_1$'s fluctuation, because $r - 1$ and $\log r$ move together.

$\KL(p \Vert q) = -H(p) - \E_p[\log q(x)]$, and only the second term depends on $q = \mathcal{N}(\mu, \sigma^2)$:

The mean enters only through $(\E_p[x] - \mu)^2$, which is minimized by $\mu = \E_p[x]$. Setting the derivative in $\sigma$ to zero gives $\sigma^2 = \Var_p[x]$. For the two-mode target, the mean is 0 and the variance is $0.6^2 + 2^2 = 4.36$, so $\sigma \approx 2.09$. The grid search in the tests finds the same values, to within its step size.

The reverse direction is $\E_q[\log q] - \E_q[\log p]$. The target’s log-density sits inside an expectation over the fitted distribution. For a mixture, that expectation has no closed form, and the objective has one local minimum per mode. Which one an optimizer finds depends on where it starts.

For $\left(\begin{smallmatrix} 0.30 & 0.10 \\ 0.15 & 0.45 \end{smallmatrix}\right)$, the marginals are $p(x) = (0.4, 0.6)$ and $p(y) = (0.45, 0.55)$:

When $Y = X$ always, for example with the diagonal table $\diag(0.2, 0.3, 0.5)$, knowing $Y$ removes all uncertainty about $X$. Then $I(X; Y) = H(X) \approx 1.030$ nats. For the product of the marginals, the joint equals the product, and the mutual information is exactly 0. All three cases are checked in the tests.

The accuracy is $248/400 = 0.62$. The plug-in standard error is $\sqrt{0.62(1-0.62)/400} = 0.024269$. A 95% normal interval is $0.62 \pm 1.96 \cdot 0.024269$, or 0.572 to 0.668. The standard error is the standard deviation we would expect for the accuracy estimate over repeated evaluation sets of the same size, not the model’s per-example error rate.

The median has no simple Bernoulli standard-error formula, but it is still a statistic of rows. The bootstrap approximates its sampling distribution by resampling rows with replacement and recomputing the median. For an LLM benchmark, the row should be the independent evaluation unit: a task, prompt, conversation, or judged comparison, not an individual token from inside one answer.

The five discordant rows contain four wins for B and one win for A, so the accuracy gap is $(4-1)/12 = 0.25$. Under the null, the five discordant rows are fair coin flips. Outcomes at least this extreme have zero or one B wins, or four or five B wins, so the p-value is $(1+5+5+1)/2^5 = 0.375$. The exact permutation test and McNemar’s test coincide here:

Rearrange $h = 1.96\sqrt{p(1-p)/N}$ to get $N = 1.96^2p(1-p)/h^2$, then round up. The worst case is $p=0.5$. With $h=0.02$, $N = 1.96^2 \cdot 0.25 / 0.02^2 = 2401$. The implementation in the chapter uses exactly that formula and ceil, because collecting a fraction of an example is impossible.

Population risk is $\E[\ell(f_\vtheta(\vx), y)]$, the average loss on future examples from the deployment distribution. Empirical risk replaces that expectation with a sample mean over the training set. The estimate can be biased if the training set was collected from a different population, filtered by a previous model, deduplicated incorrectly, or contaminated with test examples. More examples reduce sampling noise but do not fix a mismatched sampling process.

Start with $L(\vtheta)=\frac1N(\mX\vtheta-\vy)^\T(\mX\vtheta-\vy)$. Expanding gives $\frac1N(\vtheta^\T\mX^\T\mX\vtheta - 2\vy^\T\mX\vtheta + \vy^\T\vy)$. The gradient is $\frac2N\mX^\T\mX\vtheta - \frac2N\mX^\T\vy$. Setting it to zero and multiplying by $N/2$ gives $\mX^\T\mX\vtheta = \mX^\T\vy$. If $\mX^\T\mX$ is invertible, solve for $\vtheta$; otherwise use least squares or regularization.

For one example, $p=\sigma(z)$ and $z=\vx^\T\vw+b$. The BCE derivatives are $\partial \ell/\partial p = -y/p + (1-y)/(1-p)$ and $\partial p/\partial z = p(1-p)$. Multiplying gives $\partial \ell/\partial z = p-y$. The chain rule then gives $\partial \ell/\partial \vw = \vx(p-y)$. Averaging rows stacks those row vectors into $\mX^\T(\vp-\vy)/N$, whose shape matches $\vw$ because $\mX^\T$ maps per-example errors back to feature weights.

The tested snippet builds a noisy cubic training set, fits three polynomial feature maps, and measures validation loss against the clean signal. Degree 1 underfits because it cannot bend. Degree 11 overfits because it uses extra powers to chase noise. Degree 3 best matches the data generating function in this example.

The graph has leaves $x$ and $y$. It computes $xy$, $x^2$, $e^y$, then adds those three values and applies log. The node $x$ is shared because it feeds both $xy$ and $x^2$. During the backward pass, $x.grad$ must receive the contribution from the multiply path and the square path. Missing either contribution gives the derivative of a different program.

For $z=xy$, a small change gives $\dd z = y\,\dd x + x\,\dd y$. Multiplying by the output adjoint $\bar z$ gives $\bar x \mathrel{\char"2B}= \bar z y$ and $\bar y \mathrel{\char"2B}= \bar z x$. For $\mY=\mA\mW$, use the trace identity: $\langle \bar{\mY}, \dd\mA\mW + \mA\dd\mW\rangle$ = $\langle \bar{\mY}\mW^\T, \dd\mA\rangle + \langle \mA^\T\bar{\mY}, \dd\mW\rangle$. Thus $\bar{\mA}$ has the same shape as $\mA$, and $\bar{\mW}$ has the same shape as $\mW$.

Forward mode propagates one chosen input direction to all outputs, so computing a full gradient with respect to many parameters would require many sweeps. Reverse mode starts from the scalar loss adjoint and computes all parameter adjoints in one backward sweep, so it is the natural choice for training neural networks. Forward mode is attractive when the input dimension is small and the output dimension is large, or when only a few directional derivatives are needed.

The tested helper constructs the graph, calls backward, and returns gradients for all three inputs. The chapter test flattens $X$, $W$, and $b$, recomputes the scalar loss under small central-difference perturbations, and checks the concatenated autodiff gradient.

For sigmoid, $\sigma'(x)=\sigma(x)(1-\sigma(x))$. As $x\to\infty$, $\sigma(x)\to1$; as $x\to-\infty$, $\sigma(x)\to0$. Either way the product goes to 0. Since $\tanh'(x)=1-\tanh^2(x)$ and $\tanh(x)\to\pm1$, tanh also saturates. ReLU has slope 0 for negative inputs, so a unit that is always negative never updates from the loss signal. Leaky ReLU replaces that 0 slope by $\alpha$, leaving a path for gradients.

For $s=\sigma(x)=1/(1+e^{-x})$, $s'=e^{-x}/(1+e^{-x})^2=s(1-s)$, whose largest value is $1/4$ at $s=1/2$. Using $\tanh x=(e^x-e^{-x})/(e^x+e^{-x})$ gives $\tanh'(x)=1-\tanh^2(x)$, largest at 0 with value 1. SiLU is a product, so $(x\sigma(x))'=\sigma(x)+x\sigma(x)(1-\sigma(x))$. Softplus has derivative $e^x/(1+e^x)=\sigma(x)$. The test test_elementwise_activation_derivatives_match_finite_differences checks these formulas against central differences.

The exact GELU is $g(x)=x\Phi(x)$. Since $\Phi'(x)=\phi(x)$, the product rule gives $g'(x)=\Phi(x)+x\phi(x)$. The measured approximation error is computed, not guessed:

On the grid used by the tests, the maximum absolute error is $4.7324\times10^{-4}$, and the test asserts that exact measured value.

Ignoring biases, the plain MLP has $d(4d)+(4d)d=8d^2$ weights. The gated block has two input projections and one output projection, $dh+dh+hd=3dh$. Setting $3dh=8d^2$ gives $h=8d/3$.

The backward pass in scratch.activations.gated_ffn_backward applies the product rule to the gate and is checked with finite differences for $\vx$, $\mW$, $\mV$, and $\mW_2$.

For every class, $e^{z_i+c}/\sum_j e^{z_j+c}=e^c e^{z_i}/(e^c\sum_j e^{z_j})$, so the factor cancels. Temperature changes the gaps before this normalization:

The test checks that for $(2,1,-1)$ the largest probability is highest at $T=0.5$, lower at $T=1$, and closer to uniform at $T=2$.

Let $Z=\sum_k e^{z_k}$ and $p_i=e^{z_i}/Z$. For $i=j$, quotient rule gives $(e^{z_i}Z-e^{z_i}e^{z_i})/Z^2=p_i(1-p_i)$. For $i\ne j$, only the denominator changes, giving $-e^{z_i}e^{z_j}/Z^2=-p_i p_j$. Therefore $J=\diag(\vp)-\vp\vp^\T$. A row sum is $p_i-p_i\sum_jp_j=0$, matching shift invariance.

Using log-softmax, $L=-\sum_i y_i z_i+\sum_i y_i\log\sum_j e^{z_j}$. Since $\sum_i y_i=1$, this is $-\vy^\T\vz+\log\sum_j e^{z_j}$. Differentiating gives $-y_j+p_j$. With temperature, all logits inside softmax are $z_j/T$, so the chain rule multiplies the gradient by $1/T$.

The implementation computes log-softmax by shifting logits, then returns the mean loss and the gradient with respect to logits:

The tests check three cases: the plain one-hot gradient $(\vp-\vy)/B$, a temperature and label-smoothing gradient, and the z-loss gradient added to cross-entropy.

For Gaussian noise with fixed $\sigma$, $-\log p(y\mid\hat y)=(y-\hat y)^2/(2\sigma^2)\log\sigma\tfrac12\log2\pi$. Only the squared residual depends on $\hat y$, so maximum likelihood minimizes MSE. For Laplace noise, $-\log p(y\mid\hat y)=|y-\hat y|/b+\log(2b)$, so the prediction-dependent part is MAE. The scale changes the gradient size but not the optimum when fixed.

With $r=\hat y-y$, MSE has gradient $2r$, MAE has $\sign(r)$ away from zero, and Huber has $r$ for $|r|\le\delta$ and $\delta\sign(r)$ outside. The cap is visible in the tested helper:

For BCE, substitute $\sigma(z)=1/(1+e^{-z})$ and simplify separately for $z\ge0$ and $z<0$ to get $\max(z,0)-zy+\log(1+e^{-|z|})$. Differentiating gives $\sigma(z)-y$. The stable implementation stays finite even for logits $\pm1000$:

Focal loss multiplies BCE by $(1-p_t)^\gamma$. When $p_t\approx1$, the multiplier is near zero, so easy examples contribute little. When $p_t$ is small, the multiplier is near one, so hard examples keep a cross-entropy-like signal. Forward KL, $\KL(p\Vert q)$, averages over the target and punishes missing target mass. Reverse KL, $\KL(q\Vert p)$, averages over the model and is more willing to focus on one mode. Jensen-Shannon is symmetric because it averages $\KL(p\Vert m)$ and $\KL(q\Vert m)$ with the same midpoint $m=(p+q)/2$.

For softened student probabilities $q_T=\softmax(z/T)$ and fixed teacher $p_T$, the cross-entropy gradient is $(q_T-p_T)/T$. Multiplying the loss by $T^2$ makes it $T(q_T-p_T)$. The tests use this wrapper:

They check that the scaled gradient is $T^2$ times the unscaled one for the same softened distributions, and that the analytic gradient matches finite differences.

The arrays have shapes $\mZ_1, \mH \in \R^{B \times 10}$ and $\mZ_2 \in \R^{B \times 4}$. The parameters are $\mW_1 \in \R^{5 \times 10}$, $\vb_1 \in \R^{10}$, $\mW_2 \in \R^{10 \times 4}$, and $\vb_2 \in \R^4$. The count is $5\cdot10 + 10 + 10\cdot4 + 4 = 104$, checked by the helper:

The softmax derivative is $\partial p_k/\partial z_j = p_k(\delta_{kj}-p_j)$. Therefore

For the mean loss, $L = B^{-1}\sum_i \ell_i$, so every row gradient is divided by $B$. All later matrix products are linear in that upstream gradient, so dividing again would make the gradient too small by another factor of $B$. The tests check this gradient against finite differences.

Because the terms are independent and zero mean, $\Var[w_i x_i] = \E[w_i^2]\E[x_i^2] = \Var[w]\Var[x]$. Variances of independent sums add, so $\Var[z] = n\Var[w]\Var[x]$. Choosing $\Var[w] = 1/n$ keeps the preactivation scale comparable to the input scale. A ReLU makes about half of a symmetric preactivation zero, so the second moment is roughly halved; He initialization compensates with $\Var[w] = 2/n$.

The tested implementation uses float64 copies for gradient checking and float32 for the training run. Its SGD loop lowers the synthetic-data loss and reaches high training accuracy. On a tiny set with permuted labels, it can fit the update examples better than it agrees with the original clean labels:

That gap is overfitting: optimization succeeded on the training objective, but the fitted rule matched noise instead of the data-generating pattern.

For eigenvalue $\lambda_i$, the error multiplier is $1-\eta\lambda_i$. Stability requires $|1-\eta\lambda_i| < 1$ for every eigenvalue, so $0 < \eta < 2/25$. With $\eta=1/25$, the $25$ direction goes to zero in one step, but the $1$ direction is multiplied by $24/25$ each step. The best fixed-rate worst-case factor is:

Unrolling the recurrence gives

The same geometric sum gives $v_t=(1-\beta_2^t)g^2$. Without correction, both moment estimates are biased toward zero at small $t$. Dividing by $1-\beta_1^t$ and $1-\beta_2^t$ makes the constant-gradient estimates equal to $g$ and $g^2$.

With L2 regularization, Adam sees $g+\lambda\theta$ and then divides by the adaptive $\sqrt{\hat{v}}+\epsilon$, so the decay part is coordinate-scaled like any other gradient. AdamW instead applies $\theta \leftarrow (1-\eta\lambda)\theta$ separately, then takes the adaptive gradient step. For clipping, $13 > 5$, so every tensor is scaled by $5/13$:

For a matrix $\mM = \mU\mS\mV^\T$, the nearest orthogonal polar factor in Frobenius norm is $\mU\mV^\T$. The Newton-Schulz iteration approximates that factor using only matrix multiplies, which is cheaper than an SVD in a training step. The tests compare the chapter implementation with NumPy’s SVD on tall and wide matrices. Muon is restricted here to 2-D hidden weights because biases, gains, embeddings, and output heads do not have the same interior matrix geometry.

BatchNorm estimates $\mu_j$ and $\sigma_j^2$ from the current minibatch during training, but inference may use a different batch size, or one example at a time, so it uses running statistics accumulated during training. For $B \times d$, BatchNorm reduces over the $B$ axis for each feature. LayerNorm reduces over the $d$ axis inside each row, so it uses the current example’s own statistics in both training and inference.

Let $\bar{\hat{\vx}}=\bar{\vy}\odot\boldsymbol{\gamma}$ and $\hat{\vx}=(\vx-\mu)/s$. A change in one input coordinate affects the normalized output directly, through the row mean, and through the row variance. Collecting those terms gives

The first mean removes the component caused by mean subtraction. The second removes the component caused by changing the row’s variance. The tests check this formula with finite differences.

For $a>0$ and $\epsilon=0$,

The helper checks the same invariance numerically:

For a residual block $\vx_{l+1}=\vx_l+F_l(\vx_l)$, differentiating gives an identity term $\bar{\vx}_{l+1}$ plus the gradient through $F_l$. Even if the learned branch is small or badly scaled, the identity term passes gradient backward.

Inverted dropout samples a Bernoulli keep mask and divides kept activations by the keep probability, so its expectation equals the input. The tested helper estimates that mean:

The precision demo casts large float32 values to fp16 and rounded bfloat16, showing fp16 overflow while bfloat16 remains finite. It also multiplies a tiny gradient by a scale before the fp16 cast and divides after, recovering a nonzero unscaled gradient. Float32 master weights then accumulate updates that would be lost if every step were rounded to a low-precision copy.

Characters make unknown text easy and keep the alphabet small, but they make the sequence long. Words give short, readable sequences for common text, but unhappiness and the emoji may need unknown-token handling unless both were in the vocabulary. Byte-level subwords keep exact coverage because every string is bytes first, and frequent chunks such as un or ness can be merged; their drawback is that rare text may still split into many small pieces.

Stack the one-hot rows into $\mX$, so the embedding output is $\mY = \mX\mE$. For a scalar loss,

Column $j$ of $\mX^\T$ is 1 exactly for positions whose id is $j$, so row $j$ of $\bar{\mE}$ is the sum of those upstream rows. The checked helper shows the repeated id receiving both contributions:

The trainer in Section 17.2 counts pairs after each rewrite, so later merges can combine pieces that did not exist at the start. A compact way to check the result is to compare corpus length before and after BPE: