Lean Adjoint ODE Methods

• Lean adjoint ODE is a set of techniques that avoid storing full computation graphs by using paired forward and backward ODE solves.
• It enables efficient gradient computation across domains like meta-learning, neural ODE training, and coupled PDE–ODE simulations with both discrete and continuous adjoint formulations.
• These methods balance memory efficiency and accuracy through strategies like checkpointing, recomputation, and exact discrete differentiation.

Lean adjoint ODE denotes a family of adjoint-state constructions for ODE-constrained optimization in which gradient computation is organized to avoid storing or replaying large computational graphs, while retaining either the exact discrete gradient of the implemented solver or a practically accurate continuous-time sensitivity. Across meta-learning, neural ODE training, reduced-order modeling, likelihood-based inference, stochastic optimal control, and coupled PDE–ODE simulation, the shared pattern is to replace long backpropagation chains by a forward solve and a backward adjoint solve, by a scheme-matched discrete adjoint, or, in specially structured settings, by a closed-form adjoint formula (Li et al., 2021, Zhang et al., 2022, Liu et al., 12 Jan 2026, Shin et al., 12 May 2026).

1. Scope and defining idea

In the cited literature, “lean” refers to minimizing memory while preserving gradient correctness or practical gradient fidelity. In neural ODEs, naive backpropagation must retain the entire computation graph of repeated network evaluations, so memory scales with the number of uses times the network size, $O(N\cdot s\cdot |f|)$, which is often prohibitive for accurate solvers that require small step sizes and high order (Matsubara et al., 2021). In MAML-style bilevel optimization, backpropagating through many inner-loop gradient steps likewise requires building and retaining a large computational graph; as steps increase, memory usage explodes and runtime grows dramatically, so practical implementations often restrict the inner loop to only a few steps (Li et al., 2021).

The same pressure appears outside deep learning. In reduced-order modeling, likelihood-based ODE inference, and PDE–ODE coupling, direct sensitivity methods or differentiation through full solver traces can make the computational cost grow with the number of parameters, the number of stored stages, or the number of observations. Adjoint-state methods reduce that dependence because each iteration can be organized around one forward solve and one backward solve, followed by gradient assembly (Liu et al., 12 Jan 2026, Melicher et al., 2016, Farrell et al., 2017).

This suggests that lean adjoint ODE is less a single algorithm than a design principle. The principle is consistent across continuous adjoints, discrete adjoints, symplectic or partitioned Runge–Kutta adjoints, checkpointed reverse sweeps, solver-level algorithmic differentiation, and adjoint-matching constructions in stochastic control (Matsubara et al., 2021, Matsuda et al., 2020, Zhang et al., 2022, Choi et al., 22 Jun 2025).

2. Continuous adjoint formulations

The canonical continuous-time setup considers

$\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$

with objective

$$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$

The adjoint variable $\lambda$ satisfies

$$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$

and the parameter gradient is

$$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$

This formulation is stated in neural ODE training, in operator-inference ROM training, and in rigorous analyses of adjoint derivations (Zhang et al., 2022, Liu et al., 12 Jan 2026, Hu, 2024).

A prominent continuous-time reinterpretation appears in meta-learning. Adjoint MAML treats inner-loop gradient descent as the forward Euler discretization of gradient flow,

$$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$

For the meta-objective $J(x(T))=L_{\mathrm{val}}(x(T);\tau)$, the adjoint $a(t)=\partial J/\partial x(t)$ satisfies

$$a(T)=\frac{\partial L_{\mathrm{val}}}{\partial x}(x(T);\tau), \qquad \dot{a}(t)=-\left[\frac{\partial f}{\partial x}(x(t),t)\right]^\top a(t),$$

and in pure initialization meta-learning one has $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$0 (Li et al., 2021). For gradient flow, $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$1, so the adjoint dynamics are driven by Hessian-vector products $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$2, which can be computed without forming $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$3 explicitly by Pearlmutter’s trick (Li et al., 2021).

Continuous adjoints also accommodate losses that are not purely terminal. In ROM training, the loss is an integral trajectory misfit,

