Source-linked AI summary
The Stan Math Library: Reverse-Mode Automatic Differentiation in C++
Bob Carpenter, Matthew D. Hoffman, Marcus Brubaker, Daniel Lee, Peter Li, Michael Betancourt
TL;DR
Algorithms in optimization and statistical inference increasingly require derivatives, while implementing them is a substantial user burden. The Stan Math Library provides reverse-mode autodiff through direct and functional interfaces, with performance varying across expression types and further gains indicated for larger matrix evaluations.
Problem
Many contemporary algorithms require derivatives, and implementing those derivatives is a substantial user burden.
Source: §1 · p. 2Method
The Stan Math Library uses reverse-mode automatic differentiation and offers both direct data-structure use and a gradient functional for templated C++ functors.
Source: §18 · p. 78Results
Performance comparisons varied by problem: Stan and Sacado were nearly identical for basic function and operator evaluations, while matrix evaluations indicated further performance gains for larger matrices.
Source: Power Function · p. 69, §18 · p. 79Takeaways & Limitations
The library's examples demonstrate using its data structures directly or encapsulating gradient calculations for templated C++ functors with Stan's gradient functional.
Source: §18 · p. 78Adept's expression templates can make it faster than Stan at run time for complex assignment expressions, with the advantage growing as the right-hand-side expression grows.
Source: §18 · p. 79
Abstract
from arXiv · showhide
As computational challenges in optimization and statistical inference grow ever harder, algorithms that utilize derivatives are becoming increasingly more important. The implementation of the derivatives that make these algorithms so powerful, however, is a substantial user burden and the practicality of these algorithms depends critically on tools like automatic differentiation that remove the implementation burden entirely. The Stan Math Library is a C++, reverse-mode automatic differentiation library designed to be usable, extensive and extensible, efficient, scalable, stable, portable, and redistributable in order to facilitate the construction and utilization of such algorithms. Usability is achieved through a simple direct interface and a cleanly abstracted functional interface. The extensive built-in library includes functions for matrix operations, linear algebra, differential equation solving, and most common probability functions. Extensibility derives from a straightforward object-oriented framework for expressions, allowing users to easily create custom functions. Efficiency is achieved through a combination of custom memory management, subexpression caching, traits-based metaprogramming, and expression templates. Partial derivatives for compound functions are evaluated lazily for improved scalability. Stability is achieved by taking care with arithmetic precision in algebraic expressions and providing stable, compound functions where possible. For portability, the library is standards-compliant C++ (03) and has been tested for all major compilers for Windows, Mac OS X, and Linux.
1. Reverse-Mode Automatic Differentiation
This section explains reverse-mode automatic differentiation, including how it computes gradients and how Stan exposes it through direct and functional C++ interfaces. It also describes functionals for differentiating templated functors and calculating Jacobians by applying a function once and computing a gradient for each output dimension.
1. Reverse-Mode Automatic Differentiation
- Because gradients are prevalent, reverse mode is often more efficient: one backward pass computes a single output’s derivatives with respect to multiple inputs.This gives the gradient of a single-output function or one row of a multi-output function’s Jacobian; forward mode instead computes a directional derivative in one pass.Source: §1 · p. 3
1.1. Mechanics of Reverse-Mode Automatic Differentiation
- Reverse mode builds a topologically ordered expression stack, then sweeps backward, accumulating each node’s adjoint times its local partials at its operands.The root adjoint starts at 1 and other adjoints at 0; Stan usually evaluates partials lazily during the reverse pass, and input adjoints yield the gradient.Source: §1.1 · pp. 4, 6
1.2. Comparison to Alternative Methods of Computing Gradients
- Finite differences require N + 1 evaluations (2N + 1 centered) and are slow and imprecise, while symbolic differentiation can optimize formulas but is less expressive.Finite-difference precision is usually at best 10^-7 versus roughly 10^-14 for double precision; symbolic methods can handle repeated subexpressions efficiently but require multiple passes or compilation.Source: Finite Differencing · p. 6, Symbolic Differentiation · p. 7
2. Calculating Gradients and Jacobians
- Stan’s reverse-mode library evaluates gradients for functions from R^N to R and Jacobians for differentiable functions from R^N to R^M at a specified input.Source: §2 · p. 8
2.1. Direct Calculation of Gradients of Programs
- In direct C++ programs, assigning differentiable inputs to stan::math::var lets ordinary operations build the expression graph; grad() propagates derivatives and adj() retrieves them.The resulting gradient can also be collected in a standard vector with grad().Source: §2.1 · pp. 8–9
2.2. Coding Template Functions for Automatic Differentiation
- Templating functions across argument types lets them accept primitive and stan::math::var inputs, with promoted return types becoming var whenever any input is var.Argument-dependent lookup resolves autodiff overloads alongside primitive functions, while Boost promotion traits determine return and intermediate types.Source: §2.2 · p. 9, Traits Metaprogram for Computing… · p. 10, Argument-Dependent Lookup for Function… · p. 10
2.3. Calculating the Derivatives of Functors with Functionals
- Stan provides a functional interface that represents a function as a C++ functor and applies a gradient operator to it.It uses Eigen for matrix operations and linear algebra, including dynamically allocated vectors and matrices whose elements can be stan::math::var.Source: §2.3 · p. 11, Eigen C++ Library for Matrices and… · p. 11
- stan::math::gradient() evaluates a functor at an input vector and returns both its function value and gradient.Templated operators can use stan::math::var for differentiation and double for testing.Source: Calculating Gradients with Functionals · pp. 12–13
2.4. Calculating Jacobians
- A single forward pass constructs the expression graph, after which reverse passes calculate gradients for each output dimension of the Jacobian.The direct procedure converts inputs to autodiff variables, applies the function, and uses each output's gradient to populate a Jacobian row; subsequent gradients reset adjoints to zero.Source: §2.4 · p. 13, Direct Jacobian Calculation · p. 13
- The jacobian() functional applies directly to the function and automatically resizes the Jacobian and function-value vector when needed.It provides an alternative to constructing the Jacobian directly.Source: Functional Jacobian Calculation · p. 14
3. Automatic Differentiation Variable Base Classes
The section presents `stan::math::var` as the client-facing type for autodiff and explains how implementation pointers, scoped resource management, and chainable variable implementations support reverse-mode differentiation.
3. Automatic Differentiation Variable Base Classes
- `stan::math::var` replaces primitive doubles in functions to be differentiated, while its implementation is held through a pointer to `stan::math::vari`.The Pimpl pattern keeps implementation details out of the client-facing API.Source: §3 · p. 14, §3.1 · p. 14, §3.3 · p. 15
3.1. Pointers to Implementations
- Pimpl encapsulates implementation details, while RAII manages `var` instances on the stack and frees their memory when they go out of scope.Memory is allocated when an object is constructed or resized; `vari` storage uses a custom arena rather than the general C++ heap.Source: §3.1 · p. 14, §3.2 · p. 14
3.3. The var Class
- The `var` class points to a `vari` implementation, exposes its value and adjoint, and provides constructors and assignment operators for primitive types.A default-constructed `var` uses a null implementation pointer to avoid an allocation that would be overwritten.Source: §3.3 · p. 15
3.4. The chainable Base Class
- The `chainable` base class lets function-specific subclasses implement virtual `chain()` methods for derivative propagation and allocates their instances in an arena.Virtual initialization and zero-adjoint methods support setting adjoints during reverse propagation.Source: §3.4 · pp. 15–16
3.5. The vari Class
- Each `vari` stores a value and adjoint, and its constructor adds the instance to the variable stack.Its dependent-variable initializer sets the output adjoint to 1; other adjoints are set to zero when multiple derivatives, such as Jacobians, are required.Source: §3.5 · p. 16
4. Calculating Gradients
Stan computes gradients by initializing the output adjoint and traversing the expression stack in reverse topological order. It reuses that graph for Jacobians and manages graph memory with an arena that can be reset and reused.
4. Calculating Gradients
- Reverse traversal calls each node’s chain() method after its parents, propagating derivatives through the expression graph.The dependent variable’s adjoint is initialized to 1 before the traversal.Source: §4 · pp. 16–17
4.1. Jacobians and Zeroing Adjoints
- For each Jacobian row, Stan reuses the expression graph and zeros adjoints before every gradient calculation after the first.Zeroing walks the variable stack and resets each adjoint.Source: §4.1 · p. 17
5. Memory Management
- The arena-based allocator supplies memory for variable instances, and the variable stack records their expression-graph pointers for reverse traversal.The allocator and stack are held in static storage.Source: §5 · pp. 17–18
5.2. Memory Lifecycle
- The expression graph’s memory is reclaimed together after reverse propagation, while the arena’s underlying blocks are retained for reuse.The recovery method resets the arena and pointer stack; the functionals also recover memory when an operation throws.Source: §5.1 · p. 18, §5.2 · p. 20
5.1. Byte-Level Memory Manager: stack alloc
- The arena stores autodiff instances in blocks, allocating another block when a request will not fit in the current one.A new block is twice the current block’s size; allocation is 8-byte aligned.Source: §5.1 · pp. 18–19
- Nearly 10%: compiler branch-prediction hints improved allocator throughput compared with unadorned code.The hints mark branches as likely or unlikely.Source: §5.1 · p. 19
6. Gradient Functional
The gradient functional wraps a user functor in an interface that returns the function value and gradient. It also recovers autodiff memory on both normal and exceptional exits.
6. Gradient Functional
- The functional recovers memory after evaluation and also recovers it before rethrowing a caught standard exception.The underlying memory pool is retained for reuse.Source: §6 · p. 21
- The functional converts the input doubles to autodiff variables, evaluates the functor, propagates derivatives, and extracts the gradient.It returns the function value through fx and the gradient through grad_fx.Source: §6 · pp. 20–21
- The functor must accept a vector of var values and return a var; other required data or callbacks must be accessible to it.Source: §6 · p. 21
7. Constants
Stan provides two equivalent ways to create constant autodiff variables, though constants do not propagate derivatives. The paper recommends using doubles directly where possible to reduce the variable stack.
7. Constants
- A constant var can be created either with the unary constructor or by assigning a double; both have the same effect on memory and the autodiff stack.Source: §7 · pp. 21–22
- Constructors also support integers and other base types, while the explicit copy constructor avoids ambiguity.Source: §7 · p. 22
- Constants inherit a no-op chain() method and propagate no derivatives, so using doubles directly where possible reduces the variable stack and speeds gradient calculations.Source: §7 · p. 22
8. Functions and Operators
Stan implements functions and operators with specialized autodiff variable classes that store values and operands and propagate derivatives during the reverse pass. Lazy partial evaluation and special cases reduce stored work for selected operations.
8.1. Unary Operand Storage Class
- Unary operand storage keeps the function value and operand pointer in a helper class, while unary functions specialize chain() to propagate derivatives.Source: §8.1 · p. 22, §8.2 · p. 23
8.2. Unary Function Implementation
- Lazy unary derivatives are evaluated during the reverse pass, avoiding stored partials and an extra multiplication; addition and subtraction can avoid multiplication altogether.Source: §8.2 · p. 23
8.3. Binary Functions
- Binary-function classes store operand pointers; mixed var/double operations propagate only one derivative, and pow(double exponent) uses sqrt, identity, or square special cases when applicable.Source: Binary Operands Storage Class · p. 24, Binary Function Implementation · p. 25, Variable Implementations · p. 26
- For pow(var, var), a zero base returns before evaluating log(0), which would produce a not-a-number value; otherwise chain() updates both operand adjoints.Source: Variable Implementations · pp. 25–26
8.4. Arithmetic Operators
- Addition propagates by adding the expression adjoint to each operand’s adjoint, without storing or multiplying derivative values of 1.Source: §8.4 · pp. 26–27
8.5. Boolean Operators
- Boolean operators are implemented directly without creating new variable implementation instances.Source: §8.5 · p. 27
8.6. Memory Usage
- Each expression’s memory cost is 24 + K × 8 bytes for K operands, including its value, operand pointers, vtable pointer, and variable-stack pointer.Source: §8.6 · pp. 27–28
9. Assignment and Compound Assignment
Assignment copies the implementation pointer for var values, constructs a var from doubles, and forbids assigning var expressions to doubles. Compound assignment replaces the left-hand var’s implementation pointer with a new expression node and returns the left-hand variable by reference.
9.1. Assignment Operator
- A var-to-var assignment shallow-copies its vari* pointer, while assigning a double constructs a var whose pointer is copied to the left-hand side.Assigning a var expression to a double is forbidden.Source: §9.1 · p. 28
9.2. Compound Assignment Operators
- Compound operators are member functions with the left-hand variable implicit, and their double-argument implementations parallel the variable-argument versions.The implementations are defined outside the class declaration.Source: §9.2 · p. 28
- A compound assignment creates a new multiply vari, stores its pointer in the left-hand var, and returns that var by reference.For +=, the left-hand var is assigned the result of a + b; this changes the var instance, not the variable implementation.Source: §9.2 · pp. 28–29
10. Variadic Functions
Variadic implementations represent operations over arbitrary-length inputs with operand arrays and use arena allocation. Specialized log-sum-exp and sum nodes provide compact derivative propagation, with the sum node avoiding the repeated binary additions of a running accumulation.
10. Variadic Functions
- Variadic functions such as sum use a base class that stores an operand array and its size, allocated in the autodiff memory arena.Arena allocation avoids fragmentation and the overhead of calling destructors associated with standard-heap vectors.Source: §10 · p. 29
10.1. Log Sum of Exponentials
- The log-sum-exp operation is symmetric, propagates each operand’s derivative through one chain call, and uses an algebraic trick to avoid overflow.Its derivative with respect to an operand is described in terms of that operand and the log-sum-exp result.Source: §10.1 · p. 30
10.2. Sum Accumulation
- A single N-operand sum node uses one chain call and N additions, instead of N−1 binary nodes, N−1 chain calls, and 2(N−1) additions.The library uses this approach to accumulate intermediate log-density values; the passage contrasts it with a left-associated running-sum expression tree.Source: §10.2 · pp. 31–32
- The sum node stores operand pointers and their count, and the vector sum function returns a var backed by that specialized implementation.Special cases for inputs of size zero, one, or two can be added when common.Source: §10.2 · pp. 31–32
11. Matrix Data Structures and Arithmetic
Stan uses Eigen for matrix storage, arithmetic, and linear algebra, with matrices and vectors organized around univariate autodiff variables. Typed operations infer result shapes, while specialized expression nodes reduce graph and storage costs for dot products, matrix multiplication, and self-transpose multiplication.
11. Matrix Data Structures and Arithmetic
- Stan uses Eigen for matrix data structures, arithmetic, and linear algebra operations such as factoring and determinants.Source: §11 · p. 32
11.1. Matrix Data Structures
- Matrix and vector containers hold univariate var instances, while their element type can be double or var.The design is organized around univariate automatic-differentiation variables.Source: §11.1 · p. 33
11.2. Type Inference
- Typed arithmetic infers output shapes, and adding a scalar to a matrix or vector adds it elementwise.Source: §11.2 · p. 33
11.4. Matrix Multiplication
- Specialized matrix multiplication reduces node allocations from 2N to 1 and derivative-propagation edges from 4N to 2N.For multiplying two N × N matrices, the operand-storage scheme requires O(N^2) memory rather than O(N^3) for separate storage per result entry.Source: §11.3 · p. 33, §11.4 · p. 35
11.3. Dot Products
- A custom var-by-double dot product avoids promoting constants to var nodes, using one expression node and N derivative propagations for N-vectors.Promoting the double entries would introduce N unnecessary nodes and gradient propagations to constants.Source: Mixed Type Operations · p. 35
11.5. Specialized Matrix Multiplication
- A custom matrix-by-self-transpose operation halves memory use and improves speed over transposing first and then using general multiplication.For lower-triangular inputs, only nonzero portions need to be stored and traversed; the operation reconstructs Ω = L L⊤, as used for correlation, covariance, or precision matrices.Source: §11.5 · pp. 35–36
12. Linear Algebra
For linear algebra operations, specialized derivatives can be more efficient than automatically differentiating library functions. Stan’s log-determinant implementations compute the value and derivatives using matrix decompositions and precomputed gradients.
12. Linear Algebra
- Specialized derivatives for linear algebra operations can be more efficient than automatically differentiating library functions.Source: §12 · p. 36
12.1. Log Determinants
- The general log-absolute-determinant implementation uses Householder QR decomposition to obtain the value and inverse, then stores precomputed gradients for the reverse pass.Matrix calculations are performed on double values to reduce expression-graph size and reverse-pass propagation time.Source: Precomputed Gradients Implementation · p. 36, Log Determinant Implementation · pp. 37–38
- For symmetric positive-definite matrices, Stan provides a specialized log-determinant implementation using LDL⊤ factorization to efficiently extract the determinant and inverse.Source: Specialized Log Determinant… · p. 38
13. Differential Equation Solver
The solver integrates ODE systems from specified initial states and supports derivatives of solutions with respect to parameters and initial values. It computes sensitivities through an augmented coupled system rather than differentiating the numerical integrator directly.
13.2. Sensitivities of Solution to Inputs
- Solution sensitivities support optimization, parameter estimation, and reporting variation with respect to parameters.The inputs may be parameters, initial state, or both.Source: §13.2 · p. 40
13.3. Differentiating the Integrator
- Differentiating the numerical integrator can create a large, time- and memory-intensive expression graph and offers no control over sensitivity error.Limited integrator templating may require promoting all inputs, worsening the expression-graph cost.Source: §13.3 · p. 40
13.4. Coupled System
- Adding sensitivity state variables to the original ODE yields a coupled system whose stepwise integration allows sensitivity error to be controlled.The coupled system size depends on whether derivatives are required with respect to initial states, parameters, or both.Source: §13.4 · pp. 40–41
13.5. Differentiating an ODE Solution with Respect to a Parameter
- Parameter sensitivities are added as states with zero initial values, and their time derivatives complete a system coupled to the original state.The system Jacobian with respect to parameters is needed to add these sensitivities.Source: §13.5 · pp. 41–42, §13.7 · p. 43
13.6. Differentiating an ODE Solution with Respect to the Initial State
- Initial-state sensitivities are computed by shifting the state by its initial condition and treating that condition as a parameter in a revised system.The revised system starts at zero, with derivatives of the initial state with respect to itself given by the identity mapping.Source: §13.6 · pp. 42–43
13.7. Computing Sensitivities with Nested Automatic Differentiation
- Nested reverse-mode autodiff computes the system Jacobians while reusing the stack above a recorded location without rebuilding the outer expression graph.Source: §13.7 · p. 43
13.8. Putting Results Back Together
- At requested solution times, the solver connects each output state variable to its inputs through stored sensitivities.Sensitivities are stored in a general precomputed-gradients variable structure.Source: §13.8 · p. 44
13.10. Integration Function
- The integration function accepts an initial state, increasing solution times later than the initial time, parameters, and separate data inputs; backward-time solutions require reversing signs.The system is supplied as a functor, and data inputs need not be promoted to autodiff variables.Source: Example: Simple Harmonic Oscillator · p. 39, §13.9 · p. 44, §13.10 · p. 45
14. Probability Functions and Traits
Probability functions can omit constant terms when calculations are needed only up to a proportion, and templated arguments let the library select which terms to compute. The resulting return type is promoted from the argument types.
- For optimization and MCMC, probability functions need only be calculated up to a constant term, so the library uses traits and template parameters to configure them.The normal log density example drops its constant term and can omit additional terms when their arguments are constant.Source: §14 · p. 46
- A leading boolean propto argument controls whether constant terms are dropped, while separately templated arguments support all constant/variable combinations and determine the promoted return type.For scalar arguments, the return-type traits behave like Boost’s promote_args.Source: §14 · p. 46
- The include_summand traits metaprogram includes a term when propto is false or when an argument is nonconstant.Static evaluation can remove unnecessary calculations, including the loop over sigma elements when its summand is not needed.Source: §14 · p. 47
15. Vectorization of Reductions
This section describes vectorized evaluation of reduction-based log joint likelihoods for shared-parameter observations, including a normal log-density implementation that handles scalar and vector inputs.
15.1. Argument Broadcasting
- Vectorized normal log-density arguments may be scalars or standard vectors and Eigen vectors, with containers matching in size and scalars broadcast across positions.Without vectorization, arguments may instead be scalar int, double, or stan::math::var values.Source: §15.1 · p. 48
15.2. Vector Views
- VectorView uses an expression template to treat scalar and container arguments uniformly, avoiding 27 separate normal log-density implementations.Its view exposes elements through indexing; scalar views return the same value for any index, while container views return indexed elements.Source: §15.2 · pp. 48–49
15.3. Value Extractor
- The value_of helper extracts a double-based value from primitive or autodiff arguments, supporting templated definitions across these types.The helper is overloaded for autodiff variables, with its definition in stan::math to allow argument-dependent lookup.Source: §15.3 · p. 50
15.4. Vector Builder Generic Container
- VectorBuilder holds reusable intermediate values, computing log σ and σ−1 once when σ is scalar and allocating storage only when needed.The inverse scale is always computed; log σ and its memory are included only when a scale-only term is present.Source: §15.4 · pp. 50–51
15.5. Operands and Partials
- OperandsAndPartials stores function operands and their partial derivatives in parallel arrays, enabling eager partial computation and chain-rule propagation.The operands and partials arrays are sized from non-constant arguments, and the chain rule updates each operand’s adjoint using the result’s adjoint and corresponding partial.Source: §15.5 · pp. 52–54
15.6. Example: Vectorization of the Normal Log Density
- The section proceeds to apply these components to define a vectorized normal log density for reverse-mode autodiff.The example continues the normal probability density function introduced in Section 14.Source: §15.6 · p. 54
- The return type is var when any argument is an autodiff variable and double otherwise; empty inputs or constant-only proportional calculations return zero.The result accumulator is initialized only after these early-return checks.Source: §15.6 · p. 55
- The function validates y, mu, and sigma and requires container arguments to have consistent sizes, while allowing scalars alongside vectors.For example, vector y and mu must have the same size when sigma is scalar.Source: §15.6 · p. 55
- Compile-time summand traits omit unnecessary log calculations, so the implementation computes only the needed inversions and logarithms before looping over the largest argument size.Scalar arguments have length 1; container arguments use their container length.Source: §15.6 · p. 56
- Analytical partial derivatives accumulate through vector views, so scalar arguments reuse one derivative slot while vector arguments accumulate per element.The same derivative-update code handles either scalar or vector sigma.Source: §15.6 · p. 57
16. Evaluation
The evaluation benchmarks Stan’s gradient computation against other reverse-mode systems across basic operators, standard functions, matrix products, and normal-log-density evaluations. It reports performance across problem sizes and implementations, while noting that memory and compile time are not measured.
16.2. What is being Evaluated
- Tests use retaping for gradients of multiple-input, single-output functions, while Stan’s lazy partial derivatives save memory for gradients but require recomputation for Jacobians.Adol-C and CppAD can reuse tapes to speed later reverse passes, but control flow can prevent reuse; memory and compile time were not evaluated.Source: §16.2 · p. 59, Gradients vs. Jacobians · p. 59, Retaping · p. 59
16.3. Functors to Differentiate
- A shared functor interface ensures each system differentiates the same code, with inputs filled once per size before timing.Functors implement an operator over a dynamic Eigen vector, plus a name and a fill function.Source: §16.3 · pp. 60–61
16.4. Functionals for Differentiation
- The system-specific functionals propagate gradients through Adept’s adjoint computation, Adol-C’s tape, CppAD’s reverse method, or Sacado’s Gradcomp call.Adept starts a recording and sets the output adjoint; Adol-C records dependent variables, while CppAD builds an ADFun for reverse evaluation.Source: Adept Gradient Functional · p. 62, Adol-C Gradient Functional · p. 62, CppAD Gradient Functional · p. 63
16.5. Test Harness
- The harness times gradient calculations for input sizes doubling from 1 through 16,384 to reveal speed trends by problem size.Each size is filled before timing, and larger sizes were omitted because 16,384 was considered sufficient to establish the trends.Source: §16.5 · p. 64
- The harness times 100,000 or 10,000 repeated calls using clock(), and includes a double-valued function-only baseline.A variable retains the result to prevent compilation from removing the calculation; only the Mac OS X Terminal was running as a user program during timing.Source: §16.5 · p. 66
16.6. Test Hardware and Compilation
- Tests used a 2.3 GHz Intel Core i7 MacBook Pro, clang 3.7.0, optimization and debugging flags, Eigen 3.2.4, and Boost 1.55.0.Stan and Adept used thread-unsafe stacks; thread-safe mode incurs roughly a 20% performance penalty, and GCC 3.9 results were similar and omitted.Source: §16.6 · pp. 66–67
16.7. Basic Function and Operator Evaluations
- For sums, Sacado is fastest below 8 dimensions, and for products below 16; Stan is faster on larger problems, while gradient evaluation takes 100 times double-function time at a handful of dimensions.At 200 or more dimensions, gradient evaluation takes roughly 10–15 times the double-function time.Source: Sums and Products · p. 67
- About 40% faster than Adept, Sacado and Stan perform almost identically on repeated power applications, while Adol-C is roughly 15 times slower than Stan.For problems above 8 dimensions, gradient calculations take only about 50% more time than double-based function evaluation.Source: Power Function · p. 69, Figure 5 · p. 69
- For summed powers with a fixed exponent, relative speeds are similar except that Adol-C is much faster; Stan’s gradients take about twice the function-only time.Source: Figure 5 · p. 69, Power Function · p. 70
- Adept is about 40% faster than Stan on the recursive log-sum-of-exponents evaluation, while Stan and Adept perform similarly with the direct implementation.Sacado is slightly faster than Stan on the recursive form, whereas Adol-C and CppAD are more than twice as slow; the direct form also speeds double-based evaluation by using fewer log and exp operations.Source: Figure 6 · p. 70, Logarithm and Exponentiation Functions · p. 71
- Stan is faster for small naive matrix products, while Adept is faster for larger ones; Sacado performs relatively poorly as expression connectivity increases.The same size-dependent Stan/Adept pattern is reported when gradients are taken only with respect to the first matrix.Source: Matrix Products · pp. 73–74, Figure 7 · p. 73
- About 50% faster, Eigen’s direct matrix-product implementation outperforms Stan’s custom implementation, but consumes roughly twice as much memory.The paper also reports gradient relative speed versus a pure double-based implementation as better by a factor of 50%.Source: Matrix Products · pp. 75–76, Figure 8 · p. 75
- Stan is roughly 50% faster than the next-fastest system, Adept, on direct normal-log-density evaluation.The advantage declines somewhat as problem size increases.Source: Normal Log Density · p. 76, Figure 9 · p. 77
- Stan’s built-in vectorized normal-log-density evaluation shows increasing speed relative to naive double evaluation as dimensionality grows.Figure 9 presents the built-in version in its second plot.Source: Figure 9 · p. 77, Normal Log Density · p. 77
- For ten or more evaluations, Stan’s gradient calculation takes a bit more than twenty times as long as double-based evaluation; Sacado falls behind as gradients increase, while CppAd approaches half Stan’s speed.Sacado is roughly as fast as Stan on small problems; CppAd starts far behind.Source: Normal Log Density · p. 77
- The custom normal-log-density version is roughly a factor of six faster than the naive implementation.Source: Normal Log Density · p. 78
17. Previous Work
Stan’s pointer-to-implementation pattern and arena-based memory design draw on RAD, which also underlies Sacado. Stan distinguishes itself with lazy partial evaluation and value-based returns.
17. Previous Work
- Stan’s pointer-to-implementation pattern and arena-based memory design were based on RAD, which also underlies Sacado.Gay and Aiken (2001) describe the arena-based memory pattern in more detail.Source: §17 · p. 78
- Stan uses vari specializations for lazy partial evaluation, saving memory and reducing assignments.Source: §17 · p. 78
- Stan’s client class var returns values, unlike RAD’s reference-based return approach for function and operator overloads.Source: §17 · p. 78
18. Summary and Future Work
The paper reports Stan’s usability and variable performance across evaluated problems, alongside memory-saving techniques and future feature plans. Its comparisons also identify cases where other approaches may be advantageous.
18. Summary and Future Work
- Examples demonstrate direct use of Stan’s data structures and gradient calculations encapsulated through its gradient functional for templated C++ functors.Source: §18 · p. 78
- Lazy evaluation and vectorization of operations such as log density functions were used to reduce memory usage.The paper also describes memory usage in detail.Source: §18 · pp. 78–79
- Stan’s library is extensively tested for accuracy and instantiability, maintains a stable interface across dozens of Stan releases, and is distributed under a BSD license compatible with licenses such as GPL.The paper attributes expected continued growth to Stan’s open-source development community.Source: §18 · p. 79
- Stan’s relative performance varied across evaluated problems, and Adept was more efficient at run time for complex right-hand-side expressions.Adept’s advantage grows with expression size; matrix evaluations also suggest larger matrices may offer further performance gains, presumably through memory locality.Source: §18 · p. 79
- CppAD and Adol-C tape reuse may help applications without general loops or conditionals that evaluate differently across runs.Source: §18 · p. 79
- Planned features include inverse cumulative distribution functions and stiff differential equation solvers, the latter requiring second-order automatic differentiation for Jacobians coupled with sensitivities.Source: §18 · p. 79
A. Stan Functions
The appendix inventories Stan’s supported functions as of version 2.6, covering type conventions, operators, mathematical and matrix operations, linear algebra, and probability functions.
A. Stan Functions
- The appendix lists supported functions as of Stan version 2.6 and explains how real arguments and return types accommodate double, int, and var values.When any argument contains a var, the result type is var; otherwise it is double.Source: A. Stan Functions · p. 80, Scalar Functions · p. 80, Vector or Array Functions · p. 80
A.2. C++ Built-in Arithmetic Operators
- Built-in arithmetic and relational operators support real-valued expressions, with relational operations returning integers.Source: §A.2 · p. 80, §A.3 · p. 80
A.4. C++ cmath Library Functions
- The listed mathematical functions include C++ cmath operations, special mathematical functions, and control or test functions such as if_else, step, is_inf, and is_nan.Source: §A.4 · p. 81, §A.5 · p. 82, §A.6 · p. 82
A.7. Matrix Arithmetic Functions
- Matrix arithmetic functions cover elementwise and matrix products, division, subtraction, dot products, and quadratic forms across matrices, vectors, and row vectors.Source: §A.7 · pp. 83–84
A.8. Special Matrix Functions
- Special matrix and array functions include summaries and transformations such as sums, means, sorting, slicing, reshaping, and conversion among arrays, matrices, vectors, and row vectors.Source: §A.8 · pp. 85–86, §A.9 · pp. 87–88
A.10. Linear Algebra Functions
- Linear algebra support includes Cholesky decomposition, determinants, eigenvalues, inverses, matrix solves, QR factors, singular values, and traces.Source: §A.10 · p. 89
A.11. Probability Functions
- Probability functions span distribution log densities and cumulative or complementary cumulative distribution functions, including multivariate normal, Student t, and Wishart functions.Source: §A.11 · pp. 90–93