Source-linked AI summary
A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions
Yingbo Ma, Vaibhav Dixit, Mike Innes, Xingjian Guo, Christopher Rackauckas
TL;DR
Computing differential-equation sensitivities is important for parameter estimation, neural differential equations, and model diagnostics, but practitioners face many competing methods. The paper benchmarks DSAAD against continuous sensitivity approaches across biological, pharmacometric, and PDE-derived models, finding that forward-mode DSAAD favors small systems while continuous adjoints scale better for sufficiently large ones. It also identifies memory, numerical stability, and reverse-mode AD implementation as important practical considerations.
Problem
Many sensitivity-analysis choices lack comprehensive comparative evidence, making it difficult to select an effective method for a particular differential-equation application.
Method
The study benchmarks DSAAD and continuous sensitivity analysis across stiff and non-stiff biological and pharmacometric models, including a PDE discretization and multiple AD modes.
Results
Forward-mode DSAAD has a strong performance advantage on sufficiently small systems, whereas continuous adjoint sensitivity analysis has an advantage on sufficiently large systems; tape-based reverse-mode AD performs poorly and scales poorly.
Takeaways & Limitations
Method choice should account for system size, memory usage, reverse-mode AD implementation, and compatibility with the solver ecosystem.
Takeaways & Limitations
DSAAD requires a general AD tool and an ODE solver implemented compatibly with that AD system, limiting direct application in ecosystems that call external Fortran solvers.
Abstract
from arXiv · showhide
Derivatives of differential equation solutions are commonly for parameter estimation, fitting neural differential equations, and as model diagnostics. However, with a litany of choices and a Cartesian product of potential methods, it can be difficult for practitioners to understand which method is likely to be the most effective on their particular application. In this manuscript we investigate the performance characteristics of Discrete Local Sensitivity Analysis implemented via Automatic Differentiation (DSAAD) against continuous adjoint sensitivity analysis. Non-stiff and stiff biological and pharmacometric models, including a PDE discretization, are used to quantify the performance of sensitivity analysis methods. Our benchmarks show that on small systems of ODEs (approximately $<100$ parameters+ODEs), forward-mode DSAAD is more efficient than both reverse-mode and continuous forward/adjoint sensitivity analysis. The scalability of continuous adjoint methods is shown to be more efficient than discrete adjoints and forward methods after crossing this size range. These comparative studies demonstrate a trade-off between memory usage and performance in the continuous adjoint methods that should be considered when choosing the technique, while numerically unstable backsolve techniques from the machine learning literature are demonstrated as unsuitable for most scientific models. The performance of adjoint methods is shown to be heavily tied to the reverse-mode AD method, with tape-based AD methods shown to be 2 orders of magnitude slower on nonlinear partial differential equations than static AD techniques. These results also demonstrate the applicability of DSAAD to differential-algebraic equations, delay differential equations, and hybrid differential equation systems, showcasing an ease of implementation advantage for DSAAD approaches.
I. INTRODUCTION
The paper frames sensitivity computation as essential for parameter estimation and related differential-equation applications, then compares discrete AD-based and continuous sensitivity approaches across their computational trade-offs. It motivates choosing methods according to system size, parameter count, memory, and differentiation strategy.
- Motivation: Model sensitivities support biological parameter tuning, neural differential-equation training, and gradient-based parameter estimation.Efficient and accurate sensitivities matter because parameter estimation requires repeated numerical solutions.
- Motivation: Numerical differentiation requires two ODE solutions per parameter and is sensitive to step-size and floating-point cancellation errors.Large perturbations increase approximation error, while small perturbations can increase cancellation error.
- Continuous sensitivity analysis: Forward-mode CSA augments the ODE system with sensitivity equations, computing sensitivities alongside the state in one solver call.The sensitivities depend on the current state, so the expanded system is solved simultaneously.
- Continuous sensitivity analysis: O(np) forward-mode CSA scaling becomes impractical with many parameters, motivating continuous adjoint methods with O(n + p) cost for solution gradients.Adjoint sensitivities are obtained by solving a backward ODE and integrating contributions to the cost-function gradient.
- Continuous sensitivity analysis: Continuous adjoint implementations trade memory against performance, while backsolve integration is reported as numerically unstable.Quadrature on a continuous adjoint solution requires memory for a continuous extension, whereas backsolve can produce divergent gradients.
- Study scope: DSAAD differentiates numerical integrator steps directly through automatic differentiation, enabling a comprehensive comparison across discrete or continuous, forward or adjoint, and tape or static choices.The study uses DiffEqSensitivity.jl to benchmark this Cartesian product of method choices.
- Study scope: The study establishes an approximate 30-100 ODE heuristic range where continuous adjoint scalability overcomes DSAAD’s low overhead.This range is presented as guidance for selecting an effective sensitivity method for a specific problem.
II. DISCRETE SENSITIVITY ANALYSIS VIA AUTOMATIC DIFFERENTIATION (DSAAD)
DSAAD combines generic, type-specialized ODE solvers with automatic differentiation to compute discrete sensitivities, and the paper tests its correctness across representative stiff, non-stiff, small, large, and PDE-derived models. The resulting sensitivities agree with continuous sensitivity analysis at numerical tolerances.
- DSAAD approach: Forward-mode DSAAD uses Dual numbers with generic ODE solvers to propagate values and derivatives through the solver’s numerical operations.Reverse-mode DSAAD instead combines the solver with reverse-mode AD Tracker types.
- Models: The benchmark models include non-stiff Lotka-Volterra and pharmacokinetic/pharmacodynamic systems, plus stiff Brusselator PDE and pollution models.The model set covers stiff and non-stiff systems, small and large systems, and a PDE discretization.
- Correctness: DSAAD and CSA sensitivity outputs align on the Lotka-Volterra and Brusselator tests, validating the AD-based solver implementation.The comparison includes the Lotka-Volterra model and a 3 × 3 discretization of the Brusselator PDE.
- Correctness: 1.14 × 10^-5 and 3.1 × 10^-4 were the maximum sensitivity differences between DSAAD and CSA at the chosen integration tolerance.With sufficiently low ODE solver tolerance, the difference was confirmed at 1 × 10^-12.
A. Forward-Mode Sensitivity Performance Comparisons
Forward-mode DSAAD performs strongly across the tested forward sensitivity benchmarks, benefiting from optimized automatic-differentiation calculations and operation fusion.
- The forward sensitivity benchmarks used Lotka–Volterra, Brusselator, POLLU, and PK/PD models with the integrators specified in Table I.
- Forward-mode AD computes Jacobian-vector products through dual-number seeding, while reverse-mode AD computes vector-Jacobian products through tracked reals.
- The benchmarks included analytical and automatic-differentiation Jacobians, Jacobian-vector and vector-Jacobian products, and numerical differentiation.
- DSAAD was fastest or nearly fastest across all tested sensitivity-analysis cases.
- Compiler-level operation fusion improves common-subexpression elimination for forward-mode AD relative to implementations that call separate calculation functions.
B. Adjoint Sensitivity Performance Comparisons
Adjoint sensitivity benchmarks compare discrete adjoints with continuous adjoints across the tested models, while the Brusselator results provide a scaling-oriented comparison.
- The adjoint benchmarks report timings for Lotka–Volterra, Brusselator, POLLU, and PK/PD using the listed stiff and non-stiff integrators.
- The Brusselator scaling benchmark varies dimension N and plots parameter count 4N2 against runtime, both on logarithmic axes.
- Forward-mode DSAAD has a clear performance advantage over the other choices on sufficiently small models.
C. Adjoint Sensitivity Scaling
Continuous adjoint methods become more efficient as model size increases, but their performance depends on the adjoint formulation, stiffness, and reverse-mode AD architecture.
- CASA becomes more efficient than DSAAD as Brusselator size increases, with crossover points ranging from about 50 to 150 ODEs+parameters by configuration.
- Around 100 combined ODEs+parameters marks a practical range for considering a shift from forward sensitivity approaches to adjoint methods.
- Quadrature-based CASA is the most efficient CASA method, whereas interpolating adjoints incur larger stiff-solver systems and cubic linear-solve scaling.
- Backsolve adjoints are unstable on stiff equations, and reverse-mode DSAAD scales poorly because scalar tape tracking enlarges computational graphs and memory demands.
- The reverse-mode AD architecture used for internal v′J calculations strongly affects adjoint performance.
D. Parameter Estimation Performance Comparisons
Parameter-estimation benchmarks evaluate whether sensitivity-analysis timing differences transfer to optimization applications using generated data and common convergence settings.
- Parameter estimation used generated data, an L2 loss, evenly spaced observations, and the BFGS local optimizer.
- The parameter-estimation timings are reported in Table III for the tested models.
- All methods started from the same perturbed initial condition and optimization continued until convergence to tolerance 10^-6.
- More efficient sensitivity calculations provided a performance advantage in parameter estimation, while numerical differentiation incurred disadvantages associated with increased gradient error and more optimizer iterations.
IV. DSAAD GENERALIZES TO HYBRID, DELAY, AND DIFFERENTIAL-ALGEBRAIC DIFFERENTIAL EQUATIONS
DSAAD directly generalizes to hybrid, differential-algebraic, and delay differential equations, whereas standard continuous sensitivity analysis requires specialized treatment for these systems. In hybrid systems, unmodified continuous sensitivity equations can produce incorrect derivatives when event timing or jump size depends on parameters.
- Hybrid differential equations: Automatic differentiation propagates sensitivity through event-handling code, while continuous approaches can require special handling for correctness.The discrete approach propagates Dual numbers through user-defined event code.
- Hybrid differential equations: Parameter-dependent event timing and discontinuity size are two omissions that prevent standard continuous sensitivity analysis from correctly propagating sensitivities.The event time can depend on parameters, and the discontinuity amount can also be parameter-dependent.
- Hybrid differential equations: For the analyzed hybrid ODE control problem, standard continuous sensitivity analysis produced incorrect derivatives relative to the analytical solution.The comparison evaluated sensitivities at t = 1 against derivatives obtained directly from the analytical solution.
- Delay differential equations: Continuous sensitivity equations defined for ordinary differential equations do not directly apply to delay differential equations, which require a different equation set.This limitation concerns delay differential equations, including systems with fixed time delays.
- Differential-algebraic and delay differential equations: DSAAD automatically generates sensitivity propagation at the compiler level and thereby supports differential-algebraic, delay, and hybrid differential equation systems without ODE-specific specialization.The approach operates through atomic operations inside the numerical integration scheme and supports expanded systems such as state-dependent-delay DDEs and hybrid DAEs.
V. DISCUSSION
The discussion finds that DSAAD is advantageous for small systems, while continuous adjoints can scale better for sufficiently large problems, with memory, AD-mode, and software-integration trade-offs shaping method choice.
- Performance trade-offs: Forward-mode DSAAD has a strong performance advantage on sufficiently small ODE systems, while continuous adjoints are advantageous on sufficiently large systems.The paper also reports that adjoint runtime overhead can outweigh theoretical scaling benefits on smaller problems.
- Performance trade-offs: Automatic differentiation for Jacobian-vector products improves performance on larger problems, making sensitivity-analysis integration with AD tools important for efficiency.The choice of reverse-mode AD for v′J is also a major performance factor.
- Scope and generalization: The reported package-level comparisons may not generalize across solver ecosystems, motivating comparisons with native SUNDIALS and PETSc TS adjoint techniques.The authors identify this as future research rather than an established result.
- Memory trade-offs: Adjoint implementations trade runtime for memory: quadrature CASA can be fast but stores the full forward solution and interpolant, whereas checkpointing re-solves segments.For the stated ninth-order Runge–Kutta setup, the quadrature memory cost is 26NM.
- Implementation considerations: DSAAD offers implementation advantages for hybrid, delay, and differential-algebraic integrators, but requires AD-compatible solver code and can be difficult to develop elsewhere.The paper notes that applying AD directly is difficult for ecosystems whose solvers call external Fortran implementations.
A. Models
The benchmarks cover non-stiff and stiff biological, pharmacometric, reaction-diffusion, and air-pollution models, including a spatially discretized PDE with many parameters.
- Test models: The Lotka–Volterra model is the first test problem and represents a non-stiff biological system.
- Test models: The Brusselator is a stiff two-dimensional reaction-diffusion PDE discretized into N × N × 2 ODEs.Its spatially dependent parameters produce 4N^2 parameters after discretization.
- Test models: POLLU is a stiff nonlinear air-pollution model consisting of 20 ODEs.
- Test models: The pharmacometric equations include inter-compartmental transfer and concentration-dependent terms for the central and peripheral compartments.
- Test models: The pharmacometric model includes depot dosing, central and peripheral compartments, drug elimination, and a pharmacodynamic response compartment.Additional doses of 100.0 are applied to the depot variable every 24 time units.