Source-linked AI summary
The Elements of Differentiable Programming
Mathieu Blondel, Vincent Roulet
TL;DR
Derivative-free optimization scales poorly to high-dimensional problems, motivating differentiable programming. This book reviews differentiation through programs, control flows, data structures, optimization, and probability, showing that gradients can support efficient optimization and connect soft indexing with Transformer attention.
Problem
Derivative-free optimization suffers from the curse of dimensionality, whereas derivative-based methods can scale to millions or billions of parameters.
Source: Why do we need derivatives? · p. 22Method
The book synthesizes differentiable programming through optimization and probability perspectives, covering control flows, data structures, gradient computation, and stochastic methods.
Source: §5 · pp. 5–6, §6 · p. 6, §12.3.1 · p. 312Results
Reverse-mode differentiation computes gradients with complexity at most a constant factor above function evaluation, while Gaussian-kernel soft indexing recovers Transformer softargmax attention.
Source: §8.8 · p. 231Takeaways & Limitations
Differentiable programming provides a framework for optimizing high-dimensional programs and relating continuous relaxations of data structures to attention mechanisms.
Source: §8.3.3 · p. 219
Abstract
from arXiv · showhide
Artificial intelligence has recently experienced remarkable advances, fueled by large models, vast datasets, accelerated hardware, and, last but not least, the transformative power of differentiable programming. This new programming paradigm enables end-to-end differentiation of complex computer programs (including those with control flows and data structures), making gradient-based optimization of program parameters possible. As an emerging paradigm, differentiable programming builds upon several areas of computer science and applied mathematics, including automatic differentiation, graphical models, optimization and statistics. This book presents a comprehensive review of the fundamental concepts useful for differentiable programming. We adopt two main perspectives, that of optimization and that of probability, with clear analogies between the two. Differentiable programming is not merely the differentiation of programs, but also the thoughtful design of programs intended for differentiation. By making programs differentiable, we inherently introduce probability distributions over their execution, providing a means to quantify the uncertainty associated with program outputs.
5 Control flows … Differentiation
This section surveys differentiable programming constructs for control flows and structured data, then develops the mathematical foundations and practical techniques for differentiation. It covers linear-map formulations, autodiff primitives, higher-order and manifold derivatives, optimization, and extensions to nonsmooth functions.
5 Control flows
- Control-flow coverage includes comparison, inequality, equality, logical operators, continuous extensions, if-else statements, and else-if statements.The section also treats branch and predicate variables alongside continuous relaxations.Source: §5 · pp. 5–6
- Loop constructs are developed through for loops, scan functions, and while loops, including cyclic-graph, unrolled, and Markov-chain perspectives.Source: §5 · p. 6
6 Data structures
- Data-structure coverage focuses on lists and dictionaries, including variable-length operations, soft indexing, kernel-regression relaxations, probability-distribution perspectives, and links to Transformer attention.Source: §6 · p. 6
1.4 How to read this book?
- It reviews fundamental differentiable-programming techniques rather than their applications, supports non-linear reading with prerequisite guidance, and situates the material among prior texts, tutorials, and autodiff software.Source: §1.3 · p. 25, §1.4 · p. 25, §1.5 · pp. 25–26
1.3. Intended audience
- The book is designed as a graduate-level introduction for machine learning readers, assuming familiarity with calculus, linear algebra, probability, and machine learning.Source: §1.3 · p. 25
Differentiation
- Differentiation is presented as foundational to differentiable programming, with linear maps emphasized as its central organizing concept.Source: Differentiation · p. 28
2.1.1 Derivatives
- The derivative provides local linear approximations and descent directions, while continuity is necessary for differentiability and discontinuous functions require special treatment.Source: §2.1.1 · pp. 29–31
2.1.2 Calculus rules
- The chain rule is the cornerstone of differentiation, enabling mechanical differentiation of composite expressions and motivating automatic differentiation.Source: §2.1.2 · pp. 33–34
2.1.3 Leibniz’s notation
- Leibniz notation expresses the chain rule as multiplication of derivative quotients, while infinitesimal variations and dual numbers underlie automatic differentiation algorithms.Source: §2.1.3 · pp. 34–35
2.2.1 Directional derivatives
- For multidimensional functions, directional derivatives measure variation along arbitrary directions and provide the basis for partial derivatives and differentiability.Source: §2.2.1 · p. 35
2.2.2 Gradients
- The gradient collects partial derivatives and represents directional derivatives through a linear form, relying on linearity with respect to the direction.Source: §2.2.2 · pp. 36–37
- Fréchet differentiability requires a linear first-order approximation, whereas Gateaux differentiability via directional derivatives alone is insufficient; subsequent definitions use the Fréchet notion.A function can have directional derivatives in every direction without being Fréchet differentiable because the directional derivative need not be linear or control arbitrary approaches to zero.Source: §2.2.2 · pp. 37–38
Linearity of gradients
- Linearity preserves weighted sums of gradients, and linear maps generalize this structure beyond matrices; in vector spaces, every such map has a matrix representation.A linear map satisfies l[aw + bv] = a · l[w] + b · l[v].Source: Linearity of gradients · p. 39, §2.3.3 · pp. 47–48
Why is the gradient useful?
- Gradients point toward steepest ascent, while negative gradients point toward lower function values and motivate the update wt+1 := wt −γ∇f(wt).The gradient is orthogonal to level sets when nonzero.Source: Why is the gradient useful? · pp. 39–40, Figure 2.3 · p. 40
2.2.3 Jacobians
- Jacobians collect coordinate-wise partial derivatives into a matrix representing the linear directional derivative, generalizing gradients and derivatives to multi-output functions.For element-wise activations such as softplus, the Jacobian is diagonal, so directional derivatives reduce to element-wise products.Source: §2.2.3 · pp. 40–42
Chain rule
- The chain rule composes differentiable functions and multiplies their Jacobians, while its scalar-valued form gives gradients for vector-valued models followed by scalar losses.The rule is illustrated by deriving the gradient of a linear-regression squared-error loss.Source: Chain rule · pp. 43–45
2.3 Linear maps
- Jacobian-vector products and vector-Jacobian products expose the two linear maps encoding infinitesimal input and output variations, which underlie differentiable-program implementations.These maps avoid requiring explicit Jacobian matrices.Source: §2.3 · p. 45
2.3.1 The need for linear maps
- For matrix inputs, operating directly with linear maps exploits structure: differentiating f(W) := Wx in direction V yields Vx at MD cost instead of a dense flattened-Jacobian computation costing M3D2.Flattening produces a mostly zero RM×MD Jacobian and is therefore inefficient.Source: §2.3.1 · pp. 45–46
2.3.2 Euclidean spaces
- Euclidean spaces extend the framework from vectors to matrices and tensors through bases, inner products, induced norms, distances, and convergence notions.Matrices use canonical matrix bases and the Frobenius inner product; tensors use analogous entrywise inner products.Source: §2.3.2 · pp. 46–47
2.3.3 Linear maps and their adjoints
- Linear-map adjoints avoid explicit matrices: l[v] = (b⊤v)a and l∗[u] = (a⊤u)b compute equivalent products efficiently.The adjoint corresponds to the matrix transpose when a linear map is represented by a matrix.Source: §2.3.3 · p. 48
2.3.4 Jacobian-vector products; 2.3.5 Vector-Jacobian products
- JVPs represent directional derivatives as linear maps, while VJPs are their adjoints and recover output-direction gradients.For f : E → F, the JVP maps v to ∂f(w)[v], and the VJP maps u to ∂f(w)∗[u].Source: §2.3.4 · p. 49, §2.3.5 · p. 51
2.3.6 Chain rule using linear maps
- The chain rule composes JVP linear maps, with a scalar-valued final function yielding the corresponding gradient composition.This formulation is designed to exploit implementations of JVPs and VJPs rather than explicit Jacobian matrices.Source: §2.3.6 · p. 52
2.3.7 Functions of multiple inputs (fan-in)
- For multiple inputs, JVPs sum blockwise partial derivatives, while VJPs return the tuple of adjoint partial derivatives for each input.Nested inputs can be grouped for convenience, and semicolon notation treats selected inputs as constants when differentiating.Source: §2.3.7 · pp. 53–55
2.3.8 Functions with multiple outputs (fan-out)
- For multiple outputs, JVPs decompose across output components and VJPs combine output-direction components through the corresponding adjoint maps.These decompositions integrate with the chain rule for functions whose inputs or outputs are products of spaces.Source: §2.3.8 · p. 55
2.3.9 Extensions to non-Euclidean linear spaces; 2.4. Second-order differentiation; 2.4.2 Second directional derivatives; 2.4.3 Hessians
- Directional derivatives generalize to Banach spaces, while gradients and VJPs require inner products and are well-defined in Hilbert spaces through Riesz representation.Second-order differentiation defines curvature through second directional derivatives; these form a bilinear map represented by the Hessian, which is symmetric under twice differentiability.Source: §2.3.9 · p. 56, §2.4.2 · p. 57, §2.4 · p. 59
2.4.4 Hessian-vector products
- Hessian-vector products extend Hessians as bilinear forms on Euclidean spaces and can be implemented in four autodiff ways.The HVP applies the Hessian to a prescribed vector without requiring the full Hessian matrix.Source: §2.4.4 · p. 60
2.5. Higher-order differentiation
- Higher-order derivatives are defined inductively, while higher-order directional derivatives are multilinear forms and multi-output extensions form higher-order Jacobians.Although partial derivatives can be materialized as tensors, applying the multilinear maps to prescribed directions is often more convenient.Source: §2.5.1 · p. 61, §2.5.2 · p. 62, §2.5.3 · p. 62
2.5.4 Taylor expansions
- Taylor expansions approximate differentiable functions by successive directional-derivative terms, with second-order expansions improving locally over first-order ones and supporting finite-difference computation.The framework extends to nth-order expansions for n-times differentiable functions.Source: §2.5.4 · pp. 63–64
2.6 Differential geometry
- Differential geometry extends derivatives from vector spaces to smooth manifolds, which are locally approximated by Euclidean spaces and inherit differentiability from ambient extensions.The resulting framework uses tangent and cotangent spaces to generalize directional derivatives, JVPs, and VJPs.Source: §2.6 · p. 64, §2.6.1 · pp. 64–65
2.6.2 Tangent spaces and pushforward operators
- On manifolds, tangent spaces collect curve derivatives, and the JVP becomes a pushforward mapping tangent vectors from TwM to Tf(w)N.This construction preserves the chain rule by differentiating compositions with curves.Source: §2.6.2 · pp. 65–66, Figure 2.8 · p. 66
2.6.3 Cotangent spaces and pullback operators
- Cotangent spaces contain linear forms, so manifold VJPs act as pullbacks; on Riemannian manifolds, these can be identified with adjoint operators that pull back gradients.For a sphere, the tangent space is a Euclidean space of dimension P −1 consisting of vectors orthogonal to w.Source: §2.6.3 · pp. 67–68
2.7 Generalized derivatives
- For non-continuous functions such as piecewise-constant control flows, almost-everywhere derivatives may be zero and uninformative, motivating surrogate functions.Continuous non-differentiable functions instead admit generalized directional derivatives, gradients, and Jacobians.Source: §2.7 · p. 70
2.7.1 Rademacher’s theorem
- The section also establishes Rademacher’s theorem as the starting point for characterizing locally Lipschitz continuous functions.This provides the entry point for generalized derivative analysis of continuous functions.Source: §2.7.1 · p. 70
- Rademacher’s theorem states that every locally Lipschitz function between Euclidean spaces is differentiable almost everywhere, with nondifferentiability confined to a measure-zero set.This result motivates extending gradients and Jacobians to locally Lipschitz functions.Source: §2.7.1 · p. 70
2.7.2 Clarke derivatives
- Clarke’s nonsmooth analysis (Clarke, 1975) replaces pointwise slopes with neighborhood-based limit-superior derivatives and represents generalized gradients through convex combinations of nearby gradients.Clarke generalized Jacobians extend the construction to vector-valued functions, recover ordinary derivatives for continuously differentiable functions, and support generalized chain rules.Source: §2.7.2 · pp. 71–73
2.8 Summary
- Automatic differentiation computes composition derivatives through linear JVP and VJP maps, avoiding explicit intermediate Jacobians in computation chains and graphs.For scalar loss ℓ composed with vector function f, ∇(ℓ◦f)(w) = ∂f(w)∗∇ℓ(f(w)).Source: §2.8 · pp. 74–75
Probabilistic learning
This section develops probabilistic learning through distribution fundamentals and conditional exponential-family models for classification and regression. It connects maximum entropy and likelihood methods to predictive distributions, familiar loss functions, and differentiable computation graphs for inference.
3.2.1 Negative log-likelihood
- Maximum likelihood fits pλ to i.i.d. observations by maximizing data likelihood, equivalently minimizing negative log-likelihood; for normals, estimates are the sample mean and sample variance.The normal-distribution result identifies the MLEs of µ and σ2 directly from observed data.Source: §3.2.1 · p. 78
3.2.2 Consistency w.r.t. the Kullback-Leibler divergence
- MLE is consistent in Kullback-Leibler divergence: as N →∞, the estimated parameters converge in expectation to parameters minimizing divergence from the true distribution.The argument uses the displayed relationship and the law of large numbers.Source: §3.2.2 · pp. 78–79
3.3.1 Conditional probability distributions
- Probabilistic models map inputs to distribution parameters, providing full outcome distributions and statistics rather than only deterministic predictions.This enables computation of probability, expectation, variance, mode, and related quantities.Source: §3.3.1 · p. 79, §3.3.2 · p. 80
3.3.3 Binary classification
- Bernoulli, categorical, normal, multivariate normal, and Poisson distributions provide probabilistic formulations for binary, multiclass, scalar, multivariate, and integer outcomes.Valid parameters are enforced with sigmoid, softargmax, and exponential output layers; Poisson modeling assumes dispersion ratio one, with generalized Poisson alternatives when inappropriate (Satterthwaite, 1942).Source: §3.3.3 · p. 80, Parameterization using a sigmoid · p. 81, §3.3.4 · p. 82
Recovering well-known loss functions
- Canonical likelihood parameterizations connect neural-network outputs to standard losses: Bernoulli to logistic, categorical to logistic, Gaussian to squared, and Poisson to convex Poisson loss.The Bernoulli, categorical, Gaussian, and Poisson figures explicitly describe these distribution-specific negative log-likelihood losses.Source: Figure 3.1 · p. 81, Figure 3.2 · p. 83, Figure 3.3 · p. 85
3.3.5 Regression
- For regression, probabilistic predictions support CDF-based queries and uncertainty beyond the mean, while multivariate models can represent covariance, including learned positive semidefinite matrices.A diagonal covariance corresponds to independent normal outputs; covariance can be parameterized as SS⊤ to ensure valid structure.Source: §3.3.5 · p. 84, Parameterization · pp. 85–86, §3.3.6 · p. 86
Deterministic vs. probabilistic approaches
- Probabilistic objectives remain continuous and typically differentiable because distribution parameters vary continuously, avoiding the discontinuity of deterministic losses on discrete outputs.This preserves probabilistic inference while making optimization tractable for discrete-output models.Source: Deterministic vs. probabilistic… · p. 88, Figure 3.5 · p. 89
Negative log-likelihood
- Conditional maximum likelihood estimates model parameters from input-output pairs by maximizing conditional data likelihood, equivalently minimizing the negative log-likelihood.This supplies the loss-function formulation used to learn conditional probabilistic predictors.Source: §3.3.8 · p. 88, Negative log-likelihood · p. 90
Recovering well-known loss functions
- Bernoulli, categorical, normal, and Poisson exponential-family parameterizations recover binary logistic, multiclass logistic, squared, and Poisson loss functions, respectively.The Poisson loss is convex in λi and θi for yi ≥0.Source: Recovering well-known loss functions · pp. 90–91
3.4.1 Definition
- The exponential family expresses probability distributions using natural parameters, a base measure, sufficient statistics, and a log-partition function; Bernoulli parameters satisfy θ = logit(π) and π = logistic(θ).The reviewed distributions belong to this family, though the listed examples are non-exhaustive.Source: §3.4.1 · pp. 91–92
3.4.2 The log-partition function
- The log-partition gradient equals the expected sufficient statistic, its Hessian equals the covariance matrix, and in minimal families it bijectively maps natural parameters to achievable means.Source: §3.4.2 · p. 94
3.4.3 Maximum entropy principle
- Maximum-entropy distributions satisfying first-order moment constraints necessarily belong to the exponential family, with canonical parameters equal to the constraints’ Lagrange multipliers.Source: §3.4.3 · p. 95
3.4.2 The log-partition function
- Exponential-family negative log-likelihood is convex in θ, with gradient equal to the residual between model and observed sufficient-statistic expectations, enabling first-moment matching.The log-partition function is convex, and its gradient gives the expected sufficient statistic.Source: §3.4.2 · p. 94, §3.4.4 · p. 96
Training
- Neural-network training minimizes the conditional negative log-likelihood, which is typically nonconvex in w even though the exponential-family loss is convex in θ.Source: Training · p. 97
Inference
- Inference computes expectations through ∇A(f(x,w)), probabilities through pθ(y), and other statistics by converting natural to original distribution parameters.Source: Inference · pp. 97–98
3.4.5 Probabilistic learning with exponential families
- Probabilistic supervised learning maps inputs to distribution parameters rather than outputs directly, then uses the resulting conditional distribution to compute probabilities, means, variances, and other statistics.The model function must produce valid parameters, such as a positive semidefinite covariance matrix for multivariate normals.Source: §3.4.5 · p. 96, §3.5 · p. 98
Differentiable programs · Parameterized programs
This section develops parameterized differentiable programs from computation graphs and core components through residual, recurrent, and Transformer architectures. It covers normalization, positional mechanisms, attention, computational scaling, and encoder-decoder, decoder-only, and encoder-based designs for processing sequences and generating outputs.
Parameterized programs
- Parameterized programs represent neural networks as learnable functions whose parameters are adjusted to fit data by minimizing losses requiring gradients.Feedforward networks map inputs and learnable parameters to outputs through layered intermediate representations.Source: Parameterized programs · p. 100, §4.2 · p. 106
4.1.1 Computation chains
- Computation chains compose functions sequentially, whereas directed acyclic graphs represent programs whose functions can depend on multiple earlier outputs.Chains impose a total order, while DAGs encode dependencies and support topological ordering when nodes are not totally ordered.Source: §4.1.1 · pp. 100–101, §4.1.2 · pp. 101–102
4.1.3 Computer programs as DAGs
- Under pure-function and termination assumptions, a program’s finite function and variable dependencies form a DAG with designated input, output, and parent nodes.Execution evaluates each function after its parents, using a topological order; multiple inputs or outputs can be grouped into tuples.Source: §4.1.3 · p. 103
Executing a program
- Programs execute by evaluating functions in topological order, with each function viewed equivalently as taking a tuple of parent values or multiple individual inputs.This procedure is summarized in Algorithm 4.1.Source: Executing a program · p. 104
Dealing with multiple program inputs or outputs
- The formalism can group multiple inputs or outputs into single tuple-valued variables, while a bipartite alternative separates function and variable nodes and permits explicit multiple outputs.The bipartite representation is akin to directed factor graphs, but the section adopts the simpler shared-node formalism.Source: Dealing with multiple program inputs… · p. 104, Alternative representation: bipartite… · p. 105
4.1.4 Arithmetic circuits
- Arithmetic circuits are DAGs over a field whose functions are addition or multiplication gates, and their complexity is measured by the minimum edge count representing a polynomial.Circuit size counts edges, while polynomial size is the smallest circuit size among representations of that polynomial.Source: §4.1.4 · pp. 105–106
4.2 Feedforward networks
- Feedforward networks are parameterized computation chains whose layers produce intermediate representations, with multilayer perceptrons combining affine transformations and parameter-free nonlinear activations.Convolutional layers generalize the linear transformation by replacing matrix-vector products with convolutions.Source: §4.2 · p. 106, §4.3 · p. 107
4.3.2 Link with generalized linear models
- MLPs include generalized linear models as the one-layer case and jointly learn feature representations and output models at greater depth, while activations provide nonlinear or probabilistic mappings.ReLU is piecewise linear, softplus is a smooth strictly positive approximation, log-sum-exp approaches max as γ →0, and logistic provides a differentiable probability mapping unlike the binary step.Source: §4.3.2 · pp. 107–108, §4.4.1 · p. 108
4.4.4 Probability mappings: argmax and softargmax
- Argmax concentrates probability on one coordinate but is discontinuous, whereas softargmax is an everywhere-differentiable relaxation whose outputs remain strictly inside the probability simplex.Softargmax outputs satisfy π_j ∈ (0, 1), has M − 1 degrees of freedom, and is non-invertible; it is also the gradient of log-sum-exp.Source: Argmax · p. 112, Softargmax · p. 112, §4.5 · p. 113
4.5. Normalization layers
- Normalization layers address difficult optimization caused by wide-ranging activations; batch normalization standardizes across samples, while layer normalization standardizes across features for each sample.Batch normalization needs running training-set statistics at inference, whereas layer normalization is well defined per sample without a training batch.Source: §4.5 · pp. 114–115, §4.5.1 · p. 114, §4.5.2 · pp. 115–116
4.6 Residual neural networks
- Residual blocks add a skip connection to a learned residual, making identity mappings easier to represent and helping remedy vanishing gradients in deep networks.Deeper networks can otherwise be harder to train, causing accuracy to saturate or degrade with depth; residual blocks can represent identity by setting residual parameters to zero.Source: §4.6 · p. 116
4.7 Recurrent neural networks
- Recurrent neural networks maintain a state updated recursively with parameters shared across steps, supporting vector-to-sequence, sequence-to-vector, and aligned or unaligned sequence-to-sequence mappings.These setups cover tasks such as image captioning, sequence classification, part-of-speech tagging, and machine translation.Source: §4.7 · pp. 117, 119
4.7.4 Sequence to sequence (unaligned)
- Unaligned sequence-to-sequence models use an encoder-decoder architecture with an EOS token for variable-length outputs, but the context vector can bottleneck input information.The encoder maps the input sequence to a context vector, and the decoder generates the target sequence, whose length may depend on the input.Source: §4.7.4 · p. 120
4.8.1 Attention
- Attention forms each output as a convex combination of sequence values, with query-key similarities determining weights and masks able to enforce zero or causal attention.Scaled attention commonly divides by √D_k inside softargmax to reduce variance (Vaswani et al., 2017, Footnote 4).Source: §4.8.1 · pp. 121–122
4.8.2 Self-attention
- Self-attention applies row-wise softargmax to feature similarities, while learned projections increase expressiveness and keep weight-matrix sizes independent of sequence length.This allows self-attention to operate on sequences of arbitrary length (Vaswani et al., 2017).Source: §4.8.2 · pp. 122–123
4.8.3 Multi-head attention
- Multi-head attention captures multiple attention patterns, and Transformer layers combine it with residual connections, layer normalization, and token-wise MLPs iterated into depth-K Transformer blocks.Transformer blocks use different parameters at each iteration and can optionally apply a mask to sequence interactions.Source: §4.8.3 · p. 124, §4.8.4 · pp. 124–125, §4.8.5 · p. 125
Number of parameters and computational complexity
- A Transformer layer uses O(D2) parameters, versus O(L2D2) for a standard multi-layer perceptron, with self-attention costing O(L2D) and its final MLP O(LD2).Multi-head attention, MLP, and layer-normalization parameters are depth-specific.Source: Number of parameters and computational… · p. 126, §4.8.5 · p. 126
4.8.6 Token encoding
- Tokens are mapped from vocabulary indices to D-dimensional vectors, while self-attention is permutation-equivariant and therefore requires positional information to represent sequence order.Possible approaches include absolute or relative positional encodings, modifying attention, or preprocessing with an RNN.Source: §4.8.6 · p. 126, §4.8.7 · p. 127
Sinusoidal positional encoding
- Learned positional embeddings assign trainable vectors to positions, whereas sinusoidal encodings use continuous sine and cosine waves with increasing frequencies.Sinusoidal encodings use N := 10000 in the cited construction and organize the values into P ∈ R^L×D.Source: Learned positional encoding · p. 127, Sinusoidal positional encoding · pp. 128–129
Recovering relative positional information
- Sinusoidal positional encoding makes each offset representation a linear function of the original, enabling Transformers to learn attention by relative positions (Zhang et al., 2021).The construction applies a rotation to each sine-cosine coordinate pair to obtain the representation at the shifted position.Source: Recovering relative positional… · pp. 129–130
Rotary positional encoding (RoPE)
- RoPE rotates query and key vectors by position-dependent angles before attention, combining flexibility to sequence length with capture of relative token distance despite being an absolute encoding.For multi-head attention, the same rotational approach is applied independently to each head.Source: Rotary positional encoding (RoPE) · pp. 131–132
Defining a language model
- A decoder-only Transformer defines next-token probabilities from a context window using causal masking, a linear disembedding layer, and soft-argmax probabilities.Weight tying can set W D := (W E)⊤, reducing learned parameters and working well in practice (Press and Wolf, 2016).Source: Defining a language model · pp. 132–133
Training
- Training maximizes corpus log-likelihood with token-wise cross-entropy and teacher forcing, whereas generation samples sequentially from the resulting higher-order Markov chain, using a KV cache to avoid repeated attention computation.Teacher forcing makes training highly parallelizable; sampling stops at EOS or a maximum length.Source: Training · p. 134, Sampling · pp. 134–135
4.8.9 Encoder-only architectures
- Encoder-only BERT pretrains bidirectional representations with masked-token and next-sentence prediction, then finetunes the encoder for classification, tagging, or question answering; encoder-decoder Transformers support sequence-to-sequence tasks.The encoder-decoder architecture applies to machine translation, summarization, and question answering.Source: Pretraining · pp. 135–136, Finetuning · p. 136, §4.8.10 · p. 136
Encoder
- The encoder builds a rich representation Henc by weighting relationships among all input tokens with unmasked multi-head self-attention and separate subcomponent parameters.Its block includes multi-head attention, layer normalization, and an MLP.Source: Encoder · p. 137
Decoder
- The decoder generates Y one token at a time using causal self-attention over prior outputs, followed by cross-attention to the encoded representation Henc.MLP and layer normalization layers are then applied, with decoder-specific weight matrices.Source: Decoder · p. 137
Differences with decoder-only architectures
- Encoder-decoder models support bidirectional input understanding and are suited to seq2seq tasks, whereas decoder-only models use autoregressive next-token prediction, including prompted seq2seq use despite causal masking.Decoder-only architectures prepend the input to the context and can remain remarkably strong.Source: Differences with decoder-only… · pp. 137–138
4.9 Summary
- Programs can be represented as directed acyclic graphs, neural networks are parameterized programs, and feedforward networks are parameterized computation chains.MLPs, ResNets, and CNNs are feedforward-network parametrizations, while variable-length transformer blocks require positional encoding because they are permutation-equivariant.Source: §4.9 · p. 138
Control flows
This section develops differentiable representations of control-flow primitives by relaxing discrete predicates, comparisons, and branch selection into continuous or probabilistic forms. It covers differentiation through logical operators, conditionals, for loops, scans, and while loops, including unrolled and Markov-chain representations.
5.1 Comparison operators
- Hard comparison operators use Heaviside-based Boolean outputs whose derivatives are zero almost everywhere, making them uninformative for gradient backpropagation.The operators remain differentiable almost everywhere but fail to transmit useful gradients.Source: §5.1 · p. 140
5.2.1 Heuristic definition
- Soft inequality operators replace the step function with sigmoids, including logistic and Gaussian CDFs, while sparse sigmoids can reach one at finite separation.The scaling parameter σ controls the relaxation, and sigmoidσ(µ1 − µ2) approaches one as the inputs separate indefinitely.Source: §5.2.1 · pp. 141–142
5.2.2 Stochastic process perspective
- From a stochastic-process perspective, soft inequalities equal CDFs of the difference between independently perturbed inputs, with Gaussian and logistic choices corresponding to Gaussian and Gumbel noise.This probabilistic construction connects the heuristic sigmoid relaxation to location-scale perturbations and CDFs.Source: §5.2 · p. 143, Gaussian case · p. 144, Logistic case · p. 144
Recovering hard inequality operators
- Hard inequalities are recovered by assigning each input a delta distribution, eliminating the perturbation used by soft operators.The construction uses Ui ∼ δµi, where each delta distribution assigns probability one to its corresponding input.Source: Recovering hard inequality operators · p. 144
5.3 Soft equality operators
- Soft equality replaces exact equality with normalized similarity kernels, while probabilistic derivations use the normalized PDF of the input difference and recover kernel forms.A Gaussian-process construction with kernel covariance also recovers the heuristic soft equality operator through correlation.Source: §5.3.1 · pp. 145–146, §5.3.2 · pp. 146–147, §5.3.3 · p. 150
5.5.1 Probabilistic continuous extension
- Continuous logical extensions map [0,1] inputs to [0,1] outputs, using probabilistic product, addition, and complement rules for and, or, not, all, and any.These extensions correspond to probabilities of Boolean events under independent Bernoulli variables.Source: §5.5.1 · pp. 151–152, Figure 5.3 · p. 152
5.5.2 Triangular norms and co-norms
- Triangular norms and conorms provide alternative continuous relaxations of and and or, and associativity generalizes them to all and any over vectors.The probabilistic extensions are examples of t-norms and t-conorms, alongside alternatives summarized in Tables 5.1 and 5.2.Source: §5.5.2 · p. 153, Table 5.1 · p. 154
5.6.1 Differentiating through branch variables
- In if-else statements, branch variables remain differentiable through linear Jacobians, but binary predicates cannot be differentiated directly and yield null, uninformative derivatives when produced by comparisons.As long as both branch functions are differentiable, gradients can pass through their inputs; the predicate remains the problematic variable.Source: §5.6 · p. 154, §5.6.1 · p. 156, §5.6.2 · p. 156
5.6.2 Differentiating through predicate variables
- Predicate derivatives do not influence branch-function derivatives, while branch variables remain differentiable when their constituent functions are differentiable.For fixed predicates, the conditional is linear in the selected branch value and constant in the others.Source: §5.6.2 · p. 156, §5.7.2 · p. 162, §5.7 · p. 163
5.6.3 Continuous relaxations
- Soft comparison operators make predicates differentiable, allowing gradients to backpropagate, but require evaluating both conditional branches.Continuous relaxations extend binary predicates to [0, 1], while the probabilistic view interprets the result as an expectation over the branches.Source: §5.6.3 · p. 157, Probabilistic perspective · p. 157, Figure 5.5 · p. 158
Probabilistic perspective
- Local and global smoothing differ for finite smoothing scales but approximate the original gate and coincide as the smoothing scale becomes sufficiently small.Local smoothing treats comparisons independently, whereas global smoothing takes the expectation of the entire program.Source: Probabilistic perspective · pp. 159–160, Figure 5.6 · p. 160
5.7 Else-if statements
- For K-way conditionals, one-hot vectors encode branch selection, while comparison/logical operators or argmax and argmin can construct the categorical predicate.The conditional remains nondifferentiable with respect to its discrete predicate, despite being valid with respect to branch values when the predicate is fixed.Source: §5.7 · p. 161, Combining booleans · p. 161, Argmax and argmin operators · p. 161
5.7.5 Continuous relaxations
- Argmax and argmin have zero Jacobians almost everywhere and undefined Jacobians at ties, so softargmax provides a continuous alternative with non-null derivatives.The relaxed K-branch result can be interpreted as a categorical expectation, enabling variance computation while evaluating all branches.Source: Figure 5.7 · p. 163, §5.7.5 · p. 164, §5.8 · p. 165
5.8 For loops
- For loops form differentiable computation chains when their iteration functions are differentiable, enabling automatic differentiation through iterations such as unrolled gradient descent.This supports gradient-based optimization of gradient-descent hyperparameters when the objective is twice differentiable.Source: §5.8 · p. 166
5.9 Scan functions
- Scan functions carry an operation’s result across sequence elements, form differentiable computation graphs, and encompass parameterized sequence-to-sequence recurrent neural networks.The abstraction generalizes prefix sums and can exploit studied parallel scan algorithms.Source: §5.9 · pp. 167–168
5.10.1 While loops as cyclic graphs
- While loops are cyclic graphs with an iteration count that may be unknown or infinite, so recursive implementations require lazy evaluation and bounded versions can truncate execution after T iterations.Lazy evaluation is needed to avoid infinite recursion by evaluating the recursive branch only when the predicate requires another iteration.Source: §5.10.1 · p. 169, Importance of lazy evaluation · p. 170, §5.10.2 · p. 170
5.10.2 Unrolled while loops
- A truncated while loop with at most T iterations can be represented as a conditional with T + 1 branches, weighting states s0 through sT by stopping probabilities.The output is cond(˜π, s0, . . . , sT) = ∑_{t=0}^T ˜π_t s_t.Source: §5.10.2 · pp. 171–172, Figure 5.11 · p. 172
5.10.3 Markov chain perspective
- A while loop can be modeled as a Markov chain over its successive states, with a stopping time determining the output and bounded iteration counts recovering the unrolled expression.The chain may be countably infinite when unbounded, while enforcing a maximum of T iterations yields the finite unrolled form.Source: §5.10.3 · pp. 173–174
Smoothed while loops
- Replacing a step stopping criterion with a sigmoid provides informative, generally non-zero derivatives with respect to the tolerance parameter τ.The step-based derivative with respect to τ is always 0, whereas logistic sigmoid smoothing yields non-zero derivatives and an expected smoothed output.Source: Smoothed while loops · p. 175
5.11 Summary
- Differentiating predicate variables requires differentiable relaxation, while a Markov-chain perspective provides a principled way to differentiate through while-loop stopping criteria.Unlike fixed-iteration for loops and scan, while loops may have unknown or infinite iteration counts; unrolling them nevertheless produces valid directed acyclic graphs.Source: §5.11 · pp. 175–176
Data structures
The section develops differentiable representations for lists and dictionaries, showing that data access is differentiable with respect to stored values but not discrete indices or keys. Replacing discrete selectors with probability distributions yields soft operations, including convolutional list access, kernel-regression dictionary access, and Transformer softargmax attention.
6.1.2 Operations on variable-length lists
- Variable-length lists support initialization, end insertion and removal, arbitrary insertion, and stack or queue behavior while preserving differentiability in continuous arguments.Implementations include linked lists with O(K) random access and dynamic arrays with O(1) random access but possible memory reallocations.Source: §6.1.2 · p. 179, Differentiability · p. 180, Inserting values · p. 180
Multi-dimensional indexing
- Multi-dimensional indexing follows the one-dimensional construction through flattening or product distributions, including distributions concentrated on neighboring image pixels.There is a bijection between flat indices and multi-axis indices, so the same soft-indexing framework applies to arrays and tensors.Source: Multi-dimensional indexing · p. 183
6.2 Dictionaries
- A dictionary maps unique keys to values; ordinary get and set operations are differentiable in dictionary contents but not in the queried key.dict.get returns infinity when no key matches, and dict.set leaves the dictionary unchanged when the key is absent.Source: §6.2 · p. 183, Getting values · p. 184, Setting values · p. 185
6.2.2 Continuous relaxation using kernel regression
- Kernel regression relaxes dictionary lookup into a conditional expectation, producing a convex combination of values weighted by key similarity.This is the Bayes predictor under squared-error loss and corresponds to the Nadaraya–Watson estimator; Zhang et al. (2021) previously identified this perspective.Source: Figure 6.2 · p. 185, §6.2.2 · pp. 186–187
6.2.3 Discrete probability distribution perspective
- The discrete-probability perspective maps each query key to an affinity distribution over dictionary keys, enabling soft replacement that is differentiable with respect to that distribution.Unlike dict.set, the corresponding soft operation is differentiable with respect to πk.Source: §6.2.3 · pp. 187–188, Figure 6.3 · p. 188
6.2.4 Link with attention in Transformers
- With Gaussian kernels and unit-normalized keys, the dictionary relaxation recovers softargmax attention from Transformers.The same differentiable dictionary mechanism also underlies content-based memory addressing in neural Turing machines (Graves et al., 2014).Source: §6.2.4 · pp. 188–189, §6.3 · p. 189
6.3 Summary
- Replacing integer indices or dictionary keys with probability distributions produces differentiable relaxations of otherwise nondifferentiable selection operations.Lists remain differentiable with respect to list values, and dictionaries with respect to dictionary contents, but not discrete indices or keys.Source: §6.3 · p. 189
Getting values
- For lists, softGet computes a probability-weighted value that reduces to a discrete convolution, while softSet and softInsert use element-wise expectations over candidate locations.The relaxed outputs may lie in conv(V), and list.softGet has complexity O(M · |supp(πi)|) when V = R^M.Source: Getting values · p. 181, Setting values · p. 182, Inserting values · p. 182
Finite differences · Automatic differentiation
This section covers methods for computing derivatives, from finite-difference approximations to automatic differentiation through compositional programs. It develops forward- and reverse-mode techniques, including their computational, memory, checkpointing, reversibility, and variance trade-offs.
Finite differences
- Finite differences use only function evaluations and can differentiate blackbox functions, but their cost scales with dimension, limiting practical use in high-dimensional machine learning.They are mainly useful for low-dimensional blackbox functions and testing gradient implementations.Source: Finite differences · p. 191, §7.8 · p. 198
7.1 Forward differences
- Forward and backward differences have the same O(δ) truncation error, while finite-difference accuracy must balance truncation error at large δ against round-off error at small δ.Computer evaluation introduces round-off error of order O(ε), producing the trade-off illustrated in Fig. 7.1.Source: §7.1 · p. 192, §7.2 · p. 193, Figure 7.1 · p. 193
7.3 Central differences
- Central differences achieve O(δ2) truncation error, smaller than the O(δ) error of forward and backward differences for δ < 1.The symmetric formula cancels even-order Taylor terms, and the smaller error is confirmed empirically in Fig. 7.1.Source: §7.1 · p. 192, §7.2 · p. 193, Figure 7.1 · p. 193
7.4 Higher-accuracy finite differences
- Higher-order finite differences trade additional evaluations for smaller truncation error, reaching O(δp) in generalized forward schemes and O(δ2p) in generalized central schemes.For higher derivatives, forward approximations achieve O(δp−k+1) with p + 1 evaluations; a second-order central approximation achieves O(δ2) with 3 evaluations.Source: §7.4 · p. 194, §7.5 · pp. 195–196
7.6 Complex-step derivatives
- Complex-step derivatives require one complex-valued function call and usually have smaller round-off error, but all operations in f must support complex numbers.A complex call may cost roughly twice a real call, while avoiding function subtraction reduces round-off error.Source: §7.7 · p. 197, §7.8 · p. 198
7.7. Complexity
- For f : R^P → R, directional finite differences require two function calls, while gradients require P +1 calls for forward or backward differences, 2P for central differences, and P complex calls.These evaluation counts are summarized in Table 7.1.Source: §7.7 · p. 197, Table 7.1 · p. 198
Automatic differentiation
- Automatic differentiation avoids materializing full Jacobians by computing directional products: forward mode propagates JVPs rightward, while reverse mode propagates VJPs leftward through computation chains.Full Jacobians require intermediate Jacobian storage and matrix-matrix multiplications, whereas most applications need only right- or left-multiplication.Source: §8.1 · p. 200, §8.1.1 · p. 201, §8.1.2 · p. 203
8.1.1 Forward-mode
- A JVP follows the original computation while propagating a tangent variable, costing roughly twice the function and exactly twice its memory cost.Forward mode computes each intermediate JVP alongside the corresponding state.Source: §8.1.1 · pp. 201–202
8.1.2 Reverse-mode
- Reverse-mode autodiff performs a forward pass, initializes the output cotangent, and propagates VJPs backward to obtain the input VJP.For a chain, the backward recursion computes rk−1 := ∂fk(sk−1)∗[rk] from k = K down to 1, ending with ∂f(s0)∗[u] = r0.Source: §8.1.2 · p. 204
- Reverse mode stores all intermediate inputs during the forward pass, so memory grows linearly with chain length K and can be reduced using checkpointing.Stored intermediates also allow the backward pass to be reused for n different output directions without recomputing the forward chain.Source: §8.1.2 · p. 205
Complexity of forward-mode autodiff
- Forward mode computes full Jacobians column-by-column through D JVPs, but when D = 1 one directional derivative gives the entire Jacobian.Its intermediate computations can be released when no child node depends on them, yielding memory O(max{D, M}) under the stated dimensional assumptions.Source: Complexity of forward-mode autodiff · pp. 207–208, §8.3.1 · p. 212
When to use forward-mode vs. reverse-mode autodiff?
- When M < D, reverse mode is more advantageous despite memory costs, whereas forward mode is preferable when M ≥ D.For full Jacobians, forward mode requires D JVPs with cost O(MD^2 + KD^3) and memory O(max{D, M}); reverse mode requires M VJPs with cost O(M^2D + KMD^2) and memory O(KD + M).Source: Complexity of forward-mode autodiff · p. 207, When to use forward-mode vs.… · p. 208, Complexity of reverse-mode autodiff · p. 208
8.2 Feedforward networks
- For feedforward networks with layer parameters, reverse-mode differentiation extends the chain recursion to compute parameter gradients together with the input adjoint.With u = ∇ℓ(f(x, w); y), the procedure returns ∇L(w; x, y) = (g1, . . . , gK), while the input gradient is typically unnecessary except in generative modeling.Source: §8.2 · p. 209, §8.2.1 · p. 209, §8.2.2 · p. 210
8.3.1 Forward mode
- In computation graphs, forward mode propagates intermediate values and directions in topological order, summing contributions from multiple parent inputs to produce the final JVP.The graph formulation permits multiple inputs per function and applies the fan-in rule to propagate derivatives.Source: §8.3 · p. 212, §8.3.1 · pp. 212–213
8.3.2 Reverse mode
- Reverse mode on computation graphs handles fan-out by accumulating cotangent contributions from all child nodes into each shared intermediate variable.The duplication operation is a derivational formalism; in practice, shared values are accessed repeatedly and adjoints are updated in place.Source: §8.3.2 · pp. 215–216, Implementation · p. 216
Implementation
- The generic graph implementation performs a topological forward pass that instantiates VJPs, followed by a reverse pass that accumulates parent adjoints in place.Topological ordering ensures each rk is fully computed when its function fk is processed.Source: Implementation · pp. 216–217
8.3.3 Complexity, the Baur-Strassen theorem
- Reverse-mode autodiff gives f and ∇f roughly the same computational complexity, enabling gradient descent on high-dimensional functions; the Baur-Strassen theorem formalizes this for polynomial circuits.For general computation graphs, analogous results are available beyond arithmetic circuits.Source: §8.3.3 · p. 219
8.4.1 Primitive functions
- Autodiff systems compose supported primitive functions, require JVPs and VJPs for forward and reverse modes, and can include higher-level primitives such as numerically stable log-sum-exp.Closure under composition also supports higher-order differentiation and automatic linear transposition; composite functions can be inlined into primitive operations.Source: §8.4.1 · pp. 219–220, §8.4.2 · p. 220, §8.4.3 · p. 220
8.4.3 Examples of JVPs and VJPs
- For linear and separable functions, JVPs and VJPs follow directly from matrix or diagonal Jacobian structure, and separable functions have identical JVP and VJP forms.The examples cover matrix-vector products, matrix-matrix products, and element-wise separable maps.Source: §8.4.3 · p. 220, §8.4 · p. 221
8.4.4 Automatic linear transposition
- The VJP of a VJP is the corresponding JVP, allowing forward-mode autodiff to be implemented from VJPs alone; the reverse construction from JVPs uses partial evaluation.This automatic linear transposition principle exploits linearity in the output direction.Source: §8.4.4 · pp. 221–222
8.5 Checkpointing
- Checkpointing trades computation for memory by storing selected intermediates and recomputing others, spanning full caching with C(K) = K −1 and M(K) = K to constant-memory quadratic computation.Checkpoint locations determine the memory-computation trade-off.Source: §8.5 · pp. 222–224, Figure 8.10 · p. 223
8.5.1 Recursive halving
- For even K, recursive halving roughly halves computational complexity versus Algorithm 8.6 and uses M(K) = log2 K memory, while l = k/2 recovers its complexity in recursive minimization.Its cost satisfies C(K) = 2C(K/2) + K/2. The refined scheme selects l to minimize C∗(k, s) := min 1≤l≤k−1{C∗(k −l, s −1) + C∗(l, s) + l}.Source: §8.5.1 · pp. 224–226, §8.5.2 · p. 227
Analytical formula
- Analytical checkpointing formulas avoid computing dynamic-programming pointers but assume uniform function costs, while online schemes dynamically rewrite S checkpoints when graph length is unknown, such as in while loops.Online methods are proven nearly optimal and apply the fixed-length optimal approach after computation terminates.Source: Analytical formula · p. 228, §8.5.3 · p. 228
8.6. Reversible layers
- Invertible computation chains eliminate reverse-mode storage by reconstructing each predecessor on demand, yielding optimal memory complexity as intermediate states and reverse values are released during backpropagation.This procedure applies when every f_k is invertible and f_k^-1 is readily accessible.Source: §8.6.1 · p. 228, §8.6 · p. 229
- Reversible architectures can be designed to make inverse computation accessible, including reversible residual networks, orthonormal RNNs (Helfrich et al., 2018), neural ODEs, and momentum residual neural networks.Source: §8.6 · p. 229
8.6.2 Case of orthonormal JVPs
- For orthonormal JVPs, reverse-mode autodiff of f coincides with forward-mode autodiff of f^-1.This follows because the VJP of the composed function equals the JVP of its inverse.Source: §8.6.2 · p. 229
8.7 Randomized forward-mode gradient estimator
- An unbiased randomized forward-mode estimator approximates ∇f using JVPs along isotropic Gaussian directions, requiring fewer JVPs than computing all P directions.Averaging M noise-vector contributions avoids storing intermediate activations but increases variance and can slow SGD convergence.Source: §8.7 · p. 230
8.8 Summary
- Forward mode uses JVPs, builds Jacobians by columns, and avoids storing intermediates, whereas reverse mode uses VJPs, builds rows, and stores intermediates; reverse-mode gradient cost is within a constant factor of evaluation cost.Forward mode is efficient for tall Jacobians (M ≥P), reverse mode for wide Jacobians (P ≥M), and checkpointing trades computational efficiency for memory efficiency.Source: §8.8 · p. 231
Second-order automatic differentiation
This section develops second-order automatic differentiation through Hessian-vector products and backpropagation rules, covering mode compositions, computation graphs, and tractable Hessian approximations. It also presents structured alternatives and matrix-free iterative methods for obtaining inverse-Hessian products or approximate solutions without explicitly materializing Hessians.
9.1.1 Four possible methods
- Hessian-vector products can be computed in four ways by combining reverse and forward autodiff: reverse-on-reverse, forward-on-reverse, reverse-on-forward, and forward-on-forward.The HVP is generally preferred to materializing the full Hessian.Source: §9.1 · pp. 232–233, §9.1.1 · p. 232
9.1. Hessian-vector products
- For chain functions, forward-on-reverse is a priori preferable for HVPs: the first three methods have O(KP^2) time and O(KP) space, whereas forward-on-forward has O(KP^3) time.Forward-on-reverse avoids the extra memory cost of reverse-on-reverse and the latter’s branching-related summations; for general functions, the first three methods should be benchmarked.Source: §9.1 · p. 236
9.2.1 An approximation of the Hessian
- For nonconvex compositions, the Hessian is typically indefinite, while the Gauss-Newton matrix is positive semi-definite when the outer loss ℓ is convex and equals the Hessian when f is linear.Thus Gauss-Newton yields a convex quadratic approximation in the direction v under convex ℓ.Source: §9.2.1 · pp. 237–238
9.2.3 Gauss-Newton vector product
- The Gauss-Newton vector product combines the JVP and VJP of f with the HVP of ℓ, requiring 2 forward passes and 1 backward pass through f.The Gauss-Newton chain rule is straightforward, and the full matrix need not be materialized.Source: §9.2.2 · p. 238, §9.2.3 · p. 238, §9.2 · p. 239
9.2.4 Gauss-Newton matrix factorization
- When ℓ is convex, the Gauss-Newton matrix factors as V V^⊤; forming it costs O(M^3) time plus M VJP calls, uses P × M memory, and trades computation against direct evaluation.The factorization enables GNVPs as V V^⊤v.Source: §9.2.4 · p. 239
9.3.3 Equivalence with the Gauss-Newton matrix
- The Fisher information is the covariance of score gradients and, under twice differentiability with integrable second derivatives, can also be expressed as the Hessian of the negative log-likelihood.In supervised learning with an exponential-family pθ, the Fisher information matrix is equivalent to the Gauss-Newton matrix.Source: §9.3.1 · p. 240, §9.3.2 · p. 241, §9.3.3 · pp. 241–242
9.4.2 Implementation with matrix-free linear solvers
- An inverse-Hessian vector product solves H v = u using only Hessian-vector products, enabling matrix-free solvers such as CG for positive-definite H and GMRES for any invertible H.For strictly convex functions, CG can compute the IHVP; in theory it reaches the exact solution within P iterations, though numerical errors may intervene.Source: §9.4.1 · p. 243, §9.4.2 · pp. 243–244
9.4.3 Complexity
- Matrix-free CG or GMRES costs O(P^2) per iteration and can stop at approximate solutions, whereas decomposition methods cost O(P^3) to build.Well-conditioned matrices need few iterations; poorly conditioned matrices can cause instability and require more than P iterations, while LU and QR are generally more stable when H entries are accessible.Source: §9.4.3 · p. 244, §9.5 · p. 245
9.5. Second-order backpropagation
- With only Hessian-vector products, recovering all Hessian entries requires P HVP calls, so CG or GMRES can approximate solutions sooner without storing the matrix.This makes matrix-free methods preferable when accessing Hessian entries is prohibitively expensive.Source: §9.5 · p. 245
9.5.1 Second-order Jacobian chain rule
- The Hessian chain rule contains intermediate-function curvature and a term showing how first-order variations affect second-order derivatives.The rule generalizes Hessians to vector-valued functions and supports second directional derivatives of composed twice differentiable functions.Source: §9.5.1 · p. 246
9.5.2 Computation chains
- For computation chains, Hessian backpropagation propagates both vectors and matrices using intermediate first- and second-order derivatives.For fan-in graphs, cross-product terms couple branches, whereas fan-out can treat each output function separately.Source: §9.5.2 · p. 247, §9.5.3 · p. 248
9.6 Block diagonal approximations
- Block-diagonal Hessian and Gauss–Newton approximations are easier to invert, and feedforward-network blocks can be computed in one backward pass sharing most computations.The resulting outputs provide block-wise Gauss–Newton approximations for parameter groups.Source: §9.6 · p. 249, §9.6.1 · p. 250
9.6.2 Computation graphs
- For generic computation graphs, the blockwise recursion omits cross-derivatives and therefore gives an approximation rather than an exact Hessian, while remaining easy to implement.The same framework extends blockwise Gauss–Newton computation to graph-structured models.Source: §9.6.2 · p. 251, §9.7 · p. 251
9.7 Diagonal approximations
- Diagonal Hessian backpropagation replaces full Hessians with vector surrogates, runs in linear time per output dimension, and extends from feedforward networks to computation graphs.The resulting algorithms return gradients alongside approximate Hessian diagonals, including parameter-wise diagonals for multilayer perceptrons.Source: §9.7 · pp. 252–253, §9.7.1 · p. 252, §9.7.2 · pp. 253–254
9.8 Randomized estimators
- Randomized diagonal estimators are unbiased Monte Carlo methods whose variance decreases with more samples and matvecs, while Bartlett factorizations approximate Gauss–Newton matrices using sampled gradients.The diagonal Bartlett estimator uses element-wise products of gradients evaluated at sampled labels and is exact when all labels in the support are sampled.Source: §9.8.1 · p. 254, §9.8 · p. 255, §9.8.2 · p. 255
1 S diag(∇2
Hessian-vector products provide an efficient alternative to full Hessian backpropagation, while Gauss–Newton and other approximations extend scalable second-order computation. A gradient-of-averaged-loss estimator is simpler to implement but may have higher variance.
1 S diag(∇2
- The gradient-of-averaged-loss estimator is more convenient to implement because it requires only gradients of averaged losses, but it may have higher variance.It becomes exact when all labels in the support of pθ have been sampled; Liu et al. (2023) use the one-y-per-x special case.Source: §1 · p. 256
9.9 Summary
- Hessian backpropagation is simple for computation chains but computationally prohibitive for graphs because fan-in creates cross-product terms.Source: §9.9 · p. 257
- Hessian-vector products offer four computation methods, with forward-on-reverse a priori the most efficient.Computing HVPs costs only a constant factor more than evaluating the function itself.Source: §9.9 · p. 257
- The Gauss–Newton matrix approximates the Hessian, is positive semidefinite for convex ℓ, and supports efficient Gauss–Newton-vector products.For exponential families, it is equivalent to the Fisher information matrix and yields a principled quadratic approximation.Source: §9.9 · p. 257
- Other scalable second-order approximations include block-diagonal, diagonal, and randomized estimators.Source: §9.9 · p. 257
Inference in graphical models as differentiation
This section develops graphical-model inference through its connections to differentiation, covering dependency representations, distribution factorizations, and key inference tasks. It further unifies these ideas via the log-partition function, semirings, dynamic programming, and backpropagation.
Inference in graphical models as differentiation
- Graphical models encode dependencies among random variables, thereby determining how their joint probability distributions factorize and connecting inference to computation chains and graphs.The chapter uses this connection to relate graphical-model concepts to differentiation.Source: Inference in graphical models as… · p. 258, §10.3.1 · p. 260
10.1 Chain rule of probability
- The chain rule factorizes any joint distribution into conditional probabilities without independence assumptions, while conditional independence can simplify this factorization.Conditional independence means that a variable’s distribution given one set is unchanged by another set once a conditioning set is known.Source: §10.2 · pp. 259–260
10.3. Inference problems
- The described operations include computing likelihoods and log-likelihoods, marginalizing variables by summing or integrating over possible values, and computing expectations whose mean lies in the convex hull of feature values.Source: §10.3.2 · p. 260, §10.3 · p. 261, §10.3.4 · p. 261
10.3.6 Complexity of brute force
- Brute-force marginal, mode, and expectation computation requires O(∏_{k=1}^K |S_k|) time, or O(M^K) when every variable has M values.Likelihood computation is described as trivial by comparison.Source: §10.3.6 · p. 263
10.4.1 The Markov property
- A Markov chain imposes that each S_k depends only on S_{k−1}, paralleling a computation chain whose functions map each preceding state to the next.Time-homogeneous chains reuse the same conditional distribution across k, while nth-order chains condition on the previous n variables.Source: §10.4.1 · p. 264, §10.4.2 · p. 266, §10.4.3 · p. 267
10.4.3 Higher-order Markov chains
- Higher-order Markov structure describes autoregressive models such as Transformers with context window n, but larger contexts make exact dynamic-programming inference intractable, motivating beam search or ancestral sampling.This connects model context size directly to inference tractability.Source: §10.4.3 · p. 267
10.5.1 Expressing variable dependencies using DAGs
- Bayesian networks generalize Markov chains by representing dependencies with a DAG and factorizing the joint distribution into conditional distributions over parent variables, a representation suited to causal relationships.Unlike computation graphs, Bayesian networks are specified by conditional distributions rather than deterministic functions; neural networks may parameterize root-node or autoregressive distributions.Source: §10.5.1 · pp. 267–268, §10.5.2 · p. 268
10.5.3 Ancestral sampling
- Because normalized Bayesian-network conditionals yield an automatically normalized joint distribution, ancestral sampling draws independent samples by traversing nodes in topological order.At each node, the sample is drawn conditionally on previously sampled parent values.Source: §10.5.3 · p. 269
10.7.1 The forward-backward algorithm
- Forward-backward computes chain marginals through independent forward and backward recursions, with log-domain implementations improving numerical stability and permitting parallel execution.The method uses distributivity to form recursive sums and can compute normalization, marginal, and conditional probabilities.Source: §10.7.1 · pp. 273–274
10.8 Inference on trees
- Efficient dynamic-programming inference extends from chains to trees and polytrees, where sum-product performs marginal inference and max-sum performs MAP inference.Sum-product is also called belief propagation or message passing because it propagates local messages through the graph.Source: §10.8 · p. 278
10.9.1 Inference as gradient of the log-partition
- The gradient of the log-partition function yields unary and pairwise marginal probabilities, and reverse-mode autodiff computes them with roughly the complexity of evaluating the log-partition function.With binary sufficient-statistic encodings, [∇A(θ)] directly indexes unary and pairwise marginals.Source: §10.9.1 · p. 278, §10.9 · p. 279
10.9. Inference as differentiation
- In the zero-temperature limit, the gradient-derived quantity becomes a binary encoding of the maximum a-posteriori solution.This connects marginal inference and MAP inference through the temperature-scaled log-partition function.Source: §10.9 · pp. 278–279
10.9.2 Semirings and softmax operators
- Semirings unify forward-backward and Viterbi recursions by substituting operation pairs such as (+, ×), (max, +), and (maxε, +).Associativity and distributivity enable decomposition into smaller computations, while sparsemax is not associative.Source: §10.9.2 · p. 280
10.9.3 Inference as backpropagation
- Backtracking in dynamic-programming MAP inference is recovered algorithmically as a special case of backpropagation.The construction is presented for chain states with S0 = {1} and Sk = {1, . . . , M}.Source: §10.9.3 · p. 281
Q := argmaxε … Differentiating through integration
This section develops methods for differentiating through optimization-defined, equation-defined, and integral computations, including implicit functions, nonlinear and inverse problems, and expectations. It connects soft backpropagation, adjoint and implicit differentiation, envelope theorems, reverse-mode autodiff, and Monte Carlo estimation, with attention to conditions, sensitivity, and parameter-dependent sampling.
Q := argmaxε
- Soft backpropagation replaces hard backpointers with probability distributions, and its zero-temperature limit yields a binary encoding of the backtracking solution.Inference on chains and trees is equivalent to differentiating log-partition functions; soft max corresponds to sum-product, while hard max corresponds to max-plus.Source: Q := argmaxε · p. 282, §10.10 · pp. 283–284
Differentiating through optimization
- Differentiable programming studies derivatives through optimization problems and nonlinear systems, including implicit functions that lack explicit elementary decompositions.Implicit functions can represent physical equilibria, optimization solutions, or roots of nonlinear systems, enabling parameter sensitivity analysis.Source: Differentiating through optimization · p. 285, §11.1 · p. 285, §11.1.1 · p. 286
11.1.3 Application to bilevel optimization
- Bilevel optimization differentiates an outer objective through an inner optimization solution, as required for hyperparameter optimization.The inner objective trains model parameters, the outer validation loss selects hyperparameters, and gradient descent requires ∇h(λ).Source: §11.1.3 · p. 286, Figure 11.1 · p. 287
11.2 Envelope theorems
- Envelope theorems remove the need to differentiate through an argmax when the value function is formed by maximizing the same function.Under suitable assumptions, the optimizer can be treated as constant with respect to λ, although the assumptions differ across theorem variants.Source: §11.2 · pp. 287, 289
11.2.1 Bertsekas’ theorem
- Bertsekas’ theorem assumes concavity and a unique maximizer, whereas Rockafellar’s theorem drops concavity but requires stronger differentiability assumptions.The section distinguishes these results from Danskin’s original theorem, which gives a directional derivative as the maximum directional gradient over maximizers.Source: §11.2.1 · p. 288, §11.2.2 · p. 290, §11.2.3 · p. 291
11.3.1 Univariate functions
- The implicit function theorem locally differentiates w⋆(λ) when F(w, λ) = 0 has a continuously differentiable map and an invertible partial derivative in w.In one dimension the condition is ∂1F ≠ 0; in the multivariate case, ∂1F must be invertible, but the theorem does not construct w⋆(λ).Source: §11.3.1 · pp. 291, 293, §11.3.2 · p. 294
11.3.2 Multivariate functions
- For unconstrained optimization, stationarity turns the optimizer into a root of F, while constrained optimization uses a projection fixed point before applying the implicit function theorem.The unconstrained case requires an invertible Hessian; constrained differentiation additionally requires the projection Jacobian, which is easy to compute for many sets.Source: §11.3.2 · pp. 294–295
11.3.3 JVP and VJP of implicit functions
- JVPs and VJPs integrate implicit functions into autodiff by reducing derivative computation to linear-system solves under the implicit function theorem’s assumptions.This provides the operator needed to differentiate implicitly defined solutions without requiring an explicit decomposition.Source: §11.3.3 · pp. 295–296
11.4.1 Differentiating nonlinear equations
- The adjoint state method computes gradients of explicit objectives composed with uniquely defined implicit states under equality constraints.It introduces an adjoint variable obtained from a linear system, enabling gradient-based optimization of the resulting unconstrained objective.Source: §11.4.1 · pp. 297–298
11.4.2 Relation with envelope theorems
- The adjoint state method is an envelope-theorem formulation when parameters enter both the objective and the equality constraint.Its derivation uses Lagrange multipliers or the implicit function theorem, with the multiplier identified as the adjoint state.Source: §11.4.2 · p. 298, §11.4.3 · p. 299, §11.4.4 · p. 299
11.4.5 Reverse mode as adjoint method with backsubstitution
- Reverse-mode autodiff is recovered as an adjoint-state method for feedforward networks, where the resulting upper-triangular system is solved efficiently by backsubstitution.The triangular structure follows because each activation depends only on preceding activations, and the formulation also supports VJPs of vector-valued functions.Source: §11.4.5 · p. 300, Feedforward networks as the solution… · p. 300, Solving the linear system using… · pp. 301–302
11.5.1 Differentiating inverse functions
- The inverse function theorem states that an invertible Jacobian yields a locally bijective, continuously differentiable inverse whose Jacobian is the inverse of the original Jacobian.The implicit function theorem recovers this result by choosing F(w, ω) = f(w) − ω, making inverse differentiation a special case of implicit differentiation.Source: §11.5.1 · p. 303, §11.5.2 · p. 303, §11.5.3 · p. 305
11.3.3 JVP and VJP of implicit functions
- The implicit function theorem provides existence and derivative formulas for implicitly defined solutions, with JVPs and VJPs computed through associated linear systems.The systems can be solved matrix-free using conjugate gradient, GMRES, or BiCGSTAB when appropriate.Source: §11.3.3 · p. 296, §11.3.4 · pp. 296–297, §11.6 · p. 305
12.1 Differentiation under the integral sign
- Differentiation under the integral sign is valid under sufficient measurability, integrability, absolute-continuity, and local-integrability conditions.Absolute continuity can weaken differentiability assumptions and can smooth functions differentiable almost everywhere, such as ReLU.Source: §12.1 · pp. 306–307
12.2.1 Parameter-independent distributions
- For parameter-independent distributions, differentiation commutes with expectation under the same mild conditions, so gradients can be estimated by differentiating Monte Carlo sample averages.For N i.i.d. samples, both the value and gradient estimators are unbiased and converge to the true quantities as N → +∞.Source: §12.2 · p. 307, §12.2.1 · pp. 307–308
- Higher-order derivatives are approximated by computing ∇2 bF_N(θ) with autodiff.The implementation applies the same principle used for first derivatives.Source: §12.2.1 · p. 308
- This parameter-independent approach implements an approximate objective and differentiates it with autodiff, analogous to discretize-then-optimize for ODE solutions.Source: §12.2.1 · p. 308
12.2.2 Parameter-dependent distributions
- Differentiating E(θ) = E_Y∼pθ[g(Y)] is difficult because θ parameterizes the sampling distribution, while g may be a black box without gradients.Monte Carlo estimation requires sampling from pθ, but it is unclear how to differentiate with respect to θ when θ is involved in sampling.Source: §12.2.2 · pp. 308–309
Continuous case … Pushforward measures
This section develops methods for differentiating expectations involving sampled distributions, including score function, path gradient, and reparameterization estimators, with attention to variance and autodiff-compatible implementation. It formalizes differentiable transformations through pushforward measures and change-of-variables results, and connects stochastic programs with probability distributions and stochastic computation graphs.
Continuous case
- Exact gradients are generally impractical: continuous-case integrals may be intractable, while discrete-case summation is tractable only for small or sparse supports and inexpensive g.For discrete distributions, exact gradients are therefore seldom used in practice.Source: Continuous case · p. 309, Discrete case · pp. 309–310
12.2.3 Application to expected loss functions
- Expected-loss optimization reduces risk gradients to gradients of expectations, with score-function estimators supporting discrete losses and pathwise estimators requiring differentiable losses.In experimental design, stochastic optimization can lift costly evaluations g(λ) to Eλ∼pθ[g(λ)]; SFE avoids requiring gradients of g, and normal distributions provide a natural choice for continuous hyperparameters.Source: §12.2.3 · pp. 310–311, §12.2.4 · p. 311
12.3.1 Scalar-valued functions
- SFE rewrites the gradient of an expectation as an expectation, enabling unbiased Monte Carlo estimation when sampling and the score function are available in closed form.SFE does not require differentiability of g, making it suitable for discrete losses and blackbox functions.Source: §12.3.1 · pp. 312–313
12.3.2 Variance reduction
- SFE’s high variance can require many samples, motivating constant baselines, unbiased REINFORCE leave-one-out estimators, and control variates based on tractable approximations.A running-average baseline is common but is neither optimal nor guaranteed to reduce variance; RLOO remains unbiased through sample independence.Source: Baseline · pp. 315–316, REINFORCE leave-one-out (RLOO) · p. 316, Control variates · p. 317
Differentiating through both the distribution and the function
- SFE extends to parameter-dependent functions and joint distributions through the product rule, while vector-valued outputs use VJPs, JVPs, and outer products.The joint-distribution extension follows by applying Proposition 12.1 to ρθ := pθ·qθ, and generalizes straightforwardly to more variables.Source: Differentiating through both the… · pp. 314, 318, Differentiating through joint… · p. 315, §12.3.3 · pp. 317–318
12.3.4 Second derivatives
- SFE supports first- and second-order derivatives, with its score-function moments connected to Bartlett’s first and second identities.The Hessian estimator follows by applying the parameter-dependent-function result to g(y)∇θ log pθ(θ).Source: §12.3.4 · p. 318, Link with the Bartlett identities · p. 319
12.3.5 Autodiff-friendly implementation
- The magicbox operator makes score-function estimators compatible with autodiff by using stop-gradient behavior that preserves evaluation while blocking derivative flow.Monte Carlo samples can then be drawn without differentiating through sampling while autodiff recovers first- and second-derivative expressions.Source: §12.3.5 · p. 319
12.4 Path gradient estimators, reparametrization trick
- Path gradient estimators move parameters from the sampling distribution into the function via a change of variables, enabling chain-rule differentiation; location-scale transforms instantiate this for normal and related distributions.For U = µ + σZ with Z ∼Normal(0, 1), the parameters µ and σ are differentiated through the transformed function.Source: §12.4 · p. 321, §12.4.1 · pp. 321–322
12.4.1 Location-scale transforms
- In location-scale transforms, the transformed variable’s density can be obtained by differentiating its cumulative distribution function.This connects distributional transformations with density evaluation.Source: §12.4.1 · p. 322
12.4.2 Differentiable transforms
- Differentiable transformations generalize the path gradient estimator, yielding an unbiased gradient estimator when a parameter-free source variable is transformed into the target distribution.The estimator requires differentiability of the loss and that the loss be defined on the full Euclidean space, which can fail for zero-one or ranking losses.Source: §12.4.2 · pp. 322–323
12.4.3 Inverse transforms
- Inverse transform sampling generates a target distribution by applying its quantile function to a uniform random variable, and differentiable quantiles can serve as reparametrization transforms.This approach requires access to the quantile function and its derivatives, which may be difficult for complicated distributions.Source: §12.4.3 · pp. 324–325
Pushforward distributions
- Pushforward distributions formalize reparametrization by mapping a source distribution through a continuous, potentially parameterized transformation.They can be sampled easily from the source distribution and support expectation computation without an explicit target density.Source: Pushforward distributions · pp. 326–327
Pushforward measures
- The pushforward operator maps measures through continuous functions while preserving positivity and total mass, so probability measures remain probability measures.This provides the measure-theoretic foundation for pushforward distributions.Source: Pushforward measures · p. 327
12.4.5 Change-of-variables theorem
- For an invertible differentiable transformation, the change-of-variables theorem gives an explicit density for the pushforward distribution using the inverse map and its Jacobian.Normalizing flows are designed so the inverse transformation and inverse-Jacobian computation are tractable.Source: §12.4.5 · p. 328
12.5 Stochastic programs
- A stochastic program induces probability distributions over both its final outputs and its execution trajectories because its operations involve randomness.Intermediate variables can also be random variables.Source: §12.5 · p. 329
12.5.1 Stochastic computation graphs
- Stochastic computation graphs represent stochastic programs using conditional sampling and function-evaluation operations that may produce deterministic or random variables.This exposition extends the original stochastic-computation-graph framing by explicitly supporting both intermediate operation types.Source: §12.5.1 · p. 329
Function and distribution nodes … Infimal convolution with a regularization function
This section develops core methods for differentiable programming, covering stochastic computation graphs, gradient estimation, and neural ODE discretization and adjoints. It also addresses smoothing non-differentiable functions through infimal convolution and Moreau envelopes.
Function and distribution nodes
- Stochastic computation graphs combine function and conditional distribution nodes in a directed acyclic graph, with deterministic input s0 and random output SK.Distribution nodes provide randomness, while edges encode dependencies through parent variables.Source: Function and distribution nodes · p. 329, Executing a stochastic program · p. 331
Deterministic and random variables
- Intermediate variables and node inputs are split into deterministic variables si and random variables Si, allowing functions and distributions to condition on both.Parents are partitioned as pa(k) = determ(k) ∪ random(k).Source: Deterministic and random variables · p. 330
Executing a stochastic program
- Executing nodes in topological order produces deterministic or random outputs according to their functions or conditional distributions, and one execution samples the program’s induced output distribution.A distribution node may have no parents, and function-node outputs become random whenever at least one parent is random.Source: Executing a stochastic program · pp. 330–331
Special cases / 12.5.2 Examples
- All-function graphs recover computation graphs, whereas all-distribution graphs recover Bayesian networks; examples show that stochastic gradient expressions can become sufficiently complex to motivate automatic differentiation.The examples cover SFE, pathwise, and SFE-plus-chain-rule estimators.Source: Special cases · p. 332, §12.5.2 · pp. 332–333
12.5.3 Unbiased gradient estimators
- Expected outputs convert stochastic-program differentiation into deterministic optimization: pathwise estimators apply to function-node graphs, while SFE estimators apply to distribution-node graphs and yield unbiased expectation gradients.Executing once gives an i.i.d. sample from p(·|s0), and later work extends SFE to high-order differentiation, node-specific estimators, and control variates.Source: §12.5.3 · p. 334
Converting distribution nodes into function nodes and vice-versa
- Distribution nodes can be converted to function nodes through reparametrization, while function nodes can become distribution nodes through pushforward distributions and change of variables; pathwise estimation is preferred when available because it has lower variance.The conversions are illustrated in Fig. 12.1.Source: Converting distribution nodes into… · p. 335, Figure 12.1 · p. 335
12.5.4 Local vs. global expectations
- Global smoothing replaces the final random output by its expectation, preserving execution distributions and avoiding convex-hull domain requirements, whereas local smoothing replaces intermediate variables by expectations but may require evaluating all branches.Both strategies support gradient estimation through reparametrization or SFE according to node type.Source: §12.5.4 · p. 336
From residual networks to neural ODEs / Existence of a solution / Integration methods / 12.6.2 Continuous adjoint method
- Neural ODEs model continuous-time state changes, and under mild conditions their solutions are unique; numerical integration shows residual networks are discretized non-autonomous neural ODEs via explicit Euler.Continuous adjoints compute gradients by solving a backward ODE, while optimize-then-discretize integrates both state and adjoint trajectories before discretization.Source: From residual networks to neural ODEs · p. 337, Existence of a solution · p. 337, Integration methods · p. 338
12.6.4 Gradients via reverse-mode on discretization
- Discretize-then-optimize generally outperformed optimize-then-discretize, although the former can require substantial memory while the latter incurs numerical errors from mismatched discretizations.Checkpointing can amortize reverse-mode memory, while reversible schemes offer a middle ground; discretize-then-optimize computes exact gradients of a discretized surrogate objective, not necessarily the true ODE.Source: §12.6.4 · p. 344, §12.7 · p. 350
Leapfrog method
- The symmetric reversible leapfrog method exactly recovers the original input after backward discretization, so repeated forward-backward passes produce the same gradient.Its additional context variables can be interpreted as velocities, motivating second-order dynamics and momentum neural networks with low-memory reversible differentiation.Source: Leapfrog method · p. 346, Momentum residual networks · pp. 346–347
12.6.6 Proof of the continuous adjoint method
- The continuous adjoint method introduces a continuously differentiable adjoint function satisfying a backward ODE, with integration by parts yielding vector-Jacobian products for parameters and inputs.The derivation assumes the relevant partial derivatives exist, and existence of the backward adjoint curve follows from Picard–Lindelöf’s theorem under the stated assumptions.Source: §12.6.6 · pp. 347–348
12.7 Summary
- For expectation gradients, score-function estimators are unbiased but high-variance, whereas path-gradient estimators are unbiased and low-variance when integration and differentiation can be interchanged.Score-function estimation suits explicitly available log-PDFs; path gradients suit pushforward distributions generated by injecting randomness through transformations.Source: §12.7 · p. 349
12.6.5 Reversible discretization schemes
- Reversible discretization schemes enable recomputation of intermediate states during backward passes, reducing memory costs and preserving consistency between repeated forward and backward computations.A reversible method provides a closed-form inverse for each discretization step; symmetric schemes make the backward scheme exactly the forward inverse.Source: §12.6.5 · p. 345, Leapfrog method · p. 346, §12.7 · p. 350
Smoothing by optimization
- Smoothing replaces non-differentiable functions with differentiable approximations using infimal convolution and its Legendre-Fenchel dual, including smooth ReLUs, max operators, step functions, and argmax relaxations.The primal approach operates in the function’s original space and focuses particularly on the Moreau envelope.Source: Smoothing by optimization · p. 352, §13.1 · p. 352
13.1.1 Infimal convolution
- Infimal convolution combines functions by minimizing their summed values under a decomposition constraint, analogously to ordinary convolution with integration replaced by minimization.Its infimum is finite under suitable conditions, including convexity of the summed integrand for every input, though convexity is not necessary when domains are compact.Source: §13.1.1 · p. 353, Existence · pp. 353–354
13.1.2 Moreau envelope
- Using the squared 2-norm as regularization yields the Moreau envelope, which minimizes a trade-off between proximity to the input and the original function value.When the minimizer exists and is unique, it defines the associated proximal operator.Source: Infimal convolution with a… · p. 354, §13.1.2 · pp. 354–355
Properties … The softplus
This section develops smoothing and dual-representation tools, focusing on Moreau envelopes, convex conjugates, and dual regularization. It examines their properties, extensions, computational and duality aspects, and generalized entropies such as softplus smoothing of ReLU.
Properties
- For convex f, the Moreau envelope is everywhere differentiable with Lipschitz-continuous gradients, remains convex, has a proximal-operator gradient expression, and shares f’s infimum.It is specifically 1-smooth when f is convex, and it also satisfies a Moreau decomposition involving the convex conjugate.Source: Properties · pp. 355–356
Examples
- The Moreau envelope smooths the 1-norm into the Huber loss, while numerical infimal-convolution approximations also apply to nonconvex functions.The nonconvex examples include a continuous ramp function and a discontinuous step function evaluated over a finite grid and closed interval.Source: Examples · pp. 356–357, Figure 13.2 · p. 358
13.1.3 Vector-valued functions
- Coordinatewise smoothing of vector-valued functions requires T separate minimizations, whereas directional VJP smoothing can use one optimization problem independent of T.The directional approach can support reverse-mode autodiff in principle, but the optimization may become nonconvex and its smoothing strength depends on the VJP direction; this remains an open research area.Source: §13.1.3 · pp. 358–360
13.2 Legendre–Fenchel transforms, convex conjugates
- The Legendre–Fenchel transform maps f to a convex conjugate f∗ that represents f through tight affine lower bounds and tangent slopes rather than its original graph.The intercept of the tightest affine lower bound with slope v is −f∗(v), and f∗ is convex even when f is not.Source: §13.2 · p. 360, §13.2.1 · pp. 360–361
13.2.2 Closed-form examples
- Closed-form conjugates exist for some functions, while constraints can be encoded with indicator functions, yielding support functions and linear maximization oracles.The section gives quadratic and entropy-related examples and notes that conjugate computation can otherwise require difficult maximization.Source: §13.2.2 · pp. 361–362, Constraining the domain · p. 363
13.2.3 Properties
- Convex conjugates are always convex, are differentiable when the defining supremum has a unique maximizer, and satisfy f = f∗∗ exactly for convex closed f.Without unique attainment, the conjugate is instead sub-differentiable; otherwise f∗∗≤f.Source: §13.2.3 · pp. 363–364
13.2.4 Conjugate calculus
- Conjugate calculus provides rules for scalar multiplication, affine shifts and translations, invertible linear compositions, and sums represented through infimal convolution.These rules derive new conjugates from existing ones when direct derivation is difficult.Source: §13.2.4 · p. 365
13.2.5 Fast Legendre transform
- For convex f, projected gradient descent or Frank–Wolfe can compute f∗ to arbitrary precision, whereas grid-based computation without convexity costs O(|U| · |V|) and is limited to low dimensions.With N grid points per dimension, typical grids have |U| = |V| = O(N^M), making the approach practical mainly for M ∈{1, 2, 3}.Source: §13.2.5 · p. 366
13.3.2 Smoothing by dual regularization
- Strongly convex dual regularization makes f_Ω smooth and convex, even when f is nonconvex, and yields almost-everywhere differentiable gradients.If Ω is μ-strongly convex, f_Ω is 1/μ-smooth; this property supports differentiable almost-everywhere relaxations of step and argmax operators.Source: §13.3.2 · pp. 367, 370
Dual approach
- Dual smoothing bounds the smoothed function using lower and upper bounds on the regularizer, while strongly convex regularization controls its smoothness.The construction can be applied to any closed convex f, with approximation error quantified under L_Ω ≤ Ω(v) ≤ U_Ω.Source: §13.3.2 · p. 368, Dual approach · p. 371
13.3.3 Equivalence between primal and dual regularizations
- For convex closed f and R, dual and primal regularization coincide: f_Ω = (f* + Ω)* = f □ R when Ω = R*.The equivalence makes the choice between approaches primarily a matter of mathematical or algorithmic convenience.Source: §13.3.3 · pp. 370–371
Primal approach
- For positively homogeneous functions such as norms, dual regularization recovers corresponding scaling relationships and connects to primal regularization.The supplied passages identify this homogeneous case as a special setting for the dual formulation.Source: Dual approach · p. 371, Primal approach · pp. 371–372
13.3.5 Generalized entropies
- Entropy is expected information content, and negative entropy serves as a natural dual regularizer for discrete distributions and smoothed max operators.Generalized entropies need not have the expected-information form, but entropy functions are nonnegative and uniquely maximized by the uniform distribution.Source: §13.3.5 · p. 372, Information content and entropy · p. 372, Definition and properties of… · pp. 375–376
Shannon’s entropy
- Shannon negentropy can smooth f when dom(f*) ⊆ △M, whereas Gini negentropy applies when dom(f*) ⊆ RM, covering a larger function class.Shannon’s negentropy is strongly convex with respect to ∥·∥1 on △M, while Gini’s is strongly convex with respect to ∥·∥2 on RM.Source: Shannon’s entropy · p. 373, Gini’s entropy · p. 374
Tsallis entropies
- Tsallis entropy interpolates between Shannon entropy as α →1 and Gini entropy at α = 2, while squared p-norm entropy is strongly concave for p ∈(1, 2] and usable for smoothing.These entropy families provide alternative strongly concave regularizers for smoothing functions.Source: Tsallis entropies · p. 374, Definition and properties of… · p. 376
The softplus
- Applying Shannon negentropy to the ReLU conjugate yields the softplus smoothing.The ReLU conjugate is the indicator of [0,1], and the regularizer Ω(π) = π log π + (1 −π) log(1 −π) produces the stated result.Source: §13.4 · p. 377, The softplus · p. 377
The sparseplus … Computation using a convolution
This section develops differentiable relaxations of sparse and nonsmooth operators through strongly convex regularization, gradients, and convolution-based smoothing. It also examines convolution from averaging, perturbation, computational, and transform-based perspectives, including its connection to soft infimal convolution.
13.5.1 Definition and properties
- Smoothing max operators yields smooth counterparts with bounds, permutation commutativity, and translation distributivity, enabling differentiable dynamic programs and smoothed max-pooling layers.The smoothed max is useful in neural networks, while its properties support differentiable dynamic programs.Source: §13.5.1 · pp. 378–379
13.5.2 Reduction to root finding
- For strongly convex regularization, computing maxΩ reduces to solving a root equation involving δΩ, derived using simplex duality, Danskin’s theorem, and differentiation with respect to the scalar shift.The reduction keeps non-negativity explicit, dualizes the simplex equality constraint, and obtains the root equation from the derivative of τ + δΩ(u − τ1).Source: §13.5.2 · pp. 379–380
13.5.3 The softmax
- The entropy regularizer produces softmax as logsumexp, whose distinctive associativity property holds for all scalar inputs.The softmax is obtained analytically from the root-finding formulation and satisfies softmax(softmax(a, b), c) = softmax(a, softmax(b, c)).Source: §13.5.3 · p. 381
13.5.4 The sparsemax
- Gini’s negentropy produces sparsemax, which can be computed using Euclidean projection onto the probability simplex; unlike softmax, it yields sparse outputs and is differentiable only almost everywhere.The exact threshold can be obtained from sorted inputs, while approximate computation can use bisection or gradient descent.Source: §13.5.4 · pp. 381–383, Figure 13.7 · p. 383
13.5.5 Recovering smoothed ReLU functions
- Using u = (u, 0), smoothed max operators recover smoothed ReLU functions, with Shannon and Gini regularizers yielding their corresponding entropy-based and sparse forms.For Gini’s negentropy, the induced regularizer is Ψ(π) = π(π −1), the form used to smooth ReLU.Source: §13.5.5 · p. 384
13.6 Relaxed step functions (sigmoids)
- Regularized step functions are Lipschitz continuous and therefore differentiable almost everywhere under strongly convex regularization; the logistic sigmoid is everywhere differentiable, whereas the sparse sigmoid reaches 0 or 1 but has two kinks.The sparse sigmoid is piecewise linear, while the logistic function arises from Shannon’s negentropy.Source: §13.6 · p. 384, The logistic function · p. 384, The sparse sigmoid · p. 385
13.7 Relaxed argmax operators
- Each relaxed sigmoid and argmax mapping is the gradient of its corresponding smoothed ReLU or max operator, by Danskin’s theorem; softargmax is everywhere differentiable, while sparseargmax is sparse and differentiable almost everywhere.Sparseargmax is Euclidean projection onto the probability simplex, computable in O(M) expected time and O(M log M) worst-case time.Source: Link between smoothed ReLU functions… · p. 385, §13.7 · p. 386, The softargmax · p. 386
Smoothing by integration
- Convolution smooths functions by integrating products of reflected and shifted inputs; convolution with a Gaussian kernel corresponds to smoothing with normal noise of variance σ2.Convolution is commutative, and the Gaussian kernel is the density of a zero-mean normal distribution.Source: §14.1.1 · pp. 391–392, §14.1.2 · p. 392
Averaging perspective
- Convolution with a normalized kernel averages function values, weighting points near µ more heavily while σ controls the influence of distant points.For a Gaussian kernel, the convolution is the expectation of f(U) for a Gaussian centered at µ; the averaging interpretation extends to translation-invariant location-scale kernels.Source: §14.1.2 · p. 392, Averaging perspective · p. 393
Perturbation perspective
- Gaussian convolution is equivalently Gaussian input perturbation, converges to the original function as σ →0, and can be implemented discretely by sampling and renormalizing the kernel.Because the Gaussian decays exponentially, a small interval around zero can be used for discretization.Source: Perturbation perspective · p. 393, Limit case · p. 394, §14.1.3 · pp. 394–395
14.1.3 Discrete convolution
- Convolution supports finite-sum discrete computation, differentiation under integrability conditions, and multidimensional Gaussian smoothing, including the image-processing operation called Gaussian blur.If g has finite support, the discrete convolution becomes a finite summation; compactly supported n-times differentiable kernels yield n-times differentiable convolutions.Source: §14.1.3 · p. 394, §14.1.4 · p. 395, §14.1.5 · p. 396
14.1.6 Link between convolution and infimal convolution
- Classical convolution uses sum-product algebra, whereas infimal convolution uses min-plus algebra; soft infimal convolution connects them through a soft minimum.The soft Moreau envelope is computationally useful because it can be approximated using a discrete convolution, but it is not the same as Gaussian convolution.Source: §14.1.6 · p. 396, §14.1.7 · p. 396, §14.1.8 · p. 398
Computation using a convolution
- A discrete convolution on a grid of size n computes the soft infimal convolution in O(n log n) after an exponential change of variables.The reduction maps functions to the exponential domain, applies a classical convolution, and then returns to the log domain.Source: Computation using a convolution · p. 397
14.2.1 Convolution theorem
- In the Fourier domain, convolution becomes multiplication, and the same convolution theorem also holds for Laplace and bilateral Laplace transforms.The Fourier identity is F{f ∗g} = F{f} · F{g}.Source: §14.2.1 · p. 399
14.2.2 Link between Fourier and Legendre transforms
- The Legendre-Fenchel transform is to infimal convolution what the Fourier transform is to classical convolution, and soft conjugates smooth the conjugate using log-sum-exp.The analogy between Fourier and Legendre transforms is summarized in Table 14.1.Source: §14.2.2 · p. 399, §14.2.3 · p. 400, Table 14.1 · p. 400
Computation using a convolution … With respect to arbitrary norms
This section develops convolutional, probabilistic, and Gumbel-based methods for smoothing, sampling, and differentiating functions and discrete or structured choices. It also covers gradient estimation, optimization objectives, variational updates, function regularity, and related computational and estimator limitations.
Computation using a convolution
- Replacing the continuous domain with a grid enables an O(n log n) discrete-convolution approximation of the smoothed convex conjugate, but only for low-dimensional sets because of the curse of dimensionality.The method provides a GPU-friendly alternative to the fast Legendre transform and can compute an approximate biconjugate, or convex envelope.Source: Computation using a convolution · p. 402
14.3.1 Smoothed step function
- Gaussian convolution gives analytic smoothing formulas for the Heaviside and ReLU functions, with the smoothed Heaviside expressed through a Gaussian CDF and the error function available in numerical software.For the ReLU, differentiating the smoothed function yields the smoothed Heaviside, and differentiating again yields the Gaussian kernel κσ.Source: §14.3.1 · pp. 403–404, §14.3.2 · pp. 404–405
14.4.1 Expectation in a location-scale family
- Smoothing by additive noise is equivalent to convolution, and the perturbation strength σ controls recovery of the original function as σ →0.The perturbed function can be represented using a location-scale pushforward distribution and estimated by sampling noisy inputs and averaging function values.Source: §14.4.1 · pp. 405–406, Approximation by Monte-Carlo estimation · p. 407
14.4.2 Gradient estimation by reparametrization
- Reparametrization provides a gradient estimator for smoothed functions when differentiation and integration can be interchanged, but it may fail when the function is not absolutely continuous.The formula can still apply when f is differentiable almost everywhere, as for ReLU, whereas for the Heaviside function it gives zero instead of the smoothed gradient.Source: §14.4.2 · p. 407
14.4.3 Gradient estimation by SFE, Stein’s lemma
- The score-function estimator handles cases lacking ∇f or absolute continuity; for Gaussian noise, Stein’s lemma gives an unbiased estimator of the smoothed gradient but a biased estimator of the original gradient.Smoothing can nevertheless accelerate convergence of gradient-based algorithms.Source: §14.4.3 · pp. 407–408
14.4.4 Link between reparametrization and SFE
- Integration by parts converts the reparametrization-trick estimator into the score-function estimator, using gradients with respect to the noise variable rather than distribution parameters.The identity assumes that h(z)p(z) vanishes as ∥z∥→∞.Source: §14.4.4 · pp. 408–409
14.4.5 Variance reduction and evolution strategies
- Variance-reduced score-function estimators use forward or central finite differences, while the zero-temperature limit σ →0 recovers the randomized forward-mode estimator.Antithetic sampling evaluates both Z and −Z, and the variance-reduced estimators are described as more accurate and numerically stable as σ decreases.Source: Figure 14.4 · p. 410, §14.4.5 · p. 410, §14.4.6 · p. 411
14.5.2 Perturbed comparison
- Gumbel perturbations turn binary comparisons into logistic relaxations: independent shifted Gumbel noise gives E[gt(µ1 + σZ1, µ2 + σZ2)] = 1/(1 + e−(µ1−µ2)/σ).The difference of the perturbed variables follows a Logistic(µ1 − µ2, σ) distribution.Source: §14.5.2 · pp. 413–414
14.5.3 Perturbed argmax / 14.5.4 Perturbed max
- Gumbel perturbations yield closed-form expectations for perturbed argmax and max, with the latter satisfying fσ(µ) = σ · LSE(µ/σ).For perturbed argmax, the expectation is the softargmax; for perturbed max, the expectation equals a scaled log-sum-exp.Source: §14.5.3 · p. 414, §14.5.4 · p. 416, §14.6 · p. 421
14.5.5 Gumbel trick for sampling
- The Gumbel trick samples from Categorical(softargmax(µ/σ)) in one pass by taking an arg-maximum without explicitly computing normalized probabilities.The analogous Gaussian construction lacks a closed form because it requires integrating the maximum of M −1 Gaussian distributions.Source: §14.5.5 · p. 416
14.5.6 Perturb-and-MAP
- In structured prediction, perturbed objectives generally lack closed forms and are approximated by Monte Carlo; reparameterization-based gradient estimation usually has lower variance than score-function estimation.For Jacobians, score-function and differentiation-based estimators are available, but reparameterization cannot be used because the argmax is discontinuous.Source: §14.5.6 · pp. 417–418
Link between perturbation and regularization
- Under convex-polytope and positive-density assumptions, perturbation corresponds to regularization through a strictly convex function whose gradient is a mirror map, although the converse generally fails.The associated convex-conjugate regularizer does not generally have a closed form, and not every regularization has a corresponding noise distribution.Source: Link between perturbation and… · p. 418
14.5.7 Gumbel-softargmax
- Gumbel-softargmax replaces the discontinuous argmax with a temperature-controlled softargmax, enabling unbiased gradients for the smoothed expectation while introducing temperature-dependent bias relative to the original objective.As τ →0, hσ,τ(µ) → hσ(µ), but g must be defined on the simplex; the resulting continuous random variable has a Concrete distribution.Source: §14.5.7 · pp. 419–420
14.6 Summary
- The chapter summarizes convolutional smoothing as expectation under additive noise, with gradients estimated by path-gradient or score-function methods depending on whether the original gradient is available.For Gumbel noise, closed forms arise for perturbed comparison, argmax, and max, while infimal convolution provides the optimization-side counterpart of convolution.Source: §14.6 · pp. 420–421
15.1 Objective functions / 15.2 Oracles / 15.3 Variational perspective of optimization algorithms
- Optimization seeks approximate minimizers of an objective, using iterative algorithms that access zero-, first-, or second-order oracles and update parameters from oracle information and hyperparameters.A first-order example is gradient descent, while the proximal point method balances objective reduction against staying near the current iterate but is generally impractical to solve exactly.Source: §15.1 · p. 423, §15.2 · p. 424, §15.3 · p. 425
With respect to arbitrary norms / 15.4.2 Smooth functions
- The section formalizes Lipschitz continuity with arbitrary norms and domains, and defines smoothness as Lipschitz gradients that provide local quadratic upper and lower bounds.When the norm is the 2-norm, the arbitrary-norm definition recovers the standard one; smoothness makes nearby gradient information meaningful.Source: With respect to arbitrary norms · p. 426, §15.4.2 · pp. 426–427
With respect to arbitrary norms … Strongly convex case
This section develops convergence guarantees and optimization methods across convexity regimes, from non-convex to strongly convex objectives, including arbitrary norms. It examines gradient-based, stochastic, constrained, composite, momentum, and proximal approaches, emphasizing the roles of assumptions, variance, stepsizes, projections, and proximal operators.
15.4.3 Convex functions
- A zero-gradient point is a global minimum for differentiable convex functions, and smooth convexity yields sublinear gradient-descent convergence, while strong convexity enables linear convergence.For twice differentiable functions, convexity corresponds to a positive-semidefinite Hessian, whereas strict convexity corresponds to a positive-definite Hessian.Source: §15.4.3 · p. 429, §15.4.4 · p. 430
15.4.4 Strongly-convex functions
- Strong convexity lower-bounds a differentiable function by a quadratic, providing the curvature condition underlying linearly convergent algorithms.The definition and Hessian characterization extend this condition through a positive parameter µ for differentiable and twice differentiable functions.Source: §15.4.4 · p. 430
With respect to arbitrary norms
- Strong convexity extends to arbitrary norms, and specifying the norm and domain can improve the constant µ and thus convergence rates.Strong convexity with respect to one norm implies strong convexity with respect to another because all norms are equivalent, but the constants may differ.Source: With respect to arbitrary norms · pp. 430–431
15.4.5 Nonconvex functions
- For nonconvex functions, a zero-gradient point is only stationary and may be a local maximum or minimum, so gradient methods generally guarantee convergence only to vanishing gradient norm.A zero gradient is necessary for a minimum, but without convexity it does not certify a global minimum.Source: §15.4.5 · p. 432
15.5. Performance guarantees
- Convergence performance is expressed either by iterations needed for ε-accuracy or by an accuracy function R(t), whose decay defines sublinear, linear, or superlinear rates.Sublinear rates include O(1/t^α), linear rates include O(exp(−ct)), and superlinear rates have asymptotic progress ratio 0.Source: §15.5 · pp. 433–434
Zero-order vs. first-order
- For smooth strongly convex functions, the best first-order method is P times better than the best zero-order method up to constant and logarithmic factors, especially when P ≫1.First-order rates can be dimension independent, whereas zero-order rates generally depend on dimension and are therefore unsuitable for high-dimensional problems.Source: Zero-order vs. first-order · pp. 435–436
16.1.1 Variational perspective
- Gradient descent can be written as minimizing a linearized objective plus a proximity term, balancing progress from the current gradient against staying near the current iterate.The proximity term is 1/(2γ)∥w − wt∥2_2, and sufficiently small γ keeps the update within a region where the linearization is valid.Source: §16.1.1 · pp. 437–438
16.1.2 Convergence for smooth functions
- For smooth functions, each gradient-descent step decreases the objective by a term proportional to the current gradient norm, so sufficiently many iterations reach an arbitrarily small gradient norm.This guarantee follows by summing the per-iteration decrease and telescoping the objective values.Source: §16.1.2 · p. 438, §16.1 · p. 439
Non-convex case
- Gradient descent requires O(ε−2) iterations to reach an ε-stationary point in the non-convex case.Here ε-stationarity means ∥∇L(ŵ)∥2 ≤ ε.Source: Non-convex case · p. 439
Convex case
- For β-smooth convex functions with stepsize γ ≤ 1/β, gradient descent has sublinear convergence and complexity T(ε) = O(1/ε).Convexity ensures convergence to a minimum.Source: Convex case · p. 439
Strongly convex case
- For β-smooth, µ-strongly convex functions, gradient descent converges linearly with complexity T(ε) = O(ln ε−1).A convex gradient-dominating function can also be minimized at a linear rate.Source: Strongly convex case · pp. 439–440
16.1.3 Momentum and accelerated variants
- Momentum adds velocity to reduce gradient-descent oscillations, while heavy-ball achieves a linear rate that is better than gradient descent’s O(exp(−Tµ/β)) rate for suitable parameters.For ν = 0, heavy-ball recovers gradient descent; for ν > 0, accumulated velocity provides inertia.Source: §16.1.3 · p. 440, §16.2 · p. 441
16.2.1 Stochastic gradients
- Stochastic gradients provide unbiased estimates of the expected gradient, and mini-batches preserve this expectation while reducing computation from n samples to m ≪ n.Averaging multiple stochastic gradients can also reduce approximation error.Source: §16.2.1 · p. 442, §16.2 · p. 442
16.2.2 Vanilla SGD
- Vanilla stochastic gradient descent converges to the expected objective up to a variance- and stepsize-dependent constant, while mini-batches reduce variance from σ2 to σ2/m.Decreasing stepsizes such as γt = t−1/2 further control the additional error term.Source: §16.2.2 · pp. 443–444
16.2.3 Momentum variants
- Stochastic optimization extends deterministic acceleration through stochastic momentum and Nesterov momentum, while adaptive methods use gradient moments for preconditioning and stepsize robustness.Line-search rules can also select stepsizes when the smoothness constant is unknown; Adam is a prominent adaptive optimizer.Source: §16.2.3 · p. 444, §16.2.4 · p. 445
16.3 Projected gradient descent
- Projected gradient descent preserves feasibility through projection and retains gradient descent’s convergence rates, while proximal gradient handles composite objectives with convex, potentially non-differentiable terms at O(1/T) for smooth convex functions.Proximal operators are especially useful for sparsity regularization such as LASSO, and projected methods rely on available constraint oracles.Source: §16.3.1 · p. 446, §16.3 · p. 447, §16.4 · p. 448
16.4.3 Commonly-used proximal operators
- With commonly used proximal operators, proximal gradient recovers gradient descent or projected gradient descent.Source: §16.4.3 · p. 450
- Coordinate-wise operations yield the soft-thresholding operator, while group-lasso regularization encourages group sparsity.Source: §16.4.3 · p. 450
16.5 Summary
- Variationally, gradient descent linearizes the objective with quadratic regularization, while projected gradient descent adds a constraint and Euclidean projection.Source: §16.5 · pp. 450–451
- Proximal gradient linearizes the differentiable term while retaining the non-differentiable term, and stochastic methods include vanilla SGD, momentum, and Adam.Source: §16.5 · p. 451
Second-order optimization … Properties
This section develops second-order optimization methods and efficient preconditioning strategies, alongside foundational duality and divergence concepts. It emphasizes curvature-aware updates, computationally tractable approximations, and structural relationships among conjugates, divergences, and losses.
17.1. Newton’s method
- Newton’s method minimizes a local quadratic Taylor approximation, but requires a positive-definite Hessian; otherwise the quadratic minimum may be undefined or correspond to a maximum.When the Hessian is positive definite, the quadratic minimum is unique; regularization can restore positive definiteness.Source: §17.1 · p. 453
17.1.2 Regularized Newton method
- Regularized Newton methods add a proximity term and choose ηt > 0 so that ∇2L(wt) + ηt I is positive-definite, while linesearch is used because Newton updates may not decrease objective values.Source: §17.1.2 · p. 453, §17.1 · p. 455
17.1.3 Approximate direction
- Exact Newton directions require an inverse-Hessian vector product obtained by solving a linear system that a priori takes O(P^3) time, motivating approximate IHVPs.Newton’s geometric advantage is curvature adaptation: for any strongly convex quadratic, one Newton step reaches the optimum, whereas gradient descent may require many iterations.Source: §17.1.3 · p. 454, §17.1.6 · p. 456
17.1.4 Convergence guarantees
- Newton’s method can converge quadratically when initialized sufficiently close to a minimizer of a strongly convex function with Lipschitz-continuous Hessian.Under the stated neighborhood condition, the rate is R(t) ≤ O(exp(−2t)), although Newton steps cost more computationally.Source: §17.1.4 · p. 454
17.2. Gauss-Newton method
- Gauss–Newton linearizes the inner function in a composite objective with convex outer loss, producing a convex subproblem even when the original objective is nonconvex.A quadratic outer-loss approximation yields an always-convex quadratic objective, and regularization or linesearch can provide descent directions; stochastic variants estimate gradients and Gauss–Newton matrices from samples or minibatches.Source: §17.2.1 · p. 458, §17.2.2 · p. 459, §17.2.3 · p. 460
17.3 Natural gradient descent
- Natural gradient descent replaces Euclidean proximity with KL divergence between parameterized distributions, whose quadratic approximation yields a Fisher-information-based update.Gauss–Newton and natural-gradient methods coincide when p_θ is an exponential-family distribution; stochastic versions estimate gradients and Fisher matrices from sampled data.Source: §17.3 · p. 461, §17.3.1 · p. 462, §17.3.2 · p. 462
17.4. Quasi-Newton methods
- Quasi-Newton methods approximate inverse Fisher or Gauss–Newton products to reduce deep-learning costs, while BFGS builds a quadratic model from iterates satisfying a secant condition.KFAC uses an efficient structural approximation, and BFGS constructs curvature information without directly using the Hessian.Source: §17.4 · p. 463, §17.4.1 · p. 463
17.4.2 Limited-memory BFGS
- LBFGS avoids storing the full P×P BFGS matrix by applying it as a linear map using a small history, reducing cost while retaining second-order-like benefits.Storing B_t ∈ R^{P×P} is prohibitive when P is large; LBFGS therefore stores past (s_t, y_t, ρ_t) tuples and is effective for medium-scale unstructured problems.Source: §17.4.2 · p. 464, §17.4.1 · p. 464
17.5 Approximate Hessian diagonal inverse preconditioners
- Diagonal inverse-Hessian approximations provide cheap preconditioners, and algorithms such as ADAM can be interpreted as using them.Such scaling is sufficient for the quadratic example discussed in the chapter.Source: §17.5 · p. 464
17.6 Summary
- The optimization methods reviewed differ in their approximation models, proximity terms, and trade-offs between evaluating B_t∇L(w_t) and using richer information about L.The variational viewpoint frames each next iterate as balancing an approximation of the objective with a proximity term.Source: §17.6 · pp. 464–465
18.1 Dual norms
- Dual norms are support functions of norm-induced unit balls and generalize Cauchy–Schwarz inequalities; examples include the 2-norm being self-dual and the 1- and ∞-norms being dual.Source: §18.1 · p. 466, §18.2 · p. 467
18.2 Fenchel duality
- Fenchel duality gives a concave lower-bounding dual for potentially nonlinear and nonconvex primal objectives, with convex conjugate terms that are generally difficult to compute.The result is fully general, while the conjugate-based terms are computationally challenging when f is nonlinear or the loss and regularizer are nonconvex.Source: §18.2 · pp. 467–468
- Under a linear map and convex loss and regularizer, strong duality makes primal and dual optimum values equal, enabling solution in the smaller space M when M is smaller than W.Strict convexity of R additionally makes the primal solution unique.Source: §18.2 · p. 469
Properties
- Bregman divergences measure the gap between a function and its linearization, are positive under strict convexity, convex in their first argument, and admit a dual-space form.They generalize squared Euclidean and generalized Kullback–Leibler divergences and define Bregman projections and proximal operators.Source: §18.3 · p. 471, Properties · pp. 471–472
18.4 Fenchel-Young loss functions
- Fenchel–Young losses connect convex conjugates, primal-dual Bregman divergences, and prediction targets, with convexity in θ and conditional convexity in t; softmax is a KL projection onto the probability simplex.For categorical outputs, targets can use one-hot encodings, and positivity characterizes equality through ∇Ω*(θ) = t when Ω is strictly convex.Source: §18.4 · pp. 473–474, §18.5 · p. 474