CS229: Machine Learning — Modern Reader

1. \(f = \tfrac12\sum_{i,j}A_{ij}x_ix_j + \sum_i b_ix_i\). Then \(\partial f/\partial x_k = \tfrac12\big(\sum_j A_{kj}x_j + \sum_i A_{ik}x_i\big) + b_k = (Ax)_k + b_k\), using \(A_{ik} = A_{ki}\). Differentiating \((Ax)_k\) by \(x_l\) gives \(A_{kl}\), so the Hessian is \(A\).
2. Chain rule: \(\partial f/\partial x_k = g'(a^\top x)\,a_k\), so \(\nabla f = g'(a^\top x)\,a\). Differentiating again by \(x_l\): \(g''(a^\top x)\,a_ka_l\), which is entry \((k,l)\) of \(g''(a^\top x)\,aa^\top\).
3. Take \(g(t) = \log(1+e^{-t})\) and \(a = yx\). Then \(g'(t) = -1/(1+e^{t})\), so \(\nabla_\theta f = -\dfrac{y\,x}{1+\exp(y\,\theta^\top x)}\). Also \(g''(t) = e^{t}/(1+e^{t})^2 \gt 0\) and \(aa^\top = y^2xx^\top = xx^\top\), so the Hessian is a positive number times \(xx^\top\), which is PSD by P2(a). The reference code checks the gradient numerically.

1. \(x^\top zz^\top x = (z^\top x)^2\ge0\). Every column of \(zz^\top\) is a multiple of \(z\), so the rank is 1, and \(zz^\top x = z(z^\top x) = 0\) exactly when \(x\perp z\): the null space is the set of vectors orthogonal to \(z\).
2. \(v^\top X^\top Xv = (Xv)^\top(Xv) = \|Xv\|^2\ge0\). It is zero only when \(Xv = 0\), and \(Xv=0\) has only the solution \(v=0\) exactly when the columns are independent.
3. \(x^\top(A+B)x = x^\top Ax + x^\top Bx\ge0\) and \(x^\top(cA)x = c\,x^\top Ax\ge0\). Each \(x^{(i)}x^{(i)\top}\) is PSD by (a), so the sum is too. That sum equals \(X^\top X\) for the design matrix \(X\) whose rows are the \(x^{(i)\top}\), which is the matrix in the normal equations.
4. \(4x_1^2 - 4x_1x_2 + x_2^2 = (2x_1 - x_2)^2\ge0\), so \(A\) is PSD. It is not positive definite: the form is zero along the line \(x_2 = 2x_1\). Its eigenvalues are 5 and 0.

1. \(\operatorname{tr}(Q\Lambda Q^\top) = \operatorname{tr}(\Lambda Q^\top Q) = \operatorname{tr}(\Lambda) = \sum_i\lambda_i\).
2. With \(y = Q^\top x\), \(x^\top Ax = y^\top\Lambda y = \sum_i\lambda_iy_i^2\). Since \(Q\) is invertible, \(y\) ranges over all of \(\mathbb{R}^n\) as \(x\) does. If every \(\lambda_i\ge0\) the sum is nonnegative. Conversely, taking \(x = q_i\) gives \(q_i^\top Aq_i = \lambda_i\), so PSD forces \(\lambda_i\ge0\).
3. \(Q^\top\) preserves length, so \(\|y\| = \|x\| = 1\). Then \(\sum_i\lambda_iy_i^2\le\lambda_{\max}\sum_iy_i^2 = \lambda_{\max}\), with equality when \(y\) puts all its weight on the largest eigenvalue, that is, \(x = q_{\max}\).
4. The eigenvalues are 3 and 1, with eigenvectors \((1,1)/\sqrt2\) and \((1,-1)/\sqrt2\). The trace is \(2+2 = 4 = 3+1\). The maximum of \(x^\top Ax\) over unit vectors is 3, at \((1,1)/\sqrt2\); the reference code confirms it by random search.

1. Bayes' rule: \(\dfrac{0.95\cdot0.01}{0.95\cdot0.01 + 0.05\cdot0.99} = \dfrac{0.0095}{0.059}\approx0.161\). Only about 16%, because healthy people vastly outnumber sick ones, so false positives dominate.
2. Write the count as \(\sum_i I_i\) with \(I_i = 1\) if flip \(i\) is heads. \(\mathbb{E}[I_i] = p\), so the mean is \(np\). Each \(\operatorname{Var}(I_i) = p - p^2 = p(1-p)\), and variances of independent variables add, so the variance is \(np(1-p)\).
3. \(a^\top X - a^\top\mu = a^\top(X-\mu)\), so \(\operatorname{Var}(a^\top X) = \mathbb{E}\big[a^\top(X-\mu)(X-\mu)^\top a\big] = a^\top\Sigma a\). A variance cannot be negative, so \(a^\top\Sigma a\ge0\) for every \(a\): \(\Sigma\) is PSD.
4. Expectation is linear, so \(\mathbb{E}[Y] = A\mu + b\). Then \(Y - \mathbb{E}[Y] = A(X-\mu)\), and \(\operatorname{Cov}(Y) = \mathbb{E}\big[A(X-\mu)(X-\mu)^\top A^\top\big] = A\Sigma A^\top\).

1. \(\ell(\mu,\sigma^2) = -\tfrac n2\log(2\pi\sigma^2) - \tfrac{1}{2\sigma^2}\sum_i(x^{(i)}-\mu)^2\). Setting \(\partial\ell/\partial\mu = \tfrac1{\sigma^2}\sum_i(x^{(i)}-\mu) = 0\) gives \(\hat\mu = \tfrac1n\sum_ix^{(i)}\). Setting \(\partial\ell/\partial\sigma^2 = -\tfrac{n}{2\sigma^2} + \tfrac1{2\sigma^4}\sum_i(x^{(i)}-\mu)^2 = 0\) gives \(\hat\sigma^2 = \tfrac1n\sum_i(x^{(i)}-\hat\mu)^2\).
2. \(\Sigma^{-1} = Q\Lambda^{-1}Q^\top\). With \(y = Q^\top(x-\mu)\), the contour becomes \(\sum_iy_i^2/\lambda_i = c\): an axis-aligned ellipse in \(y\)-coordinates with half-lengths \(\sqrt{c\lambda_i}\). Since \(y\) holds the coordinates of \(x-\mu\) along \(q_1,\dots,q_n\), the axes point along the eigenvectors. A large eigenvalue means a long axis: the data spreads most in that direction, which is what PCA finds.
3. By P4(d) with \(A = Q\Lambda^{1/2}\) and \(\operatorname{Cov}(z) = I\): \(\operatorname{Cov}(x) = Q\Lambda^{1/2}\Lambda^{1/2}Q^\top = Q\Lambda Q^\top = \Sigma\). The reference code below matches \(\mu\) and \(\Sigma\) to about three decimal places.

import numpy as np

def logistic_loss_grad(theta, x, y):
    """f(theta) = log(1 + exp(-y theta^T x)), y in {-1, +1}."""
    m = y * theta @ x
    return np.log1p(np.exp(-m)), -y * x / (1 + np.exp(m))

def sample_gaussian(mu, Sigma, n, rng):
    """x = mu + Q Lambda^{1/2} z with Sigma = Q Lambda Q^T and z ~ N(0, I)."""
    lam, Q = np.linalg.eigh(Sigma)
    z = rng.normal(size=(n, len(mu)))
    return mu + z @ (Q * np.sqrt(lam)).T

if __name__ == "__main__":
    rng = np.random.default_rng(0)

    # P1(c): logistic loss gradient vs finite differences
    theta, x, y = rng.normal(size=4), rng.normal(size=4), -1.0
    f, g = logistic_loss_grad(theta, x, y); eps = 1e-6
    num = np.array([(logistic_loss_grad(theta + eps * e, x, y)[0] - logistic_loss_grad(theta - eps * e, x, y)[0]) / (2 * eps) for e in np.eye(4)])
    print(f"P1 logistic gradient max error: {np.abs(g - num).max():.1e}")

    # P2(d) and P3(d): eigenvalues of the two example matrices
    print("P2 eigenvalues of [[4,-2],[-2,1]]:", np.round(np.linalg.eigvalsh([[4.0, -2], [-2, 1]]), 6))
    A = np.array([[2.0, 1], [1, 2]]); lam, Q = np.linalg.eigh(A)
    print("P3 eigenvalues of [[2,1],[1,2]]:", lam, "| trace:", np.trace(A))
    xs = rng.normal(size=(100000, 2)); xs /= np.linalg.norm(xs, axis=1, keepdims=True)
    print(f"P3 max of x^T A x over 100k random unit vectors: {np.einsum('ij,jk,ik->i', xs, A, xs).max():.4f}")

    # P4(a): Bayes' rule
    prior, sens, fpr = 0.01, 0.95, 0.05
    print(f"P4 P(disease | positive) = {sens*prior / (sens*prior + fpr*(1-prior)):.3f}")

    # P5(c): sampling with the spectral decomposition
    mu = np.array([1.0, -2.0]); Sigma = np.array([[2.0, 1.2], [1.2, 1.0]])
    X = sample_gaussian(mu, Sigma, 200000, rng)
    print("P5 sample mean:", np.round(X.mean(0), 3), " sample covariance:", np.round(np.cov(X.T), 3).tolist())

P1 logistic gradient max error: 8.4e-11
P2 eigenvalues of [[4,-2],[-2,1]]: [0. 5.]
P3 eigenvalues of [[2,1],[1,2]]: [1. 3.] | trace: 4.0
P3 max of x^T A x over 100k random unit vectors: 3.0000
P4 P(disease | positive) = 0.161
P5 sample mean: [ 0.999 -2.002]  sample covariance: [[2.006, 1.202], [1.202, 1.001]]

1. Expand \(J = \tfrac12(\theta^\top X^\top X\theta - 2\vec y^\top X\theta + \vec y^\top\vec y)\) and use \(\nabla\theta^\top A\theta = 2A\theta\) for symmetric \(A\). The Hessian \(X^\top X\) is PSD since \(v^\top X^\top Xv = \|Xv\|^2\), and positive definite when invertible, so the stationary point \(X^\top X\theta = X^\top\vec y\) is the unique minimizer.
2. Since \(X^\top X\theta^\star = X^\top\vec y\), the gradient is \(X^\top X(\theta - \theta^\star)\). Subtract \(\theta^\star\) from both sides of the update.
3. Diagonalize \(X^\top X = Q\Lambda Q^\top\). In the eigenbasis each coordinate is multiplied by \(1 - \alpha\lambda_i\) per step, so every coordinate shrinks to zero iff \(|1 - \alpha\lambda_i| \lt 1\) for all \(i\), that is, \(0 \lt \alpha \lt 2/\lambda_i\) for all \(i\). The binding constraint is \(\lambda_{\max}\). Slow convergence along small-\(\lambda\) directions is why badly scaled features make gradient descent slow.
4. See the reference code (first half). Just under the limit converges; just over it diverges.

1. \(\sigma(z) = (1 + e^{-z})^{-1}\), so \(\sigma'(z) = e^{-z}(1 + e^{-z})^{-2} = \sigma(z)\cdot\frac{e^{-z}}{1+e^{-z}} = \sigma(z)(1 - \sigma(z))\).
2. With \(z_2 = W_2h + b_2\), cross-entropy gives \(\delta_2 = p - e_y\) (A2(c)). Then \(\partial\mathcal L/\partial W_2 = \delta_2h^\top\) and \(\partial\mathcal L/\partial b_2 = \delta_2\). Backpropagating through \(W_2\) and the sigmoid, \(\delta_1 = (W_2^\top\delta_2)\odot h\odot(1-h)\), so \(\partial\mathcal L/\partial W_1 = \delta_1x^\top\) and \(\partial\mathcal L/\partial b_1 = \delta_1\).
3. See the reference code (second half); the maximum error should be around \(10^{-10}\).

import numpy as np

# ---- S1: normal equations and the gradient descent step-size limit ----
def gd(X, y, alpha, steps):
    theta = np.zeros(X.shape[1])
    for _ in range(steps):
        theta -= alpha * X.T @ (X @ theta - y)       # gradient of J = 0.5 ||X theta - y||^2
    return theta

# ---- S2: backprop for a two-layer network (sigmoid hidden, softmax output) ----
def sigmoid(z): return 1 / (1 + np.exp(-z))

def softmax(z):
    z = z - z.max(); e = np.exp(z); return e / e.sum()

def forward_backward(x, y, W1, b1, W2, b2):
    h = sigmoid(W1 @ x + b1)
    p = softmax(W2 @ h + b2)
    loss = -np.log(p[y])
    d2 = p.copy(); d2[y] -= 1                         # dL/dz2 = p - e_y
    dW2, db2 = np.outer(d2, h), d2
    d1 = (W2.T @ d2) * h * (1 - h)                    # dL/dz1
    dW1, db1 = np.outer(d1, x), d1
    return loss, (dW1, db1, dW2, db2)

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    X = rng.normal(size=(100, 3)) @ np.diag([3.0, 1.0, 0.3]); y = X @ np.array([1.0, -2.0, 0.5]) + 0.1 * rng.normal(size=100)
    theta_ne = np.linalg.solve(X.T @ X, X.T @ y)
    lam_max = np.linalg.eigvalsh(X.T @ X).max()
    print(f"2 / lambda_max = {2/lam_max:.5f}")
    for frac in [0.5, 0.99, 1.01]:
        th = gd(X, y, frac * 2 / lam_max, 5000)
        print(f"alpha = {frac:4.2f} x limit: |theta - theta_normal_eq| = {np.linalg.norm(th - theta_ne):.3e}")

    d, hdim, k = 4, 5, 3
    W1, b1 = rng.normal(size=(hdim, d)), rng.normal(size=hdim)
    W2, b2 = rng.normal(size=(k, hdim)), rng.normal(size=k)
    x, yy = rng.normal(size=d), 2
    _, grads = forward_backward(x, yy, W1, b1, W2, b2)
    params = [W1, b1, W2, b2]; eps = 1e-6; err = 0.0
    for P, G in zip(params, grads):
        for idx in np.ndindex(P.shape):
            old = P[idx]
            P[idx] = old + eps; lp = forward_backward(x, yy, *params)[0]
            P[idx] = old - eps; lm = forward_backward(x, yy, *params)[0]
            P[idx] = old
            err = max(err, abs((lp - lm) / (2 * eps) - G[idx]))
    print(f"backprop vs finite differences, max error: {err:.2e}")

2 / lambda_max = 0.00186
alpha = 0.50 x limit: |theta - theta_normal_eq| = 4.213e-15
alpha = 0.99 x limit: |theta - theta_normal_eq| = 1.897e-15
alpha = 1.01 x limit: |theta - theta_normal_eq| = 9.851e+42
backprop vs finite differences, max error: 1.85e-10

1. \(\mathbb{E}[q^\top k] = \sum_i \mathbb{E}[q_i]\mathbb{E}[k_i] = 0\). The terms \(q_i k_i\) are independent, so \(\operatorname{Var}(q^\top k) = \sum_i \mathbb{E}[q_i^2]\,\mathbb{E}[k_i^2] = d\).
2. Scores have standard deviation \(\sqrt d\), so the gap between the largest score and the rest grows with \(d\) and the softmax approaches a one-hot vector. At a one-hot \(s\), every entry of \(J = \operatorname{diag}(s) - ss^\top\) is zero (each diagonal entry is \(s_i(1-s_i) = 0\)), so almost no gradient flows back into the queries and keys. Dividing by \(\sqrt d\) restores unit variance regardless of head size.

1. With \(Z = \sum_k e^{z_k}\) and \(s_i = e^{z_i}/Z\): \(\partial s_i/\partial z_j = \mathbf{1}\{i=j\}e^{z_i}/Z - e^{z_i}e^{z_j}/Z^2 = s_i\mathbf{1}\{i=j\} - s_i s_j\).
2. \(J\mathbf{1} = s - s(s^\top\mathbf{1}) = s - s = 0\), so \(\mathbf 1\) is in the null space and the rank is at most \(n-1\). Symmetry is clear. For any \(v\), \(v^\top J v = \sum_i s_i v_i^2 - (\sum_i s_i v_i)^2\), which is the variance of \(v_I\) when the index \(I\) is drawn from \(s\), so it is nonnegative.
3. \(-\log s_y = -z_y + \log\sum_j e^{z_j}\). The gradient is \(-e_y + s\), the same \(\phi - e_y\) the notes derive for the cross-entropy loss.

1. In layer 1, position \(t\) outputs \(\sum_{j\le t} a_{tj} v_j\), where \(a_{tj}\) depends on \(q_t\) and \(k_j\) for \(j \le t\). Every input is a function of \(x_{\le t}\). The MLP, residual connections and layer norm act on each position separately. If every layer-\(\ell\) output at positions \(\le t\) depends only on \(x_{\le t}\), the same argument applies to layer \(\ell+1\).
2. Position \(t\)'s output is a valid function of the prefix only, so its prediction for \(x_{t+1}\) is a legitimate conditional model, and all positions are computed in parallel from the true sequence (teacher forcing). At generation time \(x_{t+1}\) does not exist until it has been sampled, so it cannot be fed in earlier.
3. Keys and values of earlier positions never change once computed, because later tokens cannot affect them. Each layer stores one key and one value vector of total width \(d\) (summed over heads) per token, so the cache is \(2LTdb\) bytes. For the given numbers: \(2\cdot32\cdot8192\cdot4096\cdot2 = 4{,}294{,}967{,}296\) bytes, about 4.3 GB for one sequence. This is why attention variants that share keys and values across heads (Ch. 17.4) matter in practice.

1. Each head has three \(d \times (d/h)\) projections, so the \(h\) heads together have \(3d^2\), plus \(d^2\) for \(W_O\): \(4d^2\). The MLP has \(d\cdot4d + 4d\cdot d = 8d^2\). Total \(12d^2\). Changing \(h\) only reshapes the same matrices.
2. \(12Ld^2 + Vd\). For GPT-2 small: \(12\cdot12\cdot768^2 = 84{,}934{,}656\), \(Vd = 38{,}597{,}376\), positions \(1024\cdot768 = 786{,}432\). Sum \(\approx 124.3\)M, matching the reported 124M; biases and layer norms add a little more.

\(\tilde q_m^\top \tilde k_n = q^\top R(m\omega)^\top R(n\omega) k\). Rotations satisfy \(R(\alpha)^\top = R(-\alpha)\) and \(R(\alpha)R(\beta) = R(\alpha+\beta)\), so the product is \(R((n-m)\omega)\). In \(d\) dimensions RoPE applies this to each pair of coordinates with its own frequency \(\omega_i\), and the same argument holds pair by pair.

import numpy as np

def softmax(z, axis=-1):
    z = z - z.max(axis=axis, keepdims=True)
    e = np.exp(z)
    return e / e.sum(axis=axis, keepdims=True)

def causal_self_attention(X, Wq, Wk, Wv):
    """X: (T, d_model). Returns (T, d_v)."""
    T = X.shape[0]
    Q, K, V = X @ Wq, X @ Wk, X @ Wv
    scores = Q @ K.T / np.sqrt(Q.shape[1])
    mask = np.triu(np.ones((T, T), dtype=bool), k=1)   # True above diagonal = future
    scores = np.where(mask, -np.inf, scores)
    return softmax(scores, axis=1) @ V

def multi_head(X, heads, Wo):
    """heads: list of (Wq, Wk, Wv), each projecting to d_model // h."""
    return np.concatenate([causal_self_attention(X, *w) for w in heads], axis=1) @ Wo

def softmax_jacobian(z):
    s = softmax(z)
    return np.diag(s) - np.outer(s, s)

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    T, d, h = 6, 16, 4
    dh = d // h
    heads = [tuple(rng.normal(size=(d, dh)) / np.sqrt(d) for _ in range(3)) for _ in range(h)]
    Wo = rng.normal(size=(d, d)) / np.sqrt(d)
    X = rng.normal(size=(T, d))
    Y = multi_head(X, heads, Wo)

    # (c) causality: changing token t+1.. must not change outputs at positions <= t
    t = 2
    X2 = X.copy(); X2[t + 1:] = rng.normal(size=(T - t - 1, d))
    Y2 = multi_head(X2, heads, Wo)
    print("causal ok:", np.allclose(Y[:t + 1], Y2[:t + 1]), "| future changed:", not np.allclose(Y[t + 1:], Y2[t + 1:]))

    # (d) softmax Jacobian vs finite differences
    z = rng.normal(size=5); J = softmax_jacobian(z); eps = 1e-6
    Jnum = np.stack([(softmax(z + eps * e) - softmax(z - eps * e)) / (2 * eps) for e in np.eye(5)], axis=1)
    print("jacobian max err:", np.abs(J - Jnum).max(), "| rank:", np.linalg.matrix_rank(J), "| row sums ~0:", np.allclose(J.sum(1), 0))

    # (e) why 1/sqrt(d): variance of q.k grows like d
    for dd in [16, 256, 4096]:
        q, k = rng.normal(size=(20000, dd)), rng.normal(size=(20000, dd))
        print(f"d={dd:5d}  Var(q.k)={np.var((q * k).sum(1)):9.1f}  Var(q.k/sqrt d)={np.var((q * k).sum(1) / np.sqrt(dd)):.3f}")

causal ok: True | future changed: True
jacobian max err: 4.850038426429393e-11 | rank: 4 | row sums ~0: True
d=   16  Var(q.k)=     16.1  Var(q.k/sqrt d)=1.009
d=  256  Var(q.k)=    255.6  Var(q.k/sqrt d)=0.998
d= 4096  Var(q.k)=   4183.0  Var(q.k/sqrt d)=1.021

\(\log p(x)\) does not depend on \(z\), so \(\log p(x) = \mathbb{E}_q[\log p(x,z) - \log p(z\mid x)]\). Add and subtract \(\log q(z)\) inside the expectation: the first piece is the ELBO and the second is \(\mathbb{E}_q[\log q(z) - \log p(z\mid x)] = \mathrm{KL}(q\,\|\,p(\cdot\mid x)) \ge 0\), which is zero iff \(q = p(\cdot \mid x)\).

1. \(\mathbb{E}_q\big[-\tfrac12\log(2\pi\sigma^2) - \tfrac{(z-\mu)^2}{2\sigma^2} + \tfrac12\log 2\pi + \tfrac{z^2}{2}\big] = -\tfrac12\log\sigma^2 - \tfrac12 + \tfrac12(\mu^2 + \sigma^2)\), using \(\mathbb{E}_q[z^2] = \mu^2+\sigma^2\).
2. The distribution of \(\varepsilon\) does not depend on \((\mu,\sigma)\), so the derivative moves inside the expectation and the chain rule gives \(f'(z)\cdot\partial z/\partial\mu = f'(z)\) and \(f'(z)\cdot\varepsilon\). Without the rewrite the parameters sit inside the sampling distribution; the general alternative is the score-function estimator \(\mathbb{E}[f(z)\nabla\log q(z)]\), which is unbiased but usually far noisier. Set D uses that estimator for RL, where the sampling step cannot be reparameterized.

Induct. If \(x_{t-1} = \sqrt{\bar\alpha_{t-1}}x_0 + \sqrt{1-\bar\alpha_{t-1}}\,\varepsilon'\), then \(x_t = \sqrt{\bar\alpha_t}x_0 + \sqrt{\alpha_t(1-\bar\alpha_{t-1})}\,\varepsilon' + \sqrt{\beta_t}\,\varepsilon_t\). A sum of independent Gaussians is Gaussian with variance \(\alpha_t - \bar\alpha_t + 1 - \alpha_t = 1-\bar\alpha_t\). When \(\bar\alpha_T\approx0\), \(x_T\) is close to \(\mathcal N(0,I)\) whatever \(x_0\) was, which is why sampling starts from pure noise.

1. Given \(x_0\): \(x_{t-1}\) has mean \(\sqrt{\bar\alpha_{t-1}}x_0\) and variance \(1-\bar\alpha_{t-1}\); \(x_t\) has mean \(\sqrt{\bar\alpha_t}x_0\) and variance \(1-\bar\alpha_t\); their covariance is \(\sqrt{\alpha_t}(1-\bar\alpha_{t-1})\). Gaussian conditioning gives variance \((1-\bar\alpha_{t-1}) - \alpha_t(1-\bar\alpha_{t-1})^2/(1-\bar\alpha_t)\). Since \(\alpha_t(1-\bar\alpha_{t-1}) = \alpha_t - \bar\alpha_t\), this is \((1-\bar\alpha_{t-1})\big[(1-\bar\alpha_t)-(\alpha_t-\bar\alpha_t)\big]/(1-\bar\alpha_t) = \tilde\beta_t\). The mean is \(\sqrt{\bar\alpha_{t-1}}x_0 + \frac{\sqrt{\alpha_t}(1-\bar\alpha_{t-1})}{1-\bar\alpha_t}(x_t - \sqrt{\bar\alpha_t}x_0)\); with \(\sqrt{\bar\alpha_t} = \sqrt{\alpha_t}\sqrt{\bar\alpha_{t-1}}\), the \(x_0\) coefficient becomes \(\sqrt{\bar\alpha_{t-1}}\big[1 - \frac{\alpha_t - \bar\alpha_t}{1-\bar\alpha_t}\big] = \frac{\sqrt{\bar\alpha_{t-1}}\beta_t}{1-\bar\alpha_t}\).
2. The coefficient of \(x_t\) becomes \(\frac{\beta_t/\sqrt{\alpha_t} + \sqrt{\alpha_t}(1-\bar\alpha_{t-1})}{1-\bar\alpha_t} = \frac{\beta_t + \alpha_t - \bar\alpha_t}{\sqrt{\alpha_t}(1-\bar\alpha_t)} = \frac{1}{\sqrt{\alpha_t}}\). The coefficient of \(\varepsilon\) is \(-\frac{\sqrt{\bar\alpha_{t-1}}\beta_t\sqrt{1-\bar\alpha_t}}{(1-\bar\alpha_t)\sqrt{\bar\alpha_t}} = -\frac{\beta_t}{\sqrt{\alpha_t}\sqrt{1-\bar\alpha_t}}\).

1. For any function \(g\), \(\mathbb{E}\|\varepsilon - g(x_t)\|^2 = \mathbb{E}\|\varepsilon - \mathbb{E}[\varepsilon\mid x_t]\|^2 + \mathbb{E}\|\mathbb{E}[\varepsilon\mid x_t] - g(x_t)\|^2\); the cross term vanishes by the tower rule.
2. Let \(a=\sqrt{\bar\alpha_t}\), \(s=\sqrt{1-\bar\alpha_t}\). Then \(p_t(x) = \int p(x_0)\,\mathcal N(x; a x_0, s^2 I)\,dx_0\), and differentiating under the integral gives \(\nabla\log p_t(x) = \mathbb{E}\big[-(x - a x_0)/s^2 \mid x_t = x\big] = -\mathbb{E}[\varepsilon\mid x_t = x]/s\).

import numpy as np

# Data: 1-D mixture of two Gaussians
W = np.array([0.3, 0.7]); M = np.array([-2.0, 1.5]); S2 = np.array([0.15, 0.3])

def sample_data(n, rng):
    k = rng.choice(2, size=n, p=W)
    return rng.normal(M[k], np.sqrt(S2[k]))

T = 200
betas = np.linspace(1e-4, 0.05, T)
alphas = 1 - betas
abar = np.cumprod(alphas)

def q_sample(x0, t, rng):
    eps = rng.normal(size=x0.shape)
    return np.sqrt(abar[t]) * x0 + np.sqrt(1 - abar[t]) * eps, eps

def eps_star(x, t):
    """Exact E[eps | x_t = x] for the mixture: -sqrt(1 - abar) * score of p_t."""
    a, s2 = np.sqrt(abar[t]), 1 - abar[t]
    var = a**2 * S2 + s2                          # per-component variance of x_t
    logw = np.log(W) - 0.5 * np.log(var) - 0.5 * (x[:, None] - a * M) ** 2 / var
    r = np.exp(logw - logw.max(1, keepdims=True)); r /= r.sum(1, keepdims=True)
    score = (r * (-(x[:, None] - a * M) / var)).sum(1)
    return -np.sqrt(s2) * score

def ancestral_sample(n, rng):
    x = rng.normal(size=n)
    for t in range(T - 1, -1, -1):
        mean = (x - betas[t] / np.sqrt(1 - abar[t]) * eps_star(x, t)) / np.sqrt(alphas[t])
        if t > 0:
            var = (1 - abar[t - 1]) / (1 - abar[t]) * betas[t]
            x = mean + np.sqrt(var) * rng.normal(size=n)
        else:
            x = mean
    return x

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    x0 = sample_data(200000, rng)

    # (a) closed-form marginal q(x_t | x_0) vs. running the chain step by step
    t = 120
    x = x0[:50000].copy()
    for s in range(t + 1):
        x = np.sqrt(alphas[s]) * x + np.sqrt(betas[s]) * rng.normal(size=x.shape)
    xt, _ = q_sample(x0[:50000], t, rng)
    print(f"chain mean/var {x.mean():.3f}/{x.var():.3f}   closed form {xt.mean():.3f}/{xt.var():.3f}")
    print(f"abar_T = {abar[-1]:.2e}  (x_T is ~ N(0,1))")

    # (b) the exact epsilon-predictor is the regression target: it beats any other function
    xt, eps = q_sample(x0[:50000], 60, rng)
    print(f"MSE exact eps*: {np.mean((eps - eps_star(xt, 60))**2):.4f}   MSE predict 0: {np.mean(eps**2):.4f}")

    # (c) sample and compare to the data distribution
    xs = ancestral_sample(100000, rng)
    frac_left = (xs < 0).mean()
    print(f"samples: mean {xs.mean():.3f} (data {x0.mean():.3f}), var {xs.var():.3f} (data {x0.var():.3f}), P(x<0) {frac_left:.3f} (data {(x0<0).mean():.3f})")

chain mean/var 0.181/1.270   closed form 0.182/1.298
abar_T = 6.12e-03  (x_T is ~ N(0,1))
MSE exact eps*: 0.4899   MSE predict 0: 0.9935
samples: mean 0.446 (data 0.451), var 2.810 (data 2.829), P(x<0) 0.303 (data 0.302)

1. \(r(d+k) = 8\cdot8192 = 65{,}536\) versus \(dk = 16{,}777{,}216\): about 0.39%.
2. The columns of \(BA\) are combinations of the \(r\) columns of \(B\).
3. \(\mathcal L\) depends on \(B, A\) through \(W+BA\), so \(d\mathcal L = \langle G, dB\,A + B\,dA\rangle\), giving \(GA^\top\) and \(B^\top G\). At initialization \(\Delta W = 0\), so the model starts exactly at the pretrained one, and \(\partial\mathcal L/\partial A = 0\). The first step moves only \(B\); after that both move.
4. The loss is \(\tfrac12\mathbb{E}\|(\Delta W^\star - BA)x\|^2 = \tfrac12\|\Delta W^\star - BA\|_F^2\). By the Eckart–Young theorem the minimum over rank-\(r\) matrices is the truncated SVD, and the leftover error is \(\sqrt{\sum_{i\gt r}\sigma_i^2}\).

1. The per-example loss is \(\ell_{ce}(\Theta\phi(x), y)\): a linear map followed by \(-t_y + \log\sum_j e^{t_j}\). Log-sum-exp has Hessian \(\operatorname{diag}(s) - ss^\top\), which is PSD by A2(b), and convexity is preserved under composition with a linear map and under sums.
2. Fine-tuning can reshape features the task needs but the pretrained model does not expose linearly. With few labels it can overfit or distort good pretrained features; the probe has far fewer parameters and a convex problem. LoRA sits between the two.

1. Row \(i\) is a cross-entropy with label \(i\), so A2(c) gives \(P_{i\cdot} - e_i\), scaled by \(1/n\).
2. \(S_{ij}\) depends on \(u_i\) through \(v_j/\tau\). The gradient-descent step moves \(u_i\) toward its positive \(v_i\) and away from the probability-weighted average of all candidates, so negatives the model currently confuses with the positive get pushed hardest.
3. All scores are equal, so \(P_{ii} = 1/n\) and \(\mathcal L = \log n\), the loss of random guessing. Any embedding that separates pairs does better, so collapse is not optimal. This is the role negatives play.
4. The softmax concentrates on the largest score, so the loss approaches a margin between the positive and the single hardest negative.

1. Expand: \(\|u\|^2 + \|v\|^2 - 2u^\top v = 2 - 2u^\top v\), which decreases as \(u^\top v\) increases.
2. Split the corpus into chunks and embed each; embed the query; retrieve the top-\(k\) chunks by cosine similarity; put them in the prompt and generate. If the right chunk is never retrieved (bad chunking, or the query and passage use different vocabulary), the generator cannot recover it no matter how capable it is. Retrieval recall is the ceiling.

import numpy as np

def softmax(z, axis=-1):
    z = z - z.max(axis=axis, keepdims=True); e = np.exp(z)
    return e / e.sum(axis=axis, keepdims=True)

# ---- Part 1: InfoNCE (CLIP-style, one direction) ----
def info_nce(U, V, tau):
    """U, V: (n, d) unit-norm embeddings; positives are matching rows."""
    S = U @ V.T / tau
    P = softmax(S, axis=1)
    loss = -np.mean(np.log(np.diag(P)))
    dS = (P - np.eye(len(U))) / len(U)        # dLoss/dS
    dU = dS @ V / tau
    return loss, dU

# ---- Part 2: LoRA on a linear model ----
def lora_fit(X, Y, W0, r, lr, steps, rng):
    d_out, d_in = W0.shape
    A = rng.normal(size=(r, d_in)) / np.sqrt(d_in)
    B = np.zeros((d_out, r))                    # B = 0  =>  delta W = 0 at init
    n = len(X); hist = []
    for _ in range(steps):
        R = X @ (W0 + B @ A).T - Y              # residuals (n, d_out)
        G = R.T @ X / n                          # dL/dW for L = ||R||^2 / (2n)
        gB, gA = G @ A.T, B.T @ G
        B -= lr * gB; A -= lr * gA
        hist.append(0.5 * np.mean(np.sum(R**2, 1)))
    return B, A, hist

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    n, d, tau = 8, 5, 0.1
    U = rng.normal(size=(n, d)); U /= np.linalg.norm(U, axis=1, keepdims=True)
    V = U + 0.3 * rng.normal(size=(n, d)); V /= np.linalg.norm(V, axis=1, keepdims=True)
    loss, dU = info_nce(U, V, tau)
    eps = 1e-6; num = np.zeros_like(U)
    for i in range(n):
        for j in range(d):
            Up, Um = U.copy(), U.copy(); Up[i, j] += eps; Um[i, j] -= eps
            num[i, j] = (info_nce(Up, V, tau)[0] - info_nce(Um, V, tau)[0]) / (2 * eps)
    print(f"InfoNCE loss {loss:.4f}, grad max err {np.abs(dU - num).max():.2e}, chance level log n = {np.log(n):.4f}")

    d_in, d_out, r_true = 40, 30, 3
    W0 = rng.normal(size=(d_out, d_in)) / np.sqrt(d_in)
    dW_true = rng.normal(size=(d_out, r_true)) @ rng.normal(size=(r_true, d_in)) / np.sqrt(d_in)
    X = rng.normal(size=(2000, d_in)); Y = X @ (W0 + dW_true).T
    # gradient of A is exactly zero at init (B = 0)
    G0 = (X @ W0.T - Y).T @ X / len(X); B0 = np.zeros((d_out, r_true)); A0 = rng.normal(size=(r_true, d_in))
    print("at init: |grad A| =", np.abs(B0.T @ G0).max(), " |grad B| =", round(float(np.abs(G0 @ A0.T).max()), 3))
    for r in [1, 3, 8]:
        B, A, hist = lora_fit(X, Y, W0, r, lr=0.1, steps=4000, rng=rng)
        rel = np.linalg.norm(B @ A - dW_true) / np.linalg.norm(dW_true)
        print(f"r={r}: params {r*(d_in+d_out):4d} vs full {d_in*d_out}; final loss {hist[-1]:.2e}; rel err of dW {rel:.3f}")

InfoNCE loss 0.3675, grad max err 1.01e-10, chance level log n = 2.0794
at init: |grad A| = 0.0  |grad B| = 3.691
r=1: params   70 vs full 1200; final loss 2.85e+01; rel err of dW 0.739
r=3: params  210 vs full 1200; final loss 3.11e-30; rel err of dW 0.000
r=8: params  560 vs full 1200; final loss 2.68e-30; rel err of dW 0.000

The best rank-1 approximation of this \(\Delta W^\star\) has relative error 0.739, which matches the \(r=1\) run.

\(\nabla\int p_\theta R = \int p_\theta\,(\nabla\log p_\theta)\,R\). The log of the product is a sum, and \(\log p(s_0)\) and \(\log P(s_{t+1}\mid s_t,a_t)\) do not depend on \(\theta\), so only the policy terms survive. The estimator needs sampled trajectories and \(\nabla\log\pi\), never a model of the environment.

1. \(b(s)\sum_a \pi(a\mid s)\nabla\log\pi(a\mid s) = b(s)\nabla\sum_a\pi(a\mid s) = b(s)\nabla 1 = 0\). Subtracting a baseline keeps the estimator unbiased and can greatly reduce its variance.
2. \(\log\pi(a) = \theta_a - \log\sum_j e^{\theta_j}\), so the gradient is \(e_a - \pi\). Then \(\sum_a \pi_a\mu_a(e_a - \pi)\) has \(a\)-th component \(\pi_a\mu_a - \pi_a\,\pi^\top\mu\).

1. \(A\gt0\): \(L = rA\) for \(r\le1+\epsilon\), then flat at \((1+\epsilon)A\). \(A\lt0\): \(L\) is flat at \((1-\epsilon)A\) for \(r\le1-\epsilon\), then \(L = rA\) for all larger \(r\), with no cap.
2. Exactly when \(A\gt0, r\gt1+\epsilon\) or \(A\lt0, r\lt1-\epsilon\). Once an update has moved the policy far enough in the helpful direction, PPO stops pushing; a move in the harmful direction is never clipped, so it is always corrected.

1. Mean \(p\), population standard deviation \(\sqrt{p(1-p)}\). Then \((1-p)/\sqrt{p(1-p)} = \sqrt{(1-p)/p}\) and \(-p/\sqrt{p(1-p)} = -\sqrt{p/(1-p)}\). A rare success gets a large positive weight; a rare failure gets a large negative one.
2. The rewards have zero variance, so every advantage is zero (implementations set them to 0 rather than divide by zero) and the prompt contributes no gradient. Prompts the model always or never solves are wasted; training signal comes from prompts at the edge of its ability.

import numpy as np

def softmax(z):
    z = z - z.max(); e = np.exp(z); return e / e.sum()

def grad_log_pi(theta, a):
    g = -softmax(theta); g[a] += 1.0             # d/dtheta log softmax(theta)_a = e_a - pi
    return g

def reinforce_grads(theta, means, n, baseline, rng):
    pi = softmax(theta); gs = []
    for _ in range(n):
        a = rng.choice(len(theta), p=pi)
        r = rng.normal(means[a], 1.0)
        b = pi @ means if baseline else 0.0      # value baseline V = E_pi[r]
        gs.append((r - b) * grad_log_pi(theta, a))
    return np.array(gs)

def group_advantages(rewards):
    s = rewards.std()
    return np.zeros_like(rewards) if s == 0 else (rewards - rewards.mean()) / s

def ppo_clip_grad(theta, theta_old, actions, adv, eps):
    """Gradient of mean_i min(r_i A_i, clip(r_i, 1-eps, 1+eps) A_i)."""
    pi, pi_old = softmax(theta), softmax(theta_old)
    g = np.zeros_like(theta)
    for a, A in zip(actions, adv):
        ratio = pi[a] / pi_old[a]
        clipped = (A > 0 and ratio > 1 + eps) or (A < 0 and ratio < 1 - eps)
        if not clipped:
            g += A * ratio * grad_log_pi(theta, a)   # d r/dtheta = r * grad log pi
    return g / len(actions)

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    means = np.array([1.0, 2.0, 5.0, 5.5]); theta = np.zeros(4)
    for base in [False, True]:
        G = reinforce_grads(theta, means, 20000, base, rng)
        print(f"baseline={base}: mean grad {np.round(G.mean(0), 3)}, total variance {G.var(0).sum():.2f}")
    print("exact grad:", np.round(softmax(theta) * (means - softmax(theta) @ means), 3))

    # Verifiable-reward task: 4 answers, only answer 3 is "correct" (reward 1).
    # GRPO-style loop: sample a group of 8, normalize rewards, take a few PPO-clip steps.
    theta = np.array([1.0, 0.5, 0.0, -1.0])
    print(f"start P(correct) = {softmax(theta)[3]:.3f}")
    for it in range(60):
        pi = softmax(theta); acts = rng.choice(4, size=8, p=pi)
        adv = group_advantages((acts == 3).astype(float))
        theta_old = theta.copy()
        for _ in range(4):
            theta = theta + 0.5 * ppo_clip_grad(theta, theta_old, acts, adv, eps=0.2)
    print(f"end   P(correct) = {softmax(theta)[3]:.3f}")
    for p in [0.125, 0.5, 0.875]:
        k = int(8 * p); r = np.array([1.0] * k + [0.0] * (8 - k))
        A = group_advantages(r)
        print(f"p={p}: A(correct)={A[0]:.3f} vs sqrt((1-p)/p)={np.sqrt((1-p)/p):.3f};  A(wrong)={A[-1]:.3f} vs -sqrt(p/(1-p))={-np.sqrt(p/(1-p)):.3f}")

baseline=False: mean grad [-0.606 -0.344  0.415  0.535], total variance 11.23
baseline=True: mean grad [-0.597 -0.341  0.404  0.534], total variance 2.59
exact grad: [-0.594 -0.344  0.406  0.531]
start P(correct) = 0.064
end   P(correct) = 0.988
p=0.125: A(correct)=2.646 vs sqrt((1-p)/p)=2.646;  A(wrong)=-0.378 vs -sqrt(p/(1-p))=-0.378
p=0.5: A(correct)=1.000 vs sqrt((1-p)/p)=1.000;  A(wrong)=-1.000 vs -sqrt(p/(1-p))=-1.000
p=0.875: A(correct)=0.378 vs sqrt((1-p)/p)=0.378;  A(wrong)=-2.646 vs -sqrt(p/(1-p))=-2.646

Both estimators agree with the exact gradient; the baseline cuts total variance by about 4×.

1. \(V_t(s) = \min_u\, q s^2 + r u^2 + \mathbb{E}\big[p_{t+1}(as + bu + w)^2\big] + c_{t+1}\). Since \(w\) has mean zero and is independent of \(s, u\), the expectation is \(p_{t+1}(as+bu)^2 + p_{t+1}\sigma^2\). The objective is a convex quadratic in \(u\); setting its derivative \(2ru + 2bp_{t+1}(as+bu)\) to zero gives \(u = -\frac{abp_{t+1}}{r+b^2p_{t+1}}s\). Substituting, \(\min_u\big[ru^2 + p_{t+1}(as+bu)^2\big] = \big(a^2p_{t+1} - \frac{(abp_{t+1})^2}{r+b^2p_{t+1}}\big)s^2\), which gives \(p_t\); the constant terms give \(c_t\). The base case is \(p_T = q\), \(c_T = 0\).
2. The recursions for \(p_t\) and \(K_t\) never involve \(\sigma\); noise only adds the constant \(c_t\). The controller is the same one you would use if the noise were zero.
3. See the reference solution below. The open-loop system is unstable (\(a = 1.1\)); the feedback brings the closed-loop factor \(a - bK\) down to about 0.36.

import numpy as np

# Scalar finite-horizon LQR:  s_{t+1} = a s_t + b u_t + w_t,  w_t ~ N(0, sigma^2)
# cost = sum_{t=0}^{T-1} (q s_t^2 + r u_t^2) + q s_T^2
a, b, q, r, sigma, T = 1.1, 0.5, 1.0, 0.2, 0.3, 30

def riccati(a, b, q, r, sigma, T):
    """Backward pass. V_t(s) = p[t] s^2 + c[t];  optimal u_t = -K[t] s_t."""
    p, c, K = np.zeros(T + 1), np.zeros(T + 1), np.zeros(T)
    p[T] = q
    for t in range(T - 1, -1, -1):
        K[t] = a * b * p[t + 1] / (r + b * b * p[t + 1])
        p[t] = q + a * a * p[t + 1] - (a * b * p[t + 1]) ** 2 / (r + b * b * p[t + 1])
        c[t] = c[t + 1] + p[t + 1] * sigma ** 2
    return p, c, K

def rollout_cost(K, s0, n, rng):
    s = np.full(n, s0, dtype=float); cost = np.zeros(n)
    for t in range(T):
        u = -K[t] * s
        cost += q * s * s + r * u * u
        s = a * s + b * u + sigma * rng.normal(size=n)
    return cost + q * s * s

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    p, c, K = riccati(a, b, q, r, sigma, T)
    s0 = 2.0
    print(f"K_0 = {K[0]:.4f}  (steady-state gain; a = {a} is unstable, closed loop a - bK = {a - b*K[0]:.4f})")
    print(f"predicted expected cost p_0 s0^2 + c_0 = {p[0]*s0**2 + c[0]:.4f}")
    print(f"simulated expected cost               = {rollout_cost(K, s0, 400000, rng).mean():.4f}")
    for scale in [0.9, 1.1]:
        print(f"gain x{scale}: simulated cost {rollout_cost(scale * K, s0, 400000, rng).mean():.4f}")
    _, _, K0 = riccati(a, b, q, r, 0.0, T)
    print("gains identical with and without noise:", np.allclose(K, K0))

K_0 = 1.4823  (steady-state gain; a = 1.1 is unstable, closed loop a - bK = 0.3589)
predicted expected cost p_0 s0^2 + c_0 = 10.9992
simulated expected cost               = 10.9975
gain x0.9: simulated cost 11.1052
gain x1.1: simulated cost 11.0978
gains identical with and without noise: True