AI Math (500 Questions)

This deck walks you through the core mathematical building blocks behind modern AI, starting with vectors and the ways you can add, scale, and measure them. You'll get comfortable with different norms like L1, L2, and the general Lp, as well as the dot product and its close cousins, cosine similarity and the Cauchy-Schwarz inequality. From there, it moves into matrices, covering their definition, arithmetic operations, and the rules of matrix multiplication.

The material is a great fit if you're beginning a journey into machine learning or deep learning and want a solid grip on the algebra_essentials">linear algebra that shows up everywhere, from embedding spaces to neural network weights. It's also handy as a refresher for anyone who learned these concepts a while ago and wants to quickly re-establish fluency before tackling more advanced topics like gradients, eigendecomposition, or attention mechanisms.

Because the deck is concept-heavy and full of formulas, try studying it in short, focused sessions rather than long cramming blocks. Spend a moment after each card trying to sketch the idea on paper, even briefly, since visualising vectors and matrices is what really cements these definitions. Spacing your reviews over a few days will also help the terminology and properties feel like second nature rather than memorized definitions.

Vector space structure forms the foundation of all machine learning mathematics. A vector \(\mathbf{x} = [x_1, \ldots, x_n]^\top \in \mathbb{R}^n\) is an ordered list of \(n\) numbers that represents a data point, feature embedding, or parameter set as a point in \(n\)-dimensional space. Vectors support element-wise addition \([a_1+b_1, \ldots, a_n+b_n]^\top\) and scalar multiplication \(c\mathbf{x} = [cx_1, \ldots, cx_n]^\top\), which geometrically stretches, shrinks, or flips the vector. Measuring magnitude uses norms: the L2 norm \(\|\mathbf{x}\|_2 = \sqrt{\sum_i x_i^2}\) is the Euclidean distance from the origin, the L1 norm \(\|\mathbf{x}\|_1 = \sum_i |x_i|\) (Manhattan distance) encourages sparsity in regularization, and the general Lp norm \(\|\mathbf{x}\|_p = (\sum_i |x_i|^p)^{1/p}\) interpolates between them, becoming the max norm as \(p \to \infty\). The dot product \(\mathbf{a} \cdot \mathbf{b} = \sum_i a_i b_i = \|\mathbf{a}\|_2\|\mathbf{b}\|_2 \cos\theta\) measures alignment between vectors—orthogonal iff zero, with cosine similarity \(\cos\theta = (\mathbf{a}\cdot\mathbf{b})/(\|\mathbf{a}\|_2\|\mathbf{b}\|_2)\) giving a normalized \([-1,1]\) score. The Cauchy-Schwarz inequality bounds the dot product by the product of norms: \(|\mathbf{a} \cdot \mathbf{b}| \leq \|\mathbf{a}\|_2 \|\mathbf{b}\|_2\), with equality iff the vectors are parallel.

Matrix operations extend vectors into rectangular arrays. A matrix \(A \in \mathbb{R}^{m \times n}\) has \(m\) rows and \(n\) columns, with \(A_{ij}\) indexing the entry at row \(i\) and column \(j\). Addition is element-wise \((A+B)_{ij} = A_{ij}+B_{ij}\) for same-shape matrices, and scalar multiplication scales every entry. Matrix multiplication \(C = AB\) with \(A \in \mathbb{R}^{m \times p}\) and \(B \in \mathbb{R}^{p \times n}\) produces \(C \in \mathbb{R}^{m \times n}\) via \(C_{ij} = \sum_k A_{ik}B_{kj}\)—each entry is a dot product of a row of A with a column of B, representing composition of linear transformations. This multiplication is non-commutative (\(AB \neq BA\) in general) but associative (\((AB)C = A(BC)\)). The matrix-vector product \(A\mathbf{x} = \sum_j x_j \mathbf{a}_j\) is a linear combination of A's columns weighted by \(\mathbf{x}\)'s components, which is the foundation of linear layers in neural networks.

Special matrices and transformations include the transpose \((A^\top)_{ij} = A_{ji}\) that flips rows and columns with the key identity \((AB)^\top = B^\top A^\top\) (order reverses). A symmetric matrix satisfies \(A = A^\top\) and always has real eigenvalues (covariance matrices and Hessians are symmetric). The identity matrix \(I\) has 1s on the diagonal and 0s elsewhere, serving as the multiplicative identity (\(AI = IA = A\)), while the zero matrix has all zero entries and acts as the additive identity. The Frobenius norm \(\|A\|_F = \sqrt{\sum_{ij} A_{ij}^2}\) applies the L2 norm to a flattened matrix, equivalently \(\sqrt{\mathrm{tr}(A^\top A)}\).

Scalar matrix functions complete the algebraic toolkit. The trace \(\mathrm{tr}(A) = \sum_i A_{ii}\) sums diagonal entries, equals the sum of eigenvalues, and obeys the cyclic property \(\mathrm{tr}(ABC) = \mathrm{tr}(BCA) = \mathrm{tr}(CAB)\) (non-cyclic permutations are not equal in general). The determinant for a 2×2 matrix \(\det\begin{pmatrix}a & b \\ c & d\end{pmatrix} = ad - bc\) is the signed area of the parallelogram formed by column vectors—zero meaning singular. Its absolute value is the volume scaling factor of the linear map, the product rule \(\det(AB) = \det(A)\det(B)\) makes volume scaling multiplicative, and the determinant equals the product of eigenvalues. The inverse \(A^{-1}\) satisfying \(AA^{-1} = A^{-1}A = I\) exists iff \(\det(A) \neq 0\), with \((AB)^{-1} = B^{-1}A^{-1}\); orthogonal matrices (those with \(Q^\top Q = I\)) have inverse equal to their transpose. The Hadamard product \((A \odot B)_{ij} = A_{ij} B_{ij}\) multiplies corresponding elements and appears in attention masks, dropout masks, and LSTM gate computations.