$\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$4

or its weighted variant, and the adjoint has terminal condition $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$5 together with the source term $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$6 (Liu et al., 12 Jan 2026). In statistical ODE models with observations only at discrete times, the continuous adjoint evolves between measurements as $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$7 and acquires jump conditions at observation times. For Gaussian noise,

$\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$8

which is how discrete data are coupled exactly to the continuous adjoint calculus (Melicher et al., 2016).

3. Memory-lean workflows and reconstruction strategies

The operational attraction of continuous adjoints is that sensitivity computation is reduced to a forward evolution and a backward evolution. In A-MAML, one solves the inner ODE forward from $\dot{x}(t)=f(x(t),t,\theta), \qquad x(t_0)=x_0,$9 to $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$0, initializes the adjoint with $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$1, and then integrates the adjoint backward to recover $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$2. No full computational graph or long backpropagation through gradient-descent steps is needed; the method reduces sensitivity computation to two ODE solves, each using only state and adjoint vectors of dimension $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$3, the number of parameters (Li et al., 2021).

The same pattern is explicit in data-driven ROM training. Each optimization iteration requires one forward reduced-order solve, one backward adjoint solve, and gradient assembly through the time integral of $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$4. The paper emphasizes that this cost is approximately two ODE solves per iteration, independent of the number of parameters $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$5, which is lean compared to primal sensitivity methods whose cost scales with $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$6 (Liu et al., 12 Jan 2026).

A central implementation issue is how the backward solve accesses the forward trajectory. Three strategies recur. The first is full storage of states or stages, which maximizes backward convenience but defeats the lean objective for long horizons. The second is recomputation: re-solve forward segments on demand during the backward sweep, which minimizes memory but increases compute. The third is checkpointing, which stores a modest number of states and reconstructs the rest. A-MAML explicitly recommends checkpointing a modest number of states $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$7 and uses a second-order modified Euler backward scheme for robustness; empirically, storing $$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$8–$$L(\theta)=\Phi(x(t_1))+\int_{t_0}^{t_1} q(x(t),t)\,dt.$$9 checkpoints often suffices and beyond that accuracy gains are minor (Li et al., 2021). PNODE uses revolve-style binomial checkpointing tailored to multistage integrators so that memory is proportional to the number of checkpoints rather than to the number of time steps (Zhang et al., 2022).

A further refinement is to change the backward solver’s error control rather than the adjoint equations themselves. In the seminorm method for neural differential equations, the augmented adjoint state is ordered as $\lambda$0, but only $\lambda$1 and $\lambda$2 affect the adjoint dynamics. The proposed seminorm therefore uses

$\lambda$3

dropping $\lambda$4, $\lambda$5, and analogous pure-integral channels such as CNF log-density from the adaptive-step acceptance test. This produces a median improvement of $\lambda$6 fewer function evaluations, and on some problems as much as $\lambda$7 fewer function evaluations (Kidger et al., 2020).

4. Discrete adjoints, reverse accuracy, and exact discrete gradients

A major development in lean adjoint ODE is the distinction between differentiating the continuous ODE and differentiating the discrete solver that was actually executed. PNODE states that the continuous adjoint is not reverse-accurate because training uses a discrete time integrator with specific step sizes, adaptivity rules, and tolerances, whereas the continuous adjoint differentiates the continuous dynamics instead of the exact discrete computation (Zhang et al., 2022). A rigorous note on neural ODE adjoints sharpens the same point: the loss gradient with respect to parameters is an integral, not an ODE, and the adjoint form gives the same results as backpropagation if and only if the discrete adjoint has the same scheme as the discrete neural ODE (Hu, 2024).

High-level discrete adjoint differentiation addresses that mismatch by differentiating the time-integration algorithm as a whole. For a one-step discrete solver $\lambda$8, the discrete adjoint recursion is

$\lambda$9

run backward from $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$0 to $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$1. Because this is the transpose of the actual solver map, it is reverse-accurate by construction (Zhang et al., 2022).

Several works show how to make that exactness compatible with low memory. The symplectic adjoint method pairs a forward Runge–Kutta discretization with a partitioned Runge–Kutta adjoint chosen to satisfy a symplecticity condition. Under that pairing, the discrete bilinear invariant is preserved, and the discrete adjoint equals reverse-mode AD through the RK solver exactly, up to rounding error, while memory scales as $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$2 rather than $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$3 (Matsubara et al., 2021). The generalization to partitioned Runge–Kutta methods proves that, when the forward IVP is solved by a PRK method, an appropriate generalized PRK adjoint gives the exact gradient of the discrete map; when the partition weights coincide, the generalized construction reduces to an ordinary PRK adjoint (Matsuda et al., 2020).

The same program appears in systems software. A C++ library for adaptive explicit Runge–Kutta methods implements the discrete adjoint of the solver program itself, stores only accepted states and times, recomputes stages on demand, and uses automatic adjoint differentiation only for local vector–Jacobian products at stage points. The resulting memory is $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$4 for stored states and times rather than $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$5 for storing all stage data, and SIMD vectorization is used to compute multiple vector–Jacobian products concurrently (Martins et al., 2024).

This body of work shifts the meaning of “lean” from merely reducing memory to enforcing gradient consistency at the solver level. In that sense, lean adjoint ODE includes both continuous forward–backward sensitivity systems and exact discrete reverse sweeps, provided the reduction in memory does not sever the relation between the computed gradient and the actual numerical method (Matsubara et al., 2021, Zhang et al., 2022, Hu, 2024).

5. Extensions to nonsmooth, constrained, and coupled systems

Lean adjoint constructions are not confined to smooth unconstrained ODEs. For vector fields that are Lipschitz and path differentiable, flows inherit path differentiability, and the forward sensitivity differential inclusion provides a conservative Jacobian for the flow. This yields a nonsmooth adjoint method for ODE-constrained integral costs: $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$6 with parameter-gradient element

$$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$7

The method is described as lean because the backward object is a single $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$8-dimensional ODE or differential inclusion rather than the full $$\frac{d\lambda}{dt}=-\left(\frac{\partial f}{\partial x}(x(t),t,\theta)\right)^\top \lambda-\frac{\partial q}{\partial x}(x(t),t), \qquad \lambda(t_1)=\frac{\partial \Phi}{\partial x}(x(t_1)),$$9 sensitivity system (Marx et al., 2022).

In coupled PDE–ODE simulation, lean adjoints arise through high-level generation of discrete adjoint solvers rather than manual differentiation of low-level kernels. The FEniCS extension for spatially distributed ODEs and PDE–ODE coupling augments UFL with a vertex measure, allows arbitrary Butcher tableaux for ODE discretization, and automatically derives adjoint and tangent-linear solvers for the coupled system. The recommended lean strategy is to recompute intermediate Runge–Kutta stage states during the backward sweep instead of storing them; the observed adjoint-to-forward runtime ratio for ODE solves is approximately $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$0, while memory remains low because stages are not retained (Farrell et al., 2017).

Structure-preserving dynamics require a different discrete calculus. For variational integrators of unconstrained and holonomically constrained ODEs, the discrete adjoint is derived directly from the discrete Lagrange–d’Alembert principle. In constrained systems, null-space projection $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$1 and nodal reparametrization $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$2 eliminate multipliers and yield backward recursions in minimal coordinates. The backward sweep uses the stored forward trajectory together with local derivatives of the discrete Lagrangian and discrete forces; the method is structurally matched to the forward variational integrator rather than to the continuous Euler–Lagrange equations (Schubert et al., 2023).

These extensions show that lean adjoint ODE is compatible with nonsmooth analysis, geometric integration, and distributed multiphysics, provided the adjoint is formulated in the same analytic or discrete category as the forward model (Marx et al., 2022, Farrell et al., 2017, Schubert et al., 2023).

6. Application domains, empirical behavior, and limitations

Meta-learning provides one of the clearest demonstrations of the memory–accuracy trade-off. In A-MAML, full MAML is accurate for short trajectories but memory explodes with the number of inner steps, whereas A-MAML is accurate along the actual trajectory without extra regularizers. In experiments, $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$3 inner steps exhaust GPU memory for MAML, while A-MAML supports $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$4 steps with $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$5 GPU usage; on Omniglot $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$6-way $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$7-shot, A-MAML achieves $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$8, and on Mini-ImageNet $$\frac{dL}{d\theta}=\frac{\partial \Phi}{\partial \theta}+\int_{t_0}^{t_1}\left(\lambda(t)^\top \frac{\partial f}{\partial \theta}(x(t),t,\theta)+\frac{\partial q}{\partial \theta}\right)dt.$$9-way $$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$0-shot it reaches about $$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$1 (Li et al., 2021).

In reduced-order modeling, the adjoint-based training framework is comparable to standard operator inference on clean data, but under sparse temporal sampling and additive Gaussian noise it delivers lower test relative state error and better roll-out stability. The robustness comes from trajectory fitting in continuous time, which avoids explicit numerical differentiation of noisy data and introduces temporal regularization through the integral loss (Liu et al., 12 Jan 2026).

For neural ODE training, the discrete-adjoint viewpoint changes both performance and solver choice. PNODE is reported to achieve the highest memory efficiency among reverse-accurate methods, to be up to two times faster than the vanilla neural ODE, and to be up to $$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$2 times faster than the best existing reverse-accurate method on image classification problems; it also enables implicit time integration needed for stiff dynamical systems (Zhang et al., 2022). The symplectic adjoint method achieves exact gradients up to rounding error, consumes much less memory than naive backpropagation and checkpointing schemes, and is more robust to rounding errors than the classical adjoint (Matsubara et al., 2021). The seminorm method reduces backward function evaluations with only minor code modifications, leaving the forward pass unchanged (Kidger et al., 2020).

A more radical notion of leanness appears in stochastic optimal control for diffusion and flow models. Efficient Adjoint Matching reformulates reward fine-tuning with a linear base drift and a modified terminal cost so that the adjoint admits a closed-form solution

$$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$3

thereby eliminating backward adjoint simulation altogether. On standard text-to-image reward fine-tuning benchmarks, EAM converges up to $$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$4 faster than AM and matches or surpasses it across PickScore, ImageReward, HPSv2.1, CLIPScore, and Aesthetics (Shin et al., 12 May 2026). NAAS uses an annealed reference process and a lean adjoint system inspired by adjoint matching; it avoids importance weighting and reports improvements of more than $$\dot{x}(t)=f(x(t),t)=-\nabla_x L_{\mathrm{train}}(x(t);\tau), \qquad x(0)=\theta.$$5 over the second-best method on MW54 and GMM40 in Sinkhorn and MMD metrics (Choi et al., 22 Jun 2025).

The limitations are equally consistent across the literature. Continuous adjoints can suffer from discretization mismatch, especially with large step sizes, adaptive steppers, or momentum-based optimizers (Li et al., 2021, Hu, 2024). Very long trajectories can cause backward divergence or stiff dynamics, so higher-order forward solvers, checkpointing, specialized stiff solvers, or smaller step sizes may be required (Li et al., 2021, Zhang et al., 2022). Hessian-vector products can amplify minibatch noise in meta-learning and stochastic optimization (Li et al., 2021). In statistical ODE inference, adjoint runtime increases with the number of observation times because each observation introduces an adjoint jump (Melicher et al., 2016). For nonsmooth systems, measurable conservative-Jacobian selections must be available, and numerical schemes should remain consistent near kinks (Marx et al., 2022). For constrained variational problems and large PDE–ODE systems, the backward sweep requires either storing substantial space–time state information or accepting additional recomputation (Farrell et al., 2017, Schubert et al., 2023).

Taken together, these results define lean adjoint ODE as a technically heterogeneous but conceptually unified area: exact or approximate sensitivity propagation is retained, while storage of long computation graphs, stage tapes, or pathwise parameter adjoints is reduced by solver-level design, checkpointing, recomputation, structure preservation, or closed-form adjoint elimination (Matsubara et al., 2021, Zhang et al., 2022, Shin et al., 12 May 2026).