Subspace concepts reveal the structure of linear transformations. The span of vectors \(\{\mathbf{v}_1, \ldots, \mathbf{v}_k\}\) consists of all linear combinations \(\sum_i c_i \mathbf{v}_i\). Vectors are linearly independent iff the only combination giving zero is the trivial one (all coefficients zero). A basis is a linearly independent set that spans the space—any basis of \(\mathbb{R}^n\) has exactly \(n\) vectors. The rank \(\mathrm{rank}(A) = \dim(\mathrm{col}(A)) = \dim(\mathrm{row}(A))\) counts linearly independent columns (or rows), equals the number of nonzero singular values, and is full iff \(\mathrm{rank} = \min(m,n)\). The column space \(\mathrm{col}(A) = \{A\mathbf{x}\}\) contains all outputs of A, and \(Ax = b\) has a solution iff \(b\) lies in \(\mathrm{col}(A)\). The null space (kernel) \(\ker(A) = \{\mathbf{x}: A\mathbf{x} = \mathbf{0}\}\) contains all inputs A maps to zero, and the rank-nullity theorem states \(\mathrm{rank}(A) + \mathrm{nullity}(A) = n\)—the input dimensions split into mapped and collapsed parts. Rank-deficient matrices have \(\mathrm{rank}(A) < \min(m,n)\), implying singular columns and \(\det(A) = 0\).

Eigenvalues and eigenvectors capture how matrices act in special directions. The eigenvalue equation \(A\mathbf{v} = \lambda\mathbf{v}\) with \(\mathbf{v} \neq \mathbf{0}\) finds directions where A only stretches (by \(\lambda\)) without rotating—positive \(\lambda\) keeps direction, negative flips, zero collapses. To find eigenvalues, solve the characteristic polynomial \(\det(A - \lambda I) = 0\), giving exactly \(n\) eigenvalues (counting multiplicity). Distinct eigenvalues always produce linearly independent eigenvectors. The eigenspace \(E_\lambda = \ker(A - \lambda I)\) collects all eigenvectors for \(\lambda\) plus zero; its dimension is the geometric multiplicity. Symmetric matrices have real eigenvalues with orthogonal eigenvectors (spectral theorem). Eigendecomposition factors \(A = Q\Lambda Q^{-1}\) where Q holds eigenvectors and \(\Lambda\) holds eigenvalues; for symmetric matrices \(A = Q\Lambda Q^\top\) with orthogonal Q. This enables easy computation of matrix powers \(A^k = Q\Lambda^k Q^{-1}\) and matrix exponentials \(e^A = Q e^\Lambda Q^{-1}\). The trace equals the sum and determinant equals the product of eigenvalues—quick sanity checks on computations.

Matrix decompositions generalize eigendecomposition. The Singular Value Decomposition (SVD) factors any matrix \(A = U\Sigma V^\top\) with orthogonal U, V and diagonal \(\Sigma\) of non-negative singular values \(\sigma_i \geq 0\). These relate to eigenvalues via \(\sigma_i(A) = \sqrt{\lambda_i(A^\top A)}\), measuring how much A stretches space in each direction. The economy SVD keeps only nonzero singular values. By the Eckart-Young theorem, the best rank-k approximation is \(A_k = \sum_{i=1}^k \sigma_i \mathbf{u}_i \mathbf{v}_i^\top\) with reconstruction error \(\|A - A_k\|_F = \sqrt{\sum_{i>k} \sigma_i^2}\)—this drives dimensionality reduction, recommendation systems via latent factors \(\hat{R} \approx U_k \Sigma_k V_k^\top\), and Principal Component Analysis (PCA), which projects onto the k directions of maximum variance. A symmetric positive definite matrix can be Cholesky-factored \(A = LL^\top\) with lower-triangular L, useful for efficient linear system solving. The QR decomposition \(A = QR\) with orthonormal columns in Q and upper-triangular R comes from the Gram-Schmidt process. Positive definite matrices \(A \succ 0\) have all \(\lambda_i > 0\) (bowl-shaped curvature, unique minimum); positive semi-definite (PSD) matrices \(A \succeq 0\) have \(\lambda_i \geq 0\) (covariance matrices are always PSD, with precision matrix \(\Lambda = \Sigma^{-1}\) encoding conditional independence). The Moore-Penrose pseudoinverse \(A^+ = V\Sigma^+ U^\top\) generalizes the inverse to non-square or singular matrices, providing minimum-norm least-squares solutions. Additional matrix norms include the spectral norm \(\|A\|_2 = \sigma_{\max}(A)\) (maximum stretching, used in GAN spectral normalization) and nuclear norm \(\|A\|_* = \sum_i \sigma_i\) (convex relaxation of rank, used in matrix completion).

Derivatives measure instantaneous rates of change, forming the foundation of optimization in machine learning. The derivative \(f'(x) = \lim_{h \to 0}(f(x+h)-f(x))/h\) is the slope of the tangent line, equivalently written as \(df/dx\) or \(\dot{f}\) in different notations. Basic rules include the constant rule (\(d/dx[c] = 0\)), power rule (\(d/dx[x^n] = nx^{n-1}\)), sum rule (\(d/dx[f+g] = f'+g'\), reflecting linearity), product rule (\(d/dx[fg] = f'g + fg'\)), and quotient rule (\(d/dx[f/g] = (f'g - fg')/g^2\)). The chain rule for composite functions \(d/dx[f(g(x))] = f'(g(x)) g'(x)\) extends to multiple compositions \(d/dx[f_1(f_2(\cdots f_n(x)))] = f_1' \cdot f_2' \cdots f_n'\)—backpropagation in deep networks is precisely this chain rule applied repeatedly. Key derivatives include \(d/dx[e^x] = e^x\) (the exponential is its own derivative), \(d/dx[\ln x] = 1/x\), \(d/dx[\sigma(x)] = \sigma(x)(1-\sigma(x))\) for sigmoid, and \(d/dx[\tanh(x)] = 1 - \tanh^2(x)\) for tanh. The log-derivative trick \(\partial \log p(x|\theta)/\partial\theta = (1/p)\partial p/\partial\theta\) underlies MLE and the REINFORCE estimator \(\nabla_\theta \mathbb{E}_{p_\theta}[f(x)] = \mathbb{E}_{p_\theta}[f(x)\nabla_\theta \log p_\theta(x)]\).

Multivariable calculus extends these ideas. Partial derivatives \(\partial f/\partial x_i\) measure change in one direction while holding others constant. The gradient \(\nabla f(\mathbf{x}) = [\partial f/\partial x_1, \ldots, \partial f/\partial x_n]^\top\) is the vector of partial derivatives, pointing in the direction of steepest ascent—so \(-\nabla f\) points downhill (the basis of gradient descent). Key gradient identities include \(\nabla(\mathbf{w}^\top\mathbf{x}) = \mathbf{w}\) for linear functions and \(\nabla(\mathbf{x}^\top A \mathbf{x}) = (A + A^\top)\mathbf{x}\), which simplifies to \(2A\mathbf{x}\) for symmetric A (analogous to the scalar derivative of \(ax^2\)). The gradient of the squared L2 norm is \(\nabla\|\mathbf{x}\|_2^2 = 2\mathbf{x}\). The directional derivative \(D_\mathbf{u} f = \nabla f \cdot \mathbf{u}\) gives the rate of change in any unit direction, maximized when aligned with \(\nabla f\). The Jacobian matrix \(J_{ij} = \partial f_i/\partial x_j\) generalizes the gradient to vector-valued functions, representing the best linear approximation near a point: \(\mathbf{f}(\mathbf{x}+\delta) \approx \mathbf{f}(\mathbf{x}) + J\delta\). The Jacobian determinant's absolute value gives the volume scaling factor of the transformation—central in change-of-variable formulas for probability densities and normalizing flows \(p_x(\mathbf{x}) = p_z(\mathbf{f}(\mathbf{x}))|\det J|\). The multivariate chain rule expresses composition via Jacobian multiplication: \(\partial \mathbf{f}/\partial \mathbf{x} = (\partial \mathbf{f}/\partial \mathbf{g})(\partial \mathbf{g}/\partial \mathbf{x})\)—exactly backpropagation.

Second-order information describes curvature. The Hessian \(H_{ij} = \partial^2 f/\partial x_i \partial x_j\) is the symmetric matrix of second partial derivatives. At a critical point \(\nabla f(\mathbf{x}^*) = \mathbf{0}\): positive definite Hessian (\(H \succ 0\), all eigenvalues positive) implies a local minimum (bowl shape), negative definite implies local maximum, and indefinite Hessian (mixed signs) implies a saddle point—common in high-dimensional loss landscapes. The Hessian condition number \(\kappa(H) = \lambda_{\max}/\lambda_{\min}\) characterizes how much curvature differs across directions; high values create narrow valleys where gradient descent oscillates. Taylor expansion provides local approximations: \(f(\mathbf{x}+\boldsymbol{\delta}) \approx f(\mathbf{x}) + \nabla f^\top \boldsymbol{\delta}\) (first-order, used to show \(\mathcal{L}(\theta-\eta\nabla\mathcal{L}) \approx \mathcal{L}(\theta) - \eta\|\nabla\mathcal{L}\|^2\)) and the second-order form \(f(\mathbf{x}+\boldsymbol{\delta}) \approx f(\mathbf{x}) + \nabla f^\top \boldsymbol{\delta} + \frac{1}{2}\boldsymbol{\delta}^\top H \boldsymbol{\delta}\) (basis for Newton's method \(\mathbf{x}_{t+1} = \mathbf{x}_t - H^{-1}\nabla f(\mathbf{x}_t)\)). A function is convex iff its Hessian is PSD everywhere—the second-order condition for convexity. L-smooth functions satisfy \(\|\nabla f(\theta_1) - \nabla f(\theta_2)\| \leq L\|\theta_1-\theta_2\|\) (bounded gradient change). Lipschitz continuous functions \(|f(x) - f(y)| \leq L|x-y|\) have bounded slope. Subgradients generalize gradients to non-smooth functions like ReLU and \(|x|\): at the kink of \(|x|\) at zero, any value in \([-1,1]\) is a valid subgradient—this is what allows L1 regularization to produce sparsity. The mean value theorem states some point achieves the average slope; L'Hôpital's rule \(\lim f/g = \lim f'/g'\) handles indeterminate forms \(0/0\) or \(\infty/\infty\).

Unconstrained optimization finds minima of \(f(\theta)\). Convex problems \(f(\lambda\mathbf{x} + (1-\lambda)\mathbf{y}) \leq \lambda f(\mathbf{x}) + (1-\lambda)f(\mathbf{y})\) have no local minima that aren't global—gradient descent is guaranteed to converge. Strongly convex functions satisfy \(\mathcal{L}(\theta_2) \geq \mathcal{L}(\theta_1) + \nabla\mathcal{L}(\theta_1)^\top(\theta_2-\theta_1) + \frac{m}{2}\|\theta_2-\theta_1\|^2\) with constant \(m > 0\), yielding a unique global minimum and linear convergence rate \(\|\theta_t-\theta^*\|^2 \leq (1-m/L)^t\|\theta_0-\theta^*\|^2\). For L-smooth convex functions, GD achieves \(\mathcal{L}(\theta_T) - \mathcal{L}(\theta^*) \leq L\|\theta_0-\theta^*\|^2/(2T)\) (sublinear \(O(1/T)\) rate). Jensen's inequality \(f(\mathbb{E}[X]) \leq \mathbb{E}[f(X)]\) for convex \(f\) underpins many information-theoretic bounds. Constrained optimization adds equality (\(g_i(\theta) = 0\)) and inequality (\(h_j(\theta) \leq 0\)) constraints. Lagrange multipliers satisfy \(\nabla f = \lambda \nabla g\) at the optimum, with the Lagrangian \(\mathcal{L}(\mathbf{x},\lambda) = f(\mathbf{x}) + \lambda g(\mathbf{x})\). The KKT conditions extend this with complementary slackness \(\mu_i h_i(\mathbf{x}^*) = 0\): either the constraint is active (\(h_i = 0\)) or the multiplier is zero (\(\mu_i = 0\)). The Lagrangian dual function \(d(\lambda,\mu) = \min_\theta \mathcal{L}(\theta,\lambda,\mu)\) provides lower bounds; strong duality holds under Slater's condition.

Gradient descent \(\theta_{t+1} = \theta_t - \eta\nabla\mathcal{L}(\theta_t)\) takes steps opposite the gradient. The learning rate \(\eta\) is critical: too large diverges, too small is slow; optimal \(\eta \approx 1/L\). Each step decreases the loss proportional to \(\eta\|\nabla\mathcal{L}\|^2\) for small \(\eta\) by the first-order Taylor approximation. Stochastic Gradient Descent (SGD) uses single random examples—much cheaper per step but noisy. Mini-batch SGD averages B examples: \(\theta_{t+1} = \theta_t - (\eta/B)\sum_{i\in\mathcal{B}} \nabla\ell_i\). Gradient noise (variance \(\sigma^2_{\mathrm{grad}}/B\)) acts as implicit regularization, helping escape sharp minima and saddle points. For convergence, SGD requires Robbins-Monro conditions \(\sum_t \eta_t = \infty\) (enough progress) and \(\sum_t \eta_t^2 < \infty\) (steps decrease). Momentum accelerates convergence in narrow valleys: \(v_{t+1} = \beta v_t + \nabla\mathcal{L}\) with \(\beta \approx 0.9\)—consistent gradient directions accumulate, with effective learning rate \(\approx \eta/(1-\beta)\). Nesterov momentum evaluates the gradient at the look-ahead position \(\theta_t - \beta\eta v_t\) for slightly faster convergence. The subgradient method extends to non-smooth objectives using \(\mathbf{g}_t \in \partial\mathcal{L}(\theta_t)\) with diminishing step sizes.

Adaptive methods scale learning rates per parameter. AdaGrad accumulates squared gradients \(G_t = \sum g_k^2\) and divides \(\eta\) by \(\sqrt{G_t}\)—great for sparse features but LR eventually vanishes to zero. RMSprop uses exponential moving average \(v_t = \rho v_{t-1} + (1-\rho)g_t^2\) with \(\rho \approx 0.9\), fixing AdaGrad's decay issue; effective learning rate \(\eta/\sqrt{v_t}\) normalizes by recent gradient scale. Adam combines first moment \(m_t = \beta_1 m_{t-1} + (1-\beta_1)g_t\) (momentum) with second moment \(v_t = \beta_2 v_{t-1} + (1-\beta_2)g_t^2\) (RMSprop), applies bias correction \(\hat{m}_t = m_t/(1-\beta_1^t)\), and updates \(\theta_{t+1} = \theta_t - \eta\hat{m}_t/(\sqrt{\hat{v}_t}+\epsilon)\) with \(\epsilon \approx 10^{-8}\) for numerical stability. AdamW applies decoupled weight decay \(\lambda\theta\) directly, improving generalization. LAMB adds per-layer trust ratios for large-batch distributed training. Learning rate schedules include cosine annealing \(\eta_t = \eta_{\min} + \frac{1}{2}(\eta_{\max}-\eta_{\min})(1+\cos(\pi t/T))\), linear warmup, step decay, and the 1-cycle policy.

Second-order methods use curvature. Newton's method \(\theta_{t+1} = \theta_t - H^{-1}\nabla\mathcal{L}\) achieves quadratic convergence but costs \(O(n^2)\) memory and \(O(n^3)\) inversion—impractical for large networks. L-BFGS approximates \(H^{-1}\) using the last m gradient differences (\(s_k = \theta_k-\theta_{k-1}\), \(y_k = \nabla\mathcal{L}_k - \nabla\mathcal{L}_{k-1}\)), reducing to \(O(mn)\) memory. The natural gradient \(\tilde\nabla_\theta \mathcal{L} = \mathcal{I}(\theta)^{-1}\nabla_\theta \mathcal{L}\) accounts for parameterization geometry via the Fisher information matrix \(\mathcal{I}(\theta) = \mathbb{E}[\nabla\log p \cdot (\nabla\log p)^\top]\) (TRPO, natural policy gradient); K-FAC approximates this with Kronecker-factored terms \(A_{l-1} \otimes G_l\). Projected gradient descent handles simple constraints via \(\theta_{t+1} = \Pi_C(\theta_t - \eta\nabla f(\theta_t))\). Proximal methods use \(\mathrm{prox}_{\eta g}(\theta) = \arg\min_u \frac{1}{2}\|u-\theta\|^2 + \eta g(u)\)—for L1 this is soft-thresholding \(\mathrm{sign}(\theta_j)\max(|\theta_j|-\eta\lambda, 0)\). Gradient clipping \(\mathbf{g} \leftarrow \tau\mathbf{g}/\max(\|\mathbf{g}\|,\tau)\) prevents exploding gradients. Stochastic gradient Langevin dynamics adds Langevin noise \(\sqrt{2\eta}\,\boldsymbol{\epsilon}_t\) to SGD for Bayesian sampling. Gradient accumulation simulates large batches; the linear scaling rule multiplies \(\eta\) by \(B'/B\) when batch size changes. The implicit bias of GD on linear networks finds minimum nuclear-norm solutions among zero-error solutions.

Probability axioms form the foundation. Kolmogorov axioms: \(P(\Omega) = 1\), \(P(A) \geq 0\), and \(P(A \cup B) = P(A) + P(B)\) for disjoint \(A, B\). The complement rule \(P(A^c) = 1 - P(A)\) and inclusion-exclusion \(P(A \cup B) = P(A) + P(B) - P(A \cap B)\) extend these. Conditional probability \(P(A|B) = P(A \cap B)/P(B)\) restricts the sample space; independence \(P(A \cap B) = P(A)P(B)\) means \(B\) tells nothing about \(A\). The law of total probability marginalizes: \(P(A) = \sum_i P(A|B_i)P(B_i)\). Bayes' theorem \(P(\theta|x) = P(x|\theta)P(\theta)/P(x)\) updates prior beliefs: posterior \(\propto\) likelihood \(\times\) prior. The evidence (marginal likelihood) \(P(x) = \int P(x|\theta)P(\theta)d\theta\) is often intractable—requiring variational inference or MCMC. With i.i.d. observations, the Bayesian update becomes \(P(\theta|x_1,\ldots,x_n) \propto (\prod_i P(x_i|\theta)) P(\theta)\). Conditional independence \(X \perp Y|Z\) means \(P(X,Y|Z) = P(X|Z)P(Y|Z)\)—fundamental in graphical models and Bayesian networks. The Markov property \(P(X_{t+1}|X_t, X_{t-1},\ldots,X_1) = P(X_{t+1}|X_t)\) states the future is independent of the past given the present.

Key distributions cover discrete and continuous cases. Bernoulli: \(P(X=1) = p\), \(P(X=0) = 1-p\), max variance at \(p=0.5\)—foundation for binary classification. Binomial: sum of n Bernoulli trials with \(E[X] = np\), \(\mathrm{Var}(X) = np(1-p)\), approaches Gaussian for large \(n\). Poisson: \(P(X=k) = \lambda^k e^{-\lambda}/k!\) with mean and variance both equal \(\lambda\)—models rare events like word counts in NLP. Geometric: number of trials to first success, memoryless \(P(X>m+n|X>m) = P(X>n)\). Gaussian (Normal): \(\mathcal{N}(x|\mu,\sigma^2) = (1/\sqrt{2\pi\sigma^2})\exp(-(x-\mu)^2/(2\sigma^2))\)—the most important distribution, with maximum entropy under fixed mean/variance and closed under linear operations. The 68-95-99.7 rule gives standard-deviation-based probability mass. The standard normal \(Z = (X-\mu)/\sigma \sim \mathcal{N}(0,1)\) normalizes any Gaussian. Multivariate Gaussian: \(\mathcal{N}(\mathbf{x}|\boldsymbol{\mu},\Sigma) = (2\pi)^{-d/2}|\Sigma|^{-1/2}\exp(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu}))\), with the squared Mahalanobis distance \(d_M(\mathbf{x},\boldsymbol{\mu})^2 = (\mathbf{x}-\boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu})\) measuring deviations. Categorical extends Bernoulli to K outcomes; Dirichlet is the conjugate prior over probability vectors; Beta is conjugate to Bernoulli/Binomial with \(\mathbb{E}[p] = \alpha/(\alpha+\beta)\). A diagonal covariance matrix implies uncorrelated features. The entropy of a Gaussian is \(\frac{1}{2}\ln(2\pi e \sigma^2)\) nats, increasing with \(\sigma\).

Estimation and inference connect distributions to data. Maximum Likelihood Estimation (MLE) finds \(\hat\theta = \arg\max_\theta \prod_i p(x_i|\theta) = \arg\max_\theta \sum_i \log p(x_i|\theta)\)—log-likelihood is more numerically stable (products cause underflow) but equivalent by monotonicity of log. For a Gaussian: \(\hat\mu = \frac{1}{N}\sum_i x_i\) (sample mean) and \(\hat\sigma^2 = \frac{1}{N}\sum_i(x_i-\hat\mu)^2\) (biased; divide by N-1 for unbiased). Maximum A Posteriori (MAP) incorporates prior: \(\hat\theta = \arg\max[\sum_i \log p(x_i|\theta) + \log p(\theta)]\). A Gaussian prior yields L2 regularization; Laplace prior yields L1—making regularization equivalent to MAP with a chosen prior. Conjugate priors keep the posterior in the same family: Beta-Bernoulli \(\theta \sim \mathrm{Beta}(\alpha,\beta), X|\theta \sim \mathrm{Bern}(\theta) \Rightarrow \theta|X \sim \mathrm{Beta}(\alpha+\mathrm{heads}, \beta+\mathrm{tails})\), Dirichlet-Categorical, Gaussian-Gaussian. Sufficient statistics capture all parameter-relevant information. The moment generating function \(M_X(t) = \mathbb{E}[e^{tX}]\) encodes all moments. The Law of Large Numbers gives \(\bar{X}_N \to \mathbb{E}[X]\); the Central Limit Theorem gives \(\sqrt{N}(\bar{X}_N - \mu) \to \mathcal{N}(0,\sigma^2)\), with \(\bar{X}_N \approx \mathcal{N}(\mu, \sigma^2/N)\) in practice. Monte Carlo estimation approximates \(\mathbb{E}_p[f(x)] \approx \frac{1}{N}\sum_i f(x_i)\) via sampling; importance sampling reweights samples from a tractable proposal \(q\): \(\mathbb{E}_p[f(x)] \approx \frac{1}{N}\sum_i (p(x_i)/q(x_i))f(x_i)\). The bias-variance decomposition \(\mathbb{E}[(\hat f(x) - f(x))^2] = \mathrm{Bias}^2 + \mathrm{Variance} + \sigma^2\) characterizes prediction error, with the classic bias-variance tradeoff (complex models: low bias, high variance; simple: vice versa). The Cramér-Rao lower bound \(\mathrm{Var}(\hat\theta) \geq 1/\mathcal{I}(\theta)\) sets the best achievable precision. Gaussian mixture models \(p(x) = \sum_k \pi_k \mathcal{N}(x|\mu_k, \Sigma_k)\) are fit by the EM algorithm, alternating E-step (posterior \(q^{(t)}(z) = p(z|x,\theta^{(t)})\)) and M-step (\(\theta^{(t+1)} = \arg\max_\theta \mathbb{E}_{q^{(t)}}[\log p(x,z|\theta)]\)). Bootstrap resampling estimates sampling distributions without analytical derivations. The predictive distribution \(p(x_{\mathrm{new}}|x_{1:N}) = \int p(x_{\mathrm{new}}|\theta)p(\theta|x_{1:N})d\theta\) averages predictions over all plausible models, more robust than point estimates.

Self-information (surprisal) \(I(x) = -\log_2 p(x)\) measures surprise: rare events carry more information. Doubling rarity adds exactly 1 bit, since log converts products to sums, making information additive for independent events: \(I(x_1 \cap x_2) = I(x_1) + I(x_2)\). Units vary by log base: \(-\log_2\) gives bits (shannons), \(-\ln\) gives nats, \(-\log_{10}\) gives hartleys, with 1 nat ≈ 1.44 bits. Shannon entropy \(H(X) = -\sum_x p(x)\log_2 p(x)\) is the expected self-information—average surprise. It is maximized at \(\log_2 K\) bits when all K outcomes are equally likely (maximum uncertainty) and is zero when the outcome is certain. Binary entropy \(H(p) = -p\log_2 p - (1-p)\log_2(1-p)\) reaches 1 bit at \(p=0.5\) and zero at \(p \in \{0,1\}\)—used in decision tree splitting criteria. Differential entropy \(h(X) = -\int p(x)\ln p(x)dx\) for continuous variables can be negative and is reparameterization-dependent; for a Gaussian: \(h = \frac{1}{2}\ln(2\pi e\sigma^2)\) nats. The maximum entropy principle chooses the least-committed distribution subject to moment constraints: with mean and variance fixed, the result is the Gaussian; with no constraints, the uniform distribution.

Joint and conditional entropy decompose uncertainty. Joint entropy \(H(X,Y) = -\sum_{x,y} p(x,y)\log p(x,y)\) decomposes via the chain rule \(H(X_1,\ldots,X_n) = \sum_i H(X_i|X_1,\ldots,X_{i-1})\), with \(H(X,Y) = H(Y) + H(X|Y)\) for two variables. Conditional entropy \(H(X|Y) = -\sum_{x,y} p(x,y)\log p(x|y) = H(X,Y) - H(Y)\) is the remaining uncertainty in X after observing Y; it is always \(\leq H(X)\) with equality iff \(X \perp Y\). Mutual information \(I(X;Y) = H(X) - H(X|Y) = H(X)+H(Y)-H(X,Y)\) measures statistical dependence—zero iff independent. Equivalently, \(I(X;Y) = D_{\mathrm{KL}}(p(x,y)\|p(x)p(y))\) measures deviation from the independence product. Conditional mutual information \(I(X;Y|Z) = H(X|Z) - H(X|Y,Z)\) is zero iff \(X \perp Y|Z\)—used in graphical models and feature selection. The data processing inequality states \(X \to Y \to Z\) implies \(I(X;Z) \leq I(X;Y)\)—processing can't add information. Shannon's source coding theorem sets the minimum average code length at \(H(X)\) bits; Huffman coding achieves within 1 bit \(\bar{L}_{\mathrm{Huffman}} < H(X) + 1\). Channel capacity \(C = \max_{p(x)} I(X;Y)\) is the maximum reliable transmission rate. Pointwise mutual information \(\mathrm{PMI}(x,y) = \log(p(x,y)/p(x)p(y))\) measures word co-occurrence (foundation of word2vec's PPMI matrices). Normalized mutual information \(\mathrm{NMI}(X,Y) = I(X;Y)/\sqrt{H(X)H(Y)} \in [0,1]\) suits comparing different variable pairs (clustering evaluation). The information bottleneck principle \(\min I(X;Z) - \beta I(Z;Y)\) seeks minimal sufficient statistics.

KL divergence \(D_{\mathrm{KL}}(P\|Q) = \sum_x p(x)\log(p(x)/q(x))\) measures the cost of using Q to encode samples from P. It is non-negative (by Jensen's inequality) and zero iff \(P = Q\), but asymmetric: \(D_{\mathrm{KL}}(P\|Q) \neq D_{\mathrm{KL}}(Q\|P)\). Forward KL (P‖Q) is mean-seeking (Q must cover all of P); reverse KL (Q‖P) is mode-seeking (Q concentrates on P's highest modes—used in VAEs, \(\mathcal{L}_{\mathrm{VAE}} = \mathbb{E}_{q_\phi}[\log p_\theta(x|z)] - D_{\mathrm{KL}}(q_\phi(z|x)\|p(z))\)). The closed-form for two Gaussians is \(\log(\sigma_2/\sigma_1) + (\sigma_1^2+(\mu_1-\mu_2)^2)/(2\sigma_2^2) - \frac{1}{2}\); in VAEs with \(q(z|x) = \mathcal{N}(\mu,\sigma^2)\) and \(p(z) = \mathcal{N}(0,I)\), the KL becomes \(-\frac{1}{2}\sum_j(1+\log\sigma_j^2 - \mu_j^2 - \sigma_j^2)\). The Jensen-Shannon divergence \(\mathrm{JSD}(P\|Q) = \frac{1}{2}D_{\mathrm{KL}}(P\|M) + \frac{1}{2}D_{\mathrm{KL}}(Q\|M)\) with \(M = (P+Q)/2\) is symmetric and bounded by \(\log 2\); the original GAN training minimizes JSD \(\min_G \max_D V(D,G) \equiv -\log 4 + 2\mathrm{JSD}(p\|q)\). Cross-entropy \(H(P,Q) = -\sum_x p(x)\log q(x)\) is the expected code length under Q when sampling from P, always \(\geq H(P)\). Cross-entropy decomposes as \(H(P,Q) = H(P) + D_{\mathrm{KL}}(P\|Q)\)—for one-hot labels, \(H(P) = 0\), so cross-entropy loss exactly equals KL divergence, making MLE equivalent to minimizing cross-entropy and KL. For multi-class classification \(\mathcal{L} = -\sum_k y_k \log \hat y_k\); with softmax outputs, gradient simplifies to \(\hat y - y\). Binary cross-entropy \(\mathcal{L} = -[y\log\hat y + (1-y)\log(1-\hat y)]\) likewise has gradient \(\hat y - y\). Perplexity \(\mathrm{PPL} = 2^{H(P,Q)}\) is the effective vocabulary size in language models. Label smoothing replaces hard labels with \(y_k^{\mathrm{smooth}} = (1-\epsilon)y_k + \epsilon/K\), preventing overconfidence. The evidence lower bound (ELBO) \(\log p(x) \geq \mathbb{E}_q[\log p(x,z)] - \mathbb{E}_q[\log q(z)]\) bounds the log-evidence from below, maximized in variational inference. The reparameterization trick \(\mathbf{z} = \mu + \sigma \odot \boldsymbol{\epsilon}\) makes sampling differentiable for VAEs.

A single neuron computes \(a = \sigma(\mathbf{w}^\top\mathbf{x} + b)\) where \(\mathbf{w}\) are weights, \(b\) is bias, and \(\sigma\) is a nonlinear activation. A layer in matrix form: \(\mathbf{z} = W\mathbf{x} + \mathbf{b}\), \(\mathbf{a} = \sigma(\mathbf{z})\), with \(W \in \mathbb{R}^{d_{\text{out}}\times d_{\text{in}}}\) containing \(d_{\text{in}} \times d_{\text{out}} + d_{\text{out}}\) parameters. The forward pass of a multi-layer network stacks such layers: \(\mathbf{a}^{(0)} = \mathbf{x}\), \(\mathbf{z}^{(l)} = W^{(l)}\mathbf{a}^{(l-1)}+\mathbf{b}^{(l)}\), \(\mathbf{a}^{(l)} = \sigma(\mathbf{z}^{(l)})\). Depth allows hierarchical feature extraction while width controls per-layer capacity. Activation functions inject nonlinearity—without them, stacking layers is equivalent to one linear layer. Sigmoid \(\sigma(x) = 1/(1+e^{-x})\) squashes to (0,1) but saturates at ±∞ causing vanishing gradients (max derivative 0.25 at x=0). Tanh \((e^x-e^{-x})/(e^x+e^{-x})\) is zero-centered (better gradient flow) but still saturates; derivative \(1 - \tanh^2(x)\). ReLU \(\max(0,x)\) is the default for hidden layers—no saturation for positive inputs (gradient = 1), cheap to compute—but suffers from dying ReLU if pre-activation stays negative. Leaky ReLU uses slope \(\alpha\) for negative inputs (\(\alpha \approx 0.01\)); PReLU learns \(\alpha\). ELU \(\alpha(e^x-1)\) for \(x \leq 0\) saturates smoothly. GELU \(x\Phi(x)\) weights input by probability of being positive under a Gaussian (used in BERT, GPT). Swish/SiLU \(x\sigma(x)\) is smooth, non-monotonic; SwiGLU \(\mathrm{Swish}(\mathbf{x}_1) \odot \mathbf{x}_2\) is used in Llama's FFN. Softmax \(\sigma(\mathbf{z})_i = e^{z_i}/\sum_j e^{z_j}\) produces a probability distribution (numerically stable version subtracts \(\max_j z_j\) before exponentiation); temperature scaling \(\sigma(\mathbf{z}/T)\) sharpens or flattens.

Backpropagation efficiently computes gradients via the chain rule. The error signal propagates as \(\delta^{(l)} = ((W^{(l+1)})^\top \delta^{(l+1)}) \odot \sigma'(z^{(l)})\), and weight gradients are outer products \(\partial\mathcal{L}/\partial W^{(l)} = \delta^{(l)}(a^{(l-1)})^\top\); bias gradients are simply \(\delta^{(l)}\). Backward cost \(\approx 2\times\) forward cost. Key gradient identities: \(\partial\|x\|_2^2/\partial x = 2x\), \(\partial\mathcal{L}/\partial z = \hat y - y\) for both MSE \(\mathcal{L}_{\mathrm{MSE}} = \frac{1}{N}\sum_i (y_i - \hat y_i)^2\) and cross-entropy + softmax \(\partial z = \hat y - y\). Loss functions include MAE \(\frac{1}{N}\sum |y_i - \hat y_i|\) (robust to outliers, non-differentiable at 0), Huber (quadratic for small residuals, linear for large), and hinge \(\max(0, 1-y\hat f)\) for SVMs. The vanishing gradient problem arises because \(\partial\mathcal{L}/\partial W^{(1)} = \prod_{l=1}^{L-1}(W^{(l)})^\top \cdot \delta^{(L)}\)—if each factor has norm \(<1\), products decay exponentially. Exploding gradients occur analogously when norms exceed 1. Solutions include residual connections \(\mathbf{y} = F(\mathbf{x},\{W_i\}) + \mathbf{x}\) (gradient flows through identity shortcut \(\partial\mathcal{L}/\partial\mathbf{x} = \partial\mathcal{L}/\partial\mathbf{y}(1 + \partial F/\partial\mathbf{x})\)), proper initialization, and normalization layers.

Weight initialization matters for signal propagation. Xavier (Glorot) initialization \(W \sim \mathcal{N}(0, 2/(n_{\text{in}}+n_{\text{out}}))\) preserves variance for tanh/sigmoid activations. He initialization \(W \sim \mathcal{N}(0, 2/n_{\text{in}})\) accounts for ReLU killing half the activations. Normalization layers further stabilize training. Batch normalization normalizes over the batch dimension \(\hat x_i = (x_i - \mu_{\mathcal{B}})/\sqrt{\sigma_{\mathcal{B}}^2+\epsilon}\), then rescales with learned \(\gamma, \beta\); at inference, uses running exponential moving averages \(\mu_{\mathrm{running}} = (1-m)\mu_{\mathrm{running}} + m\mu_{\mathcal{B}}\). Layer normalization normalizes over features within each sample (batch-independent, preferred for transformers). Group normalization divides channels into groups; instance normalization is per-sample per-channel (used in style transfer). Spectral normalization \(\bar W = W/\sigma_{\max}(W)\) enforces 1-Lipschitz layers (stabilizes GANs). Regularization prevents overfitting: L2 (weight decay) adds \(\frac{\lambda}{2}\sum w_j^2\) with gradient \(\lambda\theta\) (equivalent to Gaussian prior in MAP); L1 adds \(\lambda\sum|w_j|\) (Laplace prior, encourages sparsity). Dropout randomly zeroes activations \(\tilde a_i = a_i m_i/(1-p)\) with Bernoulli mask \(m_i\); inverted dropout scaling \(1/(1-p)\) makes \(\mathbb{E}[\tilde a] = a\), so no rescaling needed at test time—acts as ensemble training.

Convolutional networks introduce strong inductive biases. A 2D convolution slides a kernel over the input: \((I*K)[i,j] = \sum_{mn} I[i+m,j+n]K[m,n]\)—local connectivity and weight sharing (translation equivariance, parameter-independent of spatial size). Pooling reduces spatial dimensions: max pooling \(\max_{(m,n)} X_{im,jn}\) takes the maximum in each window, average pooling takes the mean. Receptive field grows linearly with depth (kernel size \(k\) gives \(k + (k-1)(l-1)\) at layer \(l\)). Residual (skip) connections \(\mathbf{y} = F(\mathbf{x},\{W_i\}) + \mathbf{x}\) let layers learn residuals; the identity shortcut guarantees non-vanishing gradients. Pre-activation design moves BN+ReLU before weights, giving a clean identity path through the skip connection. The universal approximation theorem guarantees that a single hidden layer with enough neurons can approximate any continuous function \(\sup|\hat f - f| < \epsilon\) to arbitrary precision, but doesn't address learnability or parameter efficiency—deep networks with ReLU create \(O((N/k)^{(k-1)d})\) linear regions, exponentially more than wide shallow networks with the same parameter count. Transposed convolutions provide learnable upsampling in decoders. Spectral norm regularization, knowledge distillation \(\mathcal{L}_{\mathrm{KD}} = T^2 D_{\mathrm{KL}}(\text{softmax}(z_T/T)\|\text{softmax}(z_S/T))\) with temperature \(T\), and the LSTM cell state \(C_t = f_t \odot C_{t-1} + i_t \odot \tilde C_t\) (forget, input, output gates \(f_t, i_t, o_t = \sigma(W[h,x]+b)\) with candidate \(\tilde C_t = \tanh(W_C[h,x]+b_C)\)) are additional architectural primitives. Sharpness-aware minimization (SAM) \(\min_\theta \max_{\|\epsilon\|\leq\rho} \mathcal{L}(\theta+\epsilon)\) explicitly seeks flat minima for better generalization.

The Transformer architecture relies on self-attention. For an input sequence \(X \in \mathbb{R}^{n \times d_{\text{model}}}\), three learned projections produce queries \(Q = XW_Q\), keys \(K = XW_K\), and values \(V = XW_V\), each representing what information a token is looking for, advertising what it contains, and what content to output, respectively. Scaled dot-product attention \(\text{Attention}(Q,K,V) = \text{softmax}(QK^\top/\sqrt{d_k})V\) computes compatibility scores between queries and keys (scaled to prevent softmax saturation), converts them to probabilities via softmax applied row-wise, then weight-averages the values to produce each token's context-aware output. The scaling factor \(1/\sqrt{d_k}\) is needed because \(\mathrm{Var}(q \cdot k) = d_k\)—without scaling, large \(d_k\) pushes softmax toward one-hot, zeroing gradients. The attention weight matrix \(A = \text{softmax}(QK^\top/\sqrt{d_k}) \in \mathbb{R}^{n \times m}\) has row \(i\) as a probability distribution over keys. Multi-head attention runs \(h\) heads in parallel: \(\text{head}_i = \text{Attention}(QW_Q^i, KW_K^i, VW_V^i)\), concatenates results, and projects via output matrix \(W^O \in \mathbb{R}^{hd_v \times d_{\text{model}}}\). With \(d_k = d_v = d_{\text{model}}/h\), total parameters are exactly \(4d_{\text{model}}^2\)—the same as four square matrices regardless of \(h\). Each head operates in its own subspace, capturing diverse patterns simultaneously.

Causal (masked) self-attention enforces autoregression for language models. The mask sets future positions to \(-\infty\) before softmax: \(S_{ij}^{\text{masked}} = S_{ij}\) if \(j \le i\), else \(-\infty\). Implemented as matrix addition \(S^{\text{masked}} = S + M\) with upper-triangular \(-\infty\) mask \(M\), this makes the resulting attention matrix strictly lower-triangular—each token attends only to itself and past tokens. Cross-attention differs by sourcing Q from the decoder and K, V from the encoder: \(Q = H_{\text{dec}}W_Q\), \(K = H_{\text{enc}}W_K\), \(V = H_{\text{enc}}W_V\), allowing decoder queries over encoder representations for source-target alignment. Self-attention differs from a fully connected layer by computing input-dependent weights \(A_{ij} = \text{softmax}(q_i \cdot k_j/\sqrt{d_k})\) that adapt dynamically to each input, rather than using fixed learned parameters. Computational complexity of self-attention is \(\mathcal{O}(n^2 d)\)—quadratic in sequence length is the main bottleneck for long contexts. Flash Attention tiles Q, K, V into SRAM blocks and fuses score computation, softmax, and value aggregation into one GPU kernel, avoiding materializing the full \(n \times n\) score matrix in HBM—reducing memory traffic from \(\mathcal{O}(n^2)\) to \(\mathcal{O}(n)\) and providing significant wall-clock speedups.

Transformers need explicit positional information since attention is permutation-equivariant. Sinusoidal positional encoding assigns \(\text{PE}(\text{pos}, 2i) = \sin(\text{pos}/10000^{2i/d_{\text{model}}})\) and \(\text{PE}(\text{pos}, 2i+1) = \cos(\text{pos}/10000^{2i/d_{\text{model}}})\); pairing sin and cos at each frequency forms 2D rotations enabling relative offset encoding via linear transformations. Position vectors are added to token embeddings: \(\tilde X_i = E_i + \text{PE}(i)\). Learned positional embeddings (BERT) provide flexibility but limit max training length. Relative positional encodings condition attention on offsets: \(S_{ij}^{\text{rel}} = (x_iW_Q)(x_jW_K + r_{i-j}W_K^R)^\top/\sqrt{d_k}\), generalizing better to longer sequences by capturing distance rather than absolute position. The attention masking operation \(A_{ij} = -\infty\) if \(j > i\), else \(z_{ij}\) implements causal masking cleanly—\(-\infty\) becomes 0 after softmax. Weight sharing in self-attention (same \(W_Q, W_K, W_V\) across positions) reduces parameters and enforces equivariance. Knowledge distillation transfers behavior from teacher to student using temperature-softened outputs. The neural tangent kernel \(K(\mathbf{x},\mathbf{x}') = \nabla_\theta f_\theta(\mathbf{x})^\top \nabla_\theta f_\theta(\mathbf{x}')\) describes the infinite-width regime where networks behave like kernel methods.

Each transformer block combines attention with a feed-forward sublayer and residual connections. The FFN \(\text{FFN}(x) = \max(0, xW_1 + b_1)W_2 + b_2\) applies position-wise (independent per token, same matrices reused), typically with \(d_{ff} = 4d_{\text{model}}\) expansion. Pre-LN design normalizes before the sublayer: \(y = x + \text{Sublayer}(\text{LayerNorm}(x))\); the original Post-LN normalized after—Pre-LN improves training stability and is now standard (GPT-2, modern LLMs). Layer normalization \(\text{LayerNorm}(x) = \gamma \odot (x-\mu)/(\sigma+\epsilon) + \beta\) standardizes per token across features (with \(\mu = \frac{1}{d}\sum_i x_i\), \(\sigma\) standard deviation), then rescales with learned \(\gamma, \beta\)—preferred over batch norm in transformers because it's batch-independent. The residual formula \(y = x + \text{Sublayer}(x)\) allows gradients to flow directly to earlier layers without passing through the sublayer. The encoder-decoder Transformer has N encoder layers (bidirectional self-attention + FFN) feeding K, V to M decoder layers (causal self-attention + cross-attention + FFN). The linear classifier head projects the final representation \(\hat{\mathbf{y}} = W_{\text{cls}}\mathbf{h} + \mathbf{b}_{\text{cls}}\) with \(W_{\text{cls}} \in \mathbb{R}^{K \times d}\). During autoregressive decoding, the KV cache \(K_{1:t} = [k_1, \ldots, k_t]\), \(V_{1:t} = [v_1, \ldots, v_t]\) stores past projections to avoid recomputation, reducing per-step complexity to \(\mathcal{O}(td)\) and making inference practical. Gradient checkpointing trades computation for memory—recomputing activations during backward pass reduces memory from \(O(N)\) to \(O(\sqrt{N})\) with ~33% extra compute. Mixed precision training scales the loss \(s \cdot \mathcal{L}\) before backward pass and rescales gradients \(s^{-1}\mathbf{g}\) afterward to prevent FP16 underflow.