Geometric Methods for Adjoint Systems

Adjoint systems are widely used to inform control, optimization, and design in systems described by ordinary differential equations or differential-algebraic equations. In this paper, we explore the geometric properties and develop methods for such adjoint systems. In particular, we utilize symplectic and presymplectic geometry to investigate the properties of adjoint systems associated with ordinary differential equations and differential-algebraic equations, respectively. We show that the adjoint variational quadratic conservation laws, which are key to adjoint sensitivity analysis, arise from (pre)symplecticity of such adjoint systems. We discuss various additional geometric properties of adjoint systems, such as symmetries and variational characterizations. For adjoint systems associated with a differential-algebraic equation, we relate the index of the differential-algebraic equation to the presymplectic constraint algorithm of Gotay et al. [18]. As an application of this geometric framework, we discuss how the adjoint variational quadratic conservation laws can be used to compute sensitivities of terminal or running cost functions. Furthermore, we develop structure-preserving numerical methods for such systems using Galerkin Hamiltonian variational integrators (Leok and Zhang [23]) which admit discrete analogues of these quadratic conservation laws. We additionally show that such methods are natural, in the sense that reduction, forming the adjoint system, and discretization all commute, for suitable choices of these processes. We utilize this naturality to derive a variational error analysis result for the presymplectic variational integrator that we use to discretize the adjoint DAE system. Finally, we discuss the application of adjoint systems in the context of optimal control problems, where we prove a similar naturality result.


Introduction
1.1. Applications of the Adjoint Equations. The solution of many nonlinear problems involves successive linearization, and as such variational equations and their adjoints play a critical role in a variety of applications. Adjoint equations are of particular interest when the parameter space is significantly higher dimension than that of the output or objective. In particular, the simulation of adjoint equations arise in sensitivity analysis [10; 11], adaptive mesh refinement [25], uncertainty quantification [39], automatic differentiation [19], superconvergent functional recovery [28], optimal control [31], optimal design [16], optimal estimation [27], and deep learning viewed as an optimal control problem [5].
The study of geometric aspects of adjoint systems arose from the observation that the combination of any system of differential equations and its adjoint equations are described by a formal Lagrangian [20; 21]. This naturally leads to the question of when the formation of adjoints and discretization commutes [36], and prior work on this include the Ross-Fahroo lemma [32], and the observation by Sanz-Serna [33] that the adjoints and discretization commute if and only if the discretization is symplectic.
1.2. Symplectic and Presymplectic Geometry. Throughout the paper, we will assume that all manifolds and maps are smooth, unless otherwise stated. Let (P, Ω) be a (finite-dimensional) symplectic manifold, i.e., Ω is a closed nondegenerate two-form on P . Given a Hamiltonian H : P R, the Hamiltonian system is defined by where the vector field X H is a section of the tangent bundle to P . By nondegeneracy, the vector field X H exists and is uniquely determined. For an open interval I ⊂ R, we say that a curve z : I P is a solution of Hamilton's equations if z is an integral curve of X H , i.e.,ż(t) = X H (z(t)) for all t ∈ I.
A particularly important example for our purposes is when the symplectic manifold is the cotangent bundle of a manifold, P = T * M , equipped with the canonical symplectic form Ω = dq∧dp in natural coordinates (q, p) on T * M . A Hamiltonian system has the coordinate expressioṅ q = ∂H(q, p) ∂p , p = − ∂H(q, p) ∂q .
By Darboux's theorem, any symplectic manifold is locally symplectomorphic to a cotangent bundle equipped with its canonical symplectic form. As such, any Hamiltonian system can be locally expressed in the above form (even when P is not a cotangent bundle), using Darboux coordinates.
We now consider the generalization of Hamiltonian systems where we relax the condition that Ω is nondegenerate, i.e., presymplectic geometry. Let (P, Ω) be a presymplectic manifold, i.e., Ω is a closed two-form on P with constant rank. As before, given a Hamiltonian H : P R, we define the associated Hamiltonian system as i X H Ω = dH.
Note that since Ω is now degenerate, X H is not guaranteed to exist and if it does, it need not be unique and in general is only partially defined on a submanifold of P . Again, we say a curve on P is a solution to Hamilton's equations if it is an integral curve of X H . Using Darboux coordinates (q, p, r) adapted to (P, Ω), where Ω = dq ∧ dp and ker(Ω) = span{∂/∂r}, the local expression for Hamilton's equations is given byq = ∂H(q, p, r) ∂p , p = − ∂H(q, p, r) ∂q , 0 = ∂H(q, p, r) ∂r .
The third equation above is interpreted as a constraint equation which any solution curve must satisfy. We will assume that the constraint defines a submanifold of P . It is clear that in order for a solution vector field X H to exist, it must be restricted to lie on this submanifold. However, in order for its flow to remain on the submanifold, it must be tangent to this submanifold, which further restricts where X can be defined. Alternating restriction in order to satisfy these two constraints yields the presymplectic constraint algorithm of Gotay et al. [18]. The presymplectic constraint algorithm begins with the observation that for any X satisfying the above system, so does X + Z, where Z ∈ ker(Ω). In order to obtain such a vector field X, one considers the subset P 1 of P such that Z p (H) = 0 for any Z ∈ ker(Ω), p ∈ P 1 . We will assume that the set P 1 is a submanifold of P . We refer to P 1 as the primary constraint manifold. In order for the flow of the resulting Hamiltonian vector field X to remain on P 1 , one further requires that X is tangent to P 1 . The set of points satisfying this property defines a subsequent secondary constraint submanifold P 2 .
Iterating this process, one obtains a sequence of submanifolds · · · P k · · · P 1 P 0 ≡ P, defined by (1.1) P k+1 = {p ∈ P k : Z p (H k ) = 0 for all Z ∈ ker(Ω k )}, where Ω k+1 = Ω k | P k+1 , If there exists a nontrivial fixed point in this sequence, i.e., a submanifold P k of P such that P k = P k+1 , we refer to P k as the final constraint manifold. If such a fixed point exists, we denote by ν P the minimum integer such that P ν P = P ν P +1 , i.e., ν P is the number of steps necessary for the presymplectic constraint algorithm to terminate. If such a final constraint manifold P ν P exists, there always exists a solution vector field X defined on and tangent to P ν P such that i X Ω ν P = dH ν P and X is unique up to the kernel of Ω ν P . Furthermore, such a final constraint manifold is maximal in the sense that if there exists a submanifold N of P which admits a vector field X defined on and tangent to N such that i X Ω| N = dH| N , then N ⊂ P ν P (Gotay and Nester [17]).
1.3. Main Contributions. In this paper, we explore the geometric properties of adjoint systems associated with ordinary differential equations (ODEs) and differential-algebraic equations (DAEs). For a discussion of adjoint systems associated with ODEs and DAEs, see Sanz-Serna [33] and Cao et al. [11], respectively. In particular, we utilize the machinery of symplectic and presymplectic geometry as a basis for understanding such systems.
In Section 2.1, we review the notion of adjoint equations associated with ODEs over vector spaces. We show that the quadratic conservation law, which is the key to adjoint sensitivity analysis, arises from the symplecticity of the flow of the adjoint system. In Section 2.2, we investigate the symplectic geometry of adjoint systems associated with ODEs on manifolds. We additionally discuss augmented adjoint systems, which are useful in the adjoint sensitivity of running cost functions. In Section 2.3, we investigate the presymplectic geometry of adjoint systems associated with DAEs on manifolds. We investigate the relation between the index of the base DAE and the index of the associated adjoint system, using the notions of DAE reduction and the presymplectic constraint algorithm. We additionally consider augmented systems for such adjoint DAE systems. For the various adjoint systems that we consider, we derive various quadratic conservation laws which are useful in adjoint sensitivity analysis of terminal and running cost functions. We additionally discuss symmetry properties and present variational characterizations of such systems that provide a useful perspective for constructing geometric numerical methods for these systems.
In Section 3, we discuss applications of the various adjoint systems to adjoint sensitivity and optimal control. In Section 3.1, we show how the quadratic conservation laws developed in Section 2 can be used for adjoint sensitivity analysis of running and terminal cost functions, subject to ODE or DAE constraints. In Section 3.2, we construct structure-preserving discretizations of adjoint systems using the Galerkin Hamiltonian variational integrator construction of Leok and Zhang [23]. For adjoint DAE systems, we introduce a presymplectic analogue of the Galerkin Hamiltonian variational integrator construction. We show that such discretizations admit discrete analogues of the aforementioned quadratic conservation laws and hence are suitable for the numerical computation of adjoint sensitivities. Furthermore, we show that such discretizations are natural when applied to DAE systems, in the sense that reduction, forming the adjoint system, and discretization all commute (for particular choices of these processes). As an application of this naturality, we derive a variational error analysis result for the resulting presymplectic variational integrator for adjoint DAE systems. Finally, in Section 3.3, we discuss adjoint systems in the context of optimal control problems, where we prove a similar naturality result, in that suitable choices of reduction, extremization, and discretization commute.
By developing a geometric theory for adjoint systems, the application areas that utilize such adjoint systems can benefit from the existing work on geometric and structure-preserving methods.
1.4. Main Results. In this paper, we prove that, starting with an index 1 DAE, appopriate choices of reduction, discretization, and forming the adjoint system commute. That is, the following diagram commutes.
In the process of defining the ingredients in the above diagram, we will additionally prove various properties of adjoint systems associated with ODEs and DAEs. The key properties that we will prove for such adjoint systems are the adjoint variational quadratic conservation laws, Propositions 2.3, 2.7, 2.11, 2.12. As we will show, these conservation laws can be used to compute adjoint sensitivities of running and terminal cost functions under the flow of an ODE or DAE. In order to prove these conservation laws, we will need to define the variational equations associated with an adjoint system. We will define them as the linearization of the base ODE or DAE; for the DAE case, we will show that the variational equations have the same index as the base DAE so that they have the same (local) solvability.

Adjoint Systems
2.1. Adjoint Equations on Vector Spaces. In this section, we review the notion of adjoint equations on vector spaces and their properties, as preparation for adjoint systems on manifolds.
Let Q be a finite-dimensional vector space and consider the ordinary differential equation on Q given by where f : Q Q is a differentiable vector field on Q. Let Df (q) denote the linearization of f at q ∈ Q, Df (q) ∈ L(Q, Q). Denoting its adjoint by [Df (q)] * ∈ L(Q * , Q * ), the adjoint equation associated with (2.1) is given by where p is a curve on Q * .
Let q A be coordinates for Q and let p A be the associated dual coordinates for Q * , so that the duality pairing is given by p, q = p A q A . The linearization of f at q is given in coordinates by

Its adjoint then acts on
Thus, the ODE and its adjoint can be expressed in coordinates aṡ Next, we recall that the combined system (2.1)-(2.2), which we refer to as the adjoint system, arises from a variational principle. Letting ·, · denote the duality pairing between Q * and Q, we define the Hamiltonian The associated action, defined on the space of curves on Q × Q * covering some interval (t 0 , t 1 ), is given by Proposition 2.1. The variational principle δS = 0, subject to variations (δq, δp) which fix the endpoints δq(t 0 ) = 0, δq(t 1 ) = 0, yields the adjoint system (2.1)-(2.2).
Proof. Compute the variation of S with respect to a compactly supported variation (δq, δp), and one extremizes over C 2 curves from [t 0 , t 1 ] to T * Q such that q(t 0 ) = q 0 , p(t 1 ) = p 1 . The Type II variational principle again gives the above adjoint system, but with differing boundary conditions. These boundary conditions are typical in adjoint sensitivity analysis, where one fixes the initial position and the final momenta.
The variational principle utilized above is formulated so that the stationarity condition δS = 0 is equivalent to Hamilton's equations, where we view Q × Q * ∼ = T * Q with the canonical symplectic form on the cotangent bundle Ω = dq ∧ dp and with the corresponding Hamiltonian H : T * Q R given as above. It then follows that the flow of the adjoint system is symplectic.
The symplecticity of the adjoint system is a key feature of the system. In fact, the symplecticity of the adjoint system implies that a certain quadratic invariant is preserved along the flow of the system. This quadratic invariant is the key ingredient to the use of adjoint equations for sensitivity analysis. To state the quadratic invariant, consider the variational equation associated with equation (2.1), which corresponds to the linearization of (2.1) at q ∈ Q. For solution curves p and δq to (2.2) and (2.3), respectively, over the same curve q, one has that the quantity p, δq is preserved along the flow of the system, since To see that symplecticity implies the preservation of this quantity, recall that symplecticity is the statement that, along a solution curve of the adjoint system (2.1)-(2.2), one has where V and W are first variations to the adjoint system (i.e., that the flow of V and W on solutions are again solutions). Infinitesimally, first variations V and W correspond to solutions of the linearization of the adjoint system (2.1)-(2.2). At a solution (q, p) to the adjoint system, the linearization of the system is given by Note that the first equation is just the variational equation (2.3) while the second equation is the adjoint equation (2.2), with p replaced by δp, since the adjoint equation is linear in p. The first variation vector field V corresponding to a solution (δq, δp) of this linearized system is Now, we make two choices for the first variations V and W . For W , we take the solution δq = 0, δp = p of the linearized system, which gives W = p ∂/∂p. For V , we take the solution δq = δq, δp = 0 of the linearized system, which gives V = δq ∂/∂q. Inserting these into Ω gives Thus, symplecticity d dt Ω(V, W ) = 0 with this particular choice of first variations V, W gives the preservation of the quadratic invariant p, δq .
2.2. Adjoint Systems on Manifolds. We now extend the notion of the adjoint system to the case where the configuration space of the base ODE is a manifold. We will provide a symplectic characterization of these adjoint systems, prove the associated adjoint variational quadratic conservation laws, and additionally discuss symmetries and variational principles associated with these systems.
Let M be a manifold and consider the ODE on M given by where f is a vector field on M . Letting π : T M M denote the tangent bundle projection, we recall that a vector field f is a map f : M T M which satisfies π • f = 1 M , i.e., f is a section of the tangent bundle.
Analogous to the adjoint system on vector spaces, we will define the adjoint system on a manifold as an ODE on the cotangent bundle T * M which covers (2.4), such that the time evolution of the momenta in the fibers of T * M are given by an adjoint linearization of f .
To do this, in analogy with the vector space case, consider the Hamiltonian H : T * M R given by H(q, p) = p, f (q) q where ·, · q is the duality pairing of T * q M with T q M . When there is no possibility for confusion of the base point, we simply denote this duality pairing as ·, · . Recall that the cotangent bundle T * M possesses a canonical symplectic form Ω = −dΘ where Θ is the tautological one-form on T * M . With coordinates (q, p) = (q A , p A ) on T * M , this symplectic form has the coordinate expression Ω = dq ∧ dp ≡ dq A ∧ dp A .
We define the adjoint system as the ODE on T * M given by Hamilton's equations, with the above choice of Hamiltonian H and the canonical symplectic form. Thus, the adjoint system is given by the equation i X H Ω = dH, whose solution curves on T * M are the integral curves of the Hamiltonian vector field X H . As is well-known, for the particular choice of Hamiltonian H(q, p) = p, f (q) , the Hamiltonian vector field X H is given by the cotangent lift f of f , which is a vector field on T * M that covers f (see, for example, Bullo and Lewis [8]). With coordinates z = (q, p) on T * M , the adjoint system is the ODE on T * M given by To be more explicit, recall that the cotangent lift of f is constructed as follows. Let Φ ǫ : M M denote the one-parameter family of diffeomorphisms generated by f . Then, we consider the cotangent lifted diffeomorphisms given by (Φ −ǫ ) * : T * M T * M . This covers Φ ǫ in the sense that The cotangent lift f is then defined to be the infinitesimal generator of the cotangent lifted flow, We can directly verify that f is the Hamiltonian vector field for H, which follows from where Lf Θ = 0 follows from the fact that cotangent lifted flows preserve the tautological one-form and H = i f Θ follows from a direct computation (where i f Θ is interpreted as a function on the cotangent bundle which maps (q, p) to Θ(q, p), f (q, p) ) The adjoint system (2.5) covers (2.4) in the following sense.
Proposition 2.2. Integral curves to the adjoint system (2.5) lift integral curves to the system (2.4).
Proof. Let z = (q, p) be coordinates on T * M . Let (q,ṗ) ∈ T (q,p) T * M . Then, T π T * M (q,ṗ) =q where T π T * M is the pushforward of the cotangent projection. Furthermore, Thus, the pushforward of the cotangent projection applied to (2.5) gives (2.4). It then follows that integral curves of (2.5) lift integral curves of (2.4).
Remark 2.2. This can also be seen explicitly in coordinates.
Recalling that i f Ω = dH, one has and, on the other hand, denoting f (q, Equating these two gives the coordinate expression for the cotangent lift f , Thus, the systemż = f (z) can be expressed in coordinates aṡ which clearly covers the original ODEq A = f A (q). Also, note that this coordinate expression for the adjoint system recovers the coordinate expression for the adjoint system in the vector space case.
Analogous to the vector space case, the adjoint system possesses a quadratic invariant associated with the variational equations of (2.4). The variational equation is given by considering the tangent lifted vector field on T M , f : T M T T M , which is defined in terms of the flow Φ ǫ generated by f by where (q, δq) are coordinates on T M . That is, f is the infinitesimal generator of the tangent lifted flow. The variational equation associated with (2.4) is the ODE associated with the tangent lifted vector field. In coordinates, Proposition 2.3. For integral curves (q, p) of (2.5) and (q, δq) of (2.7), which cover the same curve q, Proof. Note that (q(t), p(t)) ∈ T * q(t) M and (q(t), δq(t)) ∈ T q(t) M so the duality pairing is welldefined. Then, , so the pairing is constant.

Remark 2.3.
In the vector space case, we saw that the preservation of the quadratic invariant is implied by symplecticity. The above result is analogously implied by symplecticity, noting that the flow of the adjoint system is symplectic since f is a Hamiltonian vector field.
Another conserved quantity for the adjoint system (2.5) is the Hamiltonian, since the adjoint system corresponds to a time-independent Hamiltonian flow, d dt H = Ω(X H , X H ) = 0. Additionally, conserved quantities for adjoint systems are generated, via cotangent lift, by symmetries of the original ODE (2.4), where we say that a vector field g is a symmetry of the ODĖ Proposition 2.4. Let g be a symmetry of (2.4), i.e., [g, f ] = 0. Then, its cotangent lift g is a symmetry of (2.5) and additionally, the function Θ, g on T * M is preserved along the flow of f , i.e., under the flow of the adjoint system (2.5).
Proof. We first show that g is a symmetry of (2.5), i.e., that [ g, f ] = 0. To see this, we recall that the cotangent lift of the Lie bracket of two vector fields equals the Lie bracket of their cotangent lifts, To see that Θ, g is preserved along the flow of f , we have where we used that L f Θ = 0 since f is a cotangent lifted vector field.
Remark 2.4. The above proposition states when [f, g] = 0, the Hamiltonian for the adjoint system associated with g, Θ, g , is preserved along the Hamiltonian flow corresponding to the Hamiltonian for the adjoint system associated with f , Θ, f , and vice versa. Note, Θ, g can be interpreted as the momentum map corresponding to the action on T * M given by the flow of g.
The above proposition shows that (at least some) symmetries of the adjoint system (2.5) can be found by cotangent lifting symmetries of the original ODE (2.4). Additionally, the above proposition states that such cotangent lifted symmetries give rise to conserved quantities.
In light of the above proposition, it is natural to ask the following question. Given a symmetry G of the adjoint system (2.5) (i.e., [G, f ] = 0), does it arise from a cotangent lifted symmetry in the sense of Proposition 2.4? In general, the answer is no. However, for a projectable vector field G which is a symmetry of the adjoint system, its projection by T π T * M to a vector field on M does satisfy the assumptions of Proposition 2.4. This gives the following partial converse to the above proposition.
Proposition 2.5. Let G be a projectable vector field on the bundle π T * M : T * M M which is a symmetry of (2.5), i.e., [G, f ] = 0. Then, the pushforward vector field g = T π T * M (G) on M satisfies the assumptions of Proposition 2.4 and T π T * M g = T π T * M G.
Proof. Since G is a projectable vector field on the cotangent bundle, g = T π T * M G defines a welldefined vector field on M . Thus, so g is a symmetry of (2.4). Furthermore, we also have The preceding proposition shows that, for the class of projectable symmetries of the adjoint system (2.5), it is always possible to find an associated symmetry of the original ODE (2.4) which, by Proposition 2.4, corresponds to a Hamiltonian symmetry. Note that this implies that we can associate a conserved quantity Θ, g to G, where g = T π T * M G. Furthermore, since T π T * M g = T π T * M G and the canonical form Θ is a horizontal one-form, this implies that Θ, G equals Θ, g and hence, is conserved.
These two propositions show that symmetries of an ODE can be identified with equivalence classes of projectable symmetries of the associated adjoint system, where two projectable symmetries are equivalent if their difference lies in the kernel of T π T * M .
We also recall that the adjoint system (2.5) formally arises from a variational principle. To do so, let Θ be the tautological one-form on T * M . The action is defined to be where ψ(t) = (q(t), p(t)) is a curve on T * M over the interval I = (t 0 , t 1 ). We consider the variational principle δS[ψ] = 0, subject to variations which fix the endpoints q(t 0 ), q(t 1 ).
Proposition 2.6. Let ψ be a curve on T * M over the interval I. Then, ψ is a stationary point of S with respect to variations which fix q(t 0 ), q(t 1 ) if and only if (2.5) holds.
The proof of the above proposition is standard in the literature, so we will omit it.
Remark 2.5. In coordinates, the above action (2.9) takes the form which is the same coordinate expression as the action in the vector space case.

Adjoint Systems with Augmented Hamiltonians.
In this section, we consider a class of modified adjoint systems, where some function on the base manifold M is added to the Hamiltonian of the adjoint system. More precisely, let H : T * M R, H(q, p) = p, f (q) be the Hamiltonian of the previous section, corresponding to the ODEq = f (q). Let L : M R be a function on M . We identify L with its pullback through π T * M : T * M M . Then, we define the augmented Hamiltonian We define the augmented adjoint system as the Hamiltonian system associated with H L relative to the canonical symplectic form Ω on T * M , Remark 2.6. The motivation for such systems arises from adjoint sensitivity analysis and optimal control. For adjoint sensitivity analysis of a running cost function, one is concerned with the sensitivity of some functional t 0 L(q)dt along the flow of the ODEq = f (q). In the setting of optimal control, the goal is to minimize such a functional, constrained to curves satisfying the ODE (see, for example, Aguiar et al. [2]). We will discuss such applications in more detail in Section 3.
In coordinates, the augmented adjoint system (2.10) takes the forṁ We now prove various properties of the augmented adjoint system, analogous to the previous section. To start, first note that we can decompose the Hamiltonian vector field X H L as follows. Let f be the cotangent lift of f . Let X L ≡ X H L − f . Then, observe that Thus, we have the decomposition X H L = f + X L , where f and X L are the Hamiltonian vector fields for H and L, respectively. In coordinates, From the coordinate expression, we see that X L is a vertical vector field over the bundle T * M M . We can also see this intrinsically, since dL is a horizontal one-form on T * M , X L satisfies i X L Ω = dL, and Ω restricts to an isomorphism from vertical vector fields on T * M to horizontal one-forms on T * M . Thus, it is immediate to see intrinsically that an analogous statement to Proposition 2.2 holds, since the flow of f lifts the flow of f , while the flow of X L is purely vertical. That is, since T π T * M X L = 0, T π T * M X H L = T π T * M f = f. We can of course also see that the augmented adjoint system lifts the original ODE from the coordinate expression for the augmented adjoint system, (2.11a)-(2.11b).
We now prove analogous statements to Propositions 2.3 and 2.4, modified appropriately for the presence of L in the augmented Hamiltonian. Proposition 2.7. Let (q, p) be an integral curve of the augmented adjoint system (2.10) and let (q, δq) be an integral curve of the variational equation (2.7), covering the same curve q. Then, d dt p, δq = − dL, δq .
Remark 2.7. Note that the variational equation associated with the above system is the same as in the nonaugmented case, equation (2.7), since augmenting L to the Hamiltonian system only shifts the Hamiltonian vector field in the vertical direction.
Proof. We will prove this in coordinates. We have the equationṡ Then, Remark 2.8. Interestingly, the above proposition states that in the augmented case, p, δq is no longer preserved but rather, its change measures the change of L with respect to the variation δq. This may at first seem contradictory since both the augmented and nonaugmented Hamiltonian vector fields, X H L and X H , preserve Ω, and as we noted previously in Remark 2.3, the preservation of the quadratic invariant is implied by symplecticity. However, upon closer inspection, there is no contradiction because the two cases have different first variations, where recall a first variation is a symmetry vector field of the Hamiltonian system and symplecticity can be stated as for first variation vector fields V and W . In the nonaugmented case, the equations satisfied by the first variation of the momenta p can be identified with p itself, since the adjoint equation for p is linear in p. On the other hand, in the augmented case, the adjoint equation for p, (2.11b), is no longer linear in p, rather, it is affine in p. Furthermore, the failure of this equation to be linear in p is given precisely by −dL. Thus, in the augmented case, first variations in p can no longer be identified with p, and this leads to the additional term − dL, δq in the above proposition.
To prove an analogous statement to Proposition 2.4, we need the additional assumption that the symmetry vector field g leaves L invariant, L g L = 0.
Proposition 2.8. Let g be a symmetry of the ODEq = f (q), i.e., [g, f ] = 0. Additionally, assume that g is a symmetry of L, i.e., L g L = 0. Then, its cotangent lift g is a symmetry of the augmented adjoint system, [ g, X H L ] = 0 and additionally, the function Θ, g on T * M is preserved along the flow of X H L .
Proof. To see that [ g, X H L ] = 0, note that with the decomposition X H L = f + X L , we have where we used that [ g, f ] = [g, f ] = 0. To see that [ g, X L ] = 0, we note that [ g, X L ] can be expressed where we interpret Ω : T (T * M ) T * (T * M ). Then, note that g preserves Ω since g is a cotangent lift and it also preserves L (where, since we identify L with its pullback through π T * M , this is equivalent to g preserving L). More precisely, since we are identifying L with its pullback (π T * M ) * L, Hence, L g (Ω −1 (dL)) = 0. One can also verify this in coordinates, and a direct computation yields where we used that L f Θ, g = 0 by Proposition 2.4. Now, we have The first term above vanishes since L g L = 0. Furthermore, d(i X L Θ), g = 0 since X L is a vertical vector field while Θ is a horizontal one-form. Hence, L X H L Θ, g = 0.

Adjoint Systems for DAEs via Presymplectic Mechanics.
In this section, we generalize the notion of adjoint system to the case where the base equation is a (semi-explicit) DAE. We will prove analogous results to the ODE case. However, more care is needed than the ODE case, since the DAE constraint introduces issues with solvability. As we will see, the adjoint system associated with a DAE is a presymplectic system, so we will approach the solvability of such systems through the presymplectic constraint algorithm.
We consider the following setup for a differential-algebraic equation. Let M d and M a be two manifolds, where we regard M d as the configuration space of the "dynamical" or "differential" variables and M a as the configuration space of the "algebraic" variables. Let M d be the projection onto the first factor and let π T M d : Then, a (semi-explicit) DAE is specified by a section f ∈ Γ(T M d ) and a section φ ∈ Γ(Φ), via the systeṁ where (q, u) are coordinates on M d × M a . We refer to T M d as the differential tangent bundle, with coordinates (q, u, v) and to Φ as the constraint bundle. Remark 2.9. For the local solvability of (2.12a)-(2.12b), regard φ locally as a map R rank(Φ) . If ∂φ/∂u is an isomorphism at a point (q 0 , u 0 ) where Φ(q 0 , u 0 ) = 0, then by the implicit function theorem, one can locally solve u = u(q) about (q 0 , u 0 ) such that φ(q, u(q)) = 0, and subsequently solve the unconstrained differential equationq = f (q, u(q)) locally. This is the case for semi-explicit index 1 DAEs.
In order for the rank(Φ) × dim(M a ) matrix ∂φ/∂u(q 0 , u 0 ) to be an isomorphism, it is necessary that rank(Φ) = dim(M a ). However, we will make no such assumption, so as to treat the theory in full generality, allowing for, e.g., nonunique solutions. Now, let T * M d be the pullback bundle of the cotangent bundle T * M d by π d , with coordinates (q, u, p), which we refer to as the differential cotangent bundle. Furthermore, let Φ * be the dual vector bundle to Φ, with coordinates (q, u, λ). Let T * M d ⊕ Φ * be the Whitney sum of these two vector bundles over M d ×M a with coordinates (q, u, p, λ), which we refer to as the generalized phase space bundle. We define a Hamiltonian on the generalized phase space, Let Ω d denote the canonical symplectic form on T * M d , with coordinate expression Ω d = dq ∧ dp.
We define a presymplectic form Ω 0 on T * M d ⊕ Φ * as follows: the pullback bundle admits the mapπ d : T * M d T * M d which covers π d and acts as the identity on fibers and furthermore, the generalized phase space bundle admits the projection Π : the Whitney sum has the structure of a double vector bundle. Hence, we can pullback Ω d along the sequence of maps on the generalized phase space bundle. Clearly, Ω 0 is closed as the pullback of a closed form. In general, Ω 0 will be degenerate except in the trivial case where M a is empty and the fibers of Φ are the zero vector space. Hence, Ω 0 is a presymplectic form. Note that since Π acts by projection andπ d acts as the identity on fibers, the coordinate expression for Ω 0 on T * M d ⊕Φ * with coordinates (q, u, p, λ) is the same as the coordinate expression for Ω d , Ω 0 = dq ∧ dp. The various spaces and their coordinates are summarized in the diagram below.
We now define the adjoint system associated with the DAE (2.12a)-(2.12b) as the Hamiltonian system Given a (generally, partially defined) vector field X on the generalized phase space satisfying (2.13), we say a curve (q(t), u(t), p(t), λ(t)) is a solution curve of (2.13) if it is an integral curve of X.
Let us find a coordinate expression for the above system. Expressing our coordinates with indices (q i , u a , p j , λ A ), the left hand side of (2.13) along a solution curve has the expression On the other hand, the right hand side of (2.13) has the expression Equating these expressions gives the coordinate expression for the adjoint DAE system, Remark 2.10. As mentioned in Remark 2.9, in the index 1 case, one can locally solve the original DAE (2.14a) and (2.14c). Viewing such a solution (q, u) as fixed, one can subsequently locally solve for λ in equation (2.14d) as a function of p, since ∂φ/∂u is locally invertible. Substituting this into (2.14b) gives an ODE solely in the variable p, which can be solved locally.
Stated another way, if the original DAE (2.12a)-(2.12b) is an index 1 system, then the adjoint DAE system (2.14a)-(2.14d) is an index 1 system with dynamical variables (q, p) and algebraic variables (u, λ). To see this, if one denotes the constraints for the adjoint system (2.14c) and (2.14d) as then the matrix derivative ofφ with respect to the algebraic variables (u, λ) can be locally expressed in block form as where the block A has components given by the derivative of the right hand side of (2.14d) with respect to u. It is clear from the block triangular form of this matrix that it is pointwise invertible if ∂φ/∂u is.
Remark 2.11. It is clear from the coordinate expression (2.14a)-(2.14d) that a solution curve of the adjoint DAE system, if it exists, covers a solution curve of the original DAE system.
We now prove several results regarding the structure of the adjoint DAE system.
First, we show that the constraint equations (2.14c)-(2.14d) can be interpreted as the statement that the Hamiltonian H has the same time dependence as the "dynamical" Hamiltonian, when evaluated along a solution curve.
Proposition 2.9. For a solution curve (q, u, p, λ) of (2.13), Proof. For brevity, all functions below are appropriately evaluated along the solution curve. We where in the third equality, we used (2.14c) and (2.14d).
Remark 2.12. A more geometric way to view the above proposition is as follows: note that if a partially-defined vector field X exists such that i X Ω 0 = dH, then the change of H in a given direction Y , at any point where X is defined, can be computed as dH(Y ) = Ω 0 (X, Y ). Observe that the kernel of Ω 0 is locally spanned by ∂/∂u, ∂/∂λ, i.e., it is spanned by the coordinate vectors in the algebraic coordinates. Hence, the change of H in the algebraic coordinate directions is zero. This justifies referring to (u, λ) as "algebraic" variables.
We now show that the adjoint system (2.14a)-(2.14d) formally arises from a variational principle.
To do so, let Θ 0 be the pullback of the tautological one-form Θ d on the cotangent bundle T * M d by the maps Π andπ d , Θ 0 = Π * •π * d (Θ d ). Of course, one has Ω 0 = −dΘ 0 . Consider the action S defined by where ψ(t) = (q(t), u(t), p(t), λ(t)) is a curve on the generalized phase space bundle over the interval I = (t 0 , t 1 ). We consider the variational principle δS[ψ] = 0, subject to variations which fix the endpoints q(t 0 ), q(t 1 ). Proposition 2.10. Let ψ be a curve on the generalized phase space bundle over the interval I. Then, ψ is a stationary point of S with respect to variations which fix q(t 0 ), q(t 1 ) if and only if (2.14a)-(2.14d) hold.
Proof. In ψ = (q, u, p, λ) coordinates, the action has the expression The variation of the action reads δS[q, u, p, λ] · (δq, δu, δp, δλ) Remark 2.13. We will use the variational structure associated with the adjoint DAE system to construct numerical integrators in Section 3.2.
We now prove a result regarding the conservation of a quadratic invariant, analogous to the case of cotangent lifted adjoint systems in the ODE case. To do this, we define the variational equations as the linearization of the DAE (2.12a)-(2.12b). The coordinate expressions for the variational equations are obtained by taking the variation of equations (2.12a)-(2.12b) with respect to variations (δq, δu),q Proof. This follows from a direct computation, where we used (2.14b), (2.15c), (2.15d), and (2.14d).
Remark 2.14. Although we proved the previous proposition in coordinates, it can be understood intrinsically through the presymplecticity of the adjoint DAE flow. To see this, assume a partiallydefined vector field X exists such that i X Ω 0 = dH. Then, the flow of X preserves Ω 0 , which follows from The coordinate expression for the preservation of the presymplectic form Ω 0 = dq i ∧ dp i , with the appropriate choice of first variations, gives the previous proposition, analogous to the argument that we made in the symplectic (unconstrained) case.
Additionally, as we will see in Section 3.1, Proposition 2.11 will provide a method for computing adjoint sensitivities.
These two observations are interesting when constructing numerical methods to compute adjoint sensitivities, since if we can construct integrators that preserve the presymplectic form, then it will preserve the quadratic invariant and hence, be suitable for computing adjoint sensitivities efficiently.
Remark 2.15. For an index 1 DAE (2.12a)-(2.12b), since ∂φ/∂u is (pointwise) invertible for a fixed curve (q, u), one can solve for δu as a function of δq in the variational equation (2.15d) and substitute this into (2.15c) to obtain an explicit ODE for δq. Hence, in the index 1 case, given a solution (q, u) of the DAE (2.12a)-(2.12b) and an initial condition δq(0) in the tangent fiber over q(0), there is a corresponding (at least local) unique solution of the variational equations.
2.3.1. DAE Index and the Presymplectic Constraint Algorithm. In this section, we relate the index of the DAE (2.12a)-(2.12b) to the number of steps for convergence in the presymplectic constraint algorithm associated with the adjoint DAE system (2.13). In particular, we show that for an index 1 DAE, the presymplectic constraint algorithm for the associated adjoint DAE system converges after ν P = 1 step. Subsequently, we discuss how one can formally handle the more general index ν DAE case.
We consider again the presymplectic system given by the adjoint DAE system, P = T * M d ⊕ Φ * equipped with the presymplectic form Ω 0 = dq ∧ dp and Hamiltonian H(q, u, p, λ) = p, f (q, u) + λ, φ(q, u) , as discussed in the previous section. Our goal is to bound the number of steps in the presymplectic constraint algorithm ν P for this presymplectic system in terms of the index ν of the underlying DAE (2.12a)-(2.12b).
We now consider first the case when the DAE system (2.12a)-(2.12b) has index ν = 1 and subsequently, consider the general case ν ≥ 1.
The Presymplectic Constraint Algorithm for ν = 1. For the case ν = 1, we will show that the presymplectic constraint algorithm terminates after 1 step, i.e., ν P = ν = 1. Now, assume that the DAE system (2.12a)-(2.12b) has index ν = 1, i.e., for each (q, u) ∈ M d × M a such that φ(q, u) = 0, the matrix with A th row and a th column entry ∂φ A (q, u) ∂u a is invertible. Observe that the definition of the presymplectic constraint algorithm, equation (1.1), is local and hence, we seek a local coordinate expression for Ω 1 ≡ Ω 0 | P 1 and its kernel.
Let (q, u, p, λ) ∈ P 1 . In particular, φ(q, u) = 0. Since ∂φ(q, u)/∂u is invertible, by the implicit function theorem, one can locally solve for u as a function of q, which we denote u = u(q), such that φ(q, u(q)) = 0. Then, one can furthermore locally solve for λ as a function of q and p from the second constraint equation, Thus, we can coordinatize P 1 via coordinates (q ′ , p ′ ), where the inclusion i 1 : P 1 ֒ P is given by the coordinate expression ). Then, one obtains the local expression for Ω 1 , This is clearly nondegenerate, i.e., Z p = 0 for any Z ∈ ker(Ω 1 ), p ∈ P 1 , so the presymplectic constraint algorithm terminates, P 2 = P 1 . We conclude that ν P = 1.
To conclude the discussion of the index 1 case, we obtain coordinate expressions for the resulting nondegenerate Hamiltonian system. The Hamiltonian on P 1 can be expressed as

Thus, with the coordinate expression
We will now show explicitly that this Hamiltonian system solves (2.14a)-(2.14d) along the submanifold P 1 . Clearly, the latter two equations (2.14c)-(2.14d) are satisfied, by definition of P 1 . So, we want to show that the first two equations (2.14a)-(2.14b) are satisfied. Using the second constraint equation (2.14d), we have Substituting this into the equation forṗ ′ i above giveṡ Hence, the Hamiltonian system on P 1 can be equivalently expressed aṡ Thus, we have explicitly verified that (2.14a)-(2.14d) are satisfied along P 1 . Note that since the presymplectic constraint algorithm terminates at ν P = 1, X is guaranteed to be tangent to P 1 .
One can also verify this explicitly by computing the pushforward T i 1 (X) and verifying that it annihilates the constraint functions whose zero level set defines P 1 , 16. It is interesting to note that the Hamiltonian system i X Ω 1 = dH 1 , which we obtained by forming the adjoint system of the underlying index 1 DAE and subsequently, reducing the index of the adjoint DAE system through the presymplectic constraint algorithm, can be equivalently obtained (at least locally) by first reducing the index of the underlying DAE and then forming the adjoint system.
More precisely, if one locally solves φ(q, u) = 0 for u = u(q), then the index 1 DAE can be reduced to an ODE,q = f (q, u(q)). Subsequently, we can form the adjoint system to this ODE, as discussed in Section 2.2. The corresponding Hamiltonian is H(q, p) = p, f (q, u(q)) , which is the same as H 1 .
Thus, for the index 1 case, the process of forming the adjoint system and reducing the index commute.
Remark 2.17. In the language of the presymplectic constraint algorithm, Proposition 2.9 can be restated as the statement that the Hamiltonian H and its first derivatives, restricted to the primary constraint manifold, agrees with the dynamical Hamiltonian H 1 and its first derivatives.
Remark 2.18. An alternative view of the solution theory of the presymplectic adjoint DAE system (2.14a)-(2.14d) is through singular perturbation theory (see, for example, Berglund [6] and Chen and Trenn [13]). We proceed by writing (2.14a)-(2.14d) aṡ Applying a singular perturbation to the constraint equations yields the systeṁ where ǫ > 0. Observe that this is a nondegenerate Hamiltonian system with H(q, u, p, λ) as previously defined but with the modified symplectic form Ω ǫ = dq ∧ dp + ǫ du∧ dλ. Then, the above system can be expressed i X H Ω ǫ = dH. In the language of perturbation theory, the primary constraint manifold for the presymplectic system is precisely the slow manifold of the singularly perturbed system. One can utilize techniques from singular perturbation theory to develop a solution theory for this system, using Tihonov's theorem, whose assumptions for this particular system depend on the eigenvalues of the algebraic Hessian D 2 u,λ H (see, Berglund [6]). Although we will not elaborate on this here, this could be an interesting approach for the existence, stability, and approximation theory of such systems. In particular, the slow manifold integrators introduced in Burby and Klotz [9] may be relevant to their discretization. It is also interesting to note that for a solution (q ǫ , p ǫ , u ǫ , λ ǫ ) of the singularly perturbed system and a solution (δq ǫ , δu ǫ ) of the variational equations, one has the perturbed adjoint variational quadratic conservation law d dt p ǫ , δq ǫ + ǫ λ ǫ , δu ǫ = 0, which follows immediately from the preservation of Ω ǫ under the symplectic flow.
The Presymplectic Constraint Algorithm for General ν ≥ 1. Note that for the general case, we assume that the index of the DAE is finite, 1 ≤ ν < ∞.
In this case, there are two possible approaches to reduce the adjoint system: either form the adjoint system associated with the index ν DAE and then successively apply the presymplectic constraint algorithm or, alternatively, reduce the index of the DAE, form the adjoint system, and then apply the presymplectic constraint algorithm as necessary.
Since we have already worked out the presymplectic constraint algorithm for the index 1 case, we will take the latter approach. Namely, we reduce an index ν DAE to an index 1 DAE, and subsequently, apply the presymplectic constraint algorithm to the reduced index 1 DAE. Given an index ν DAE, it is generally possible to reduce the DAE to an index 1 DAE using the algorithm introduced in Mattsson and Söderlind [26]. The process of index reduction is given by differentiating the equations of the DAE to reveal hidden constraints. Geometrically, the process of index reduction can be understood as the successive jet prolongation of the DAE and subsequent projection back onto the first jet (see, Reid et al. [29]).
Thus, given an index ν DAEẋ =f (x, y),φ(x, y) = 0, we can, after ν − 1 reduction steps, transform it into an index 1 DAE of the formq = f (q, u), φ(q, u) = 0. Subsequently, we can form the adjoint DAE system and apply one iteration of the presymplectic constraint algorithm to obtain the underlying nondegenerate dynamical system. If we let the ν R,P denote the minimum number of DAE index reduction steps plus presymplectic constraint algorithm iterations necessary to take an index ν DAE and obtain the underlying nondegenerate Hamiltonian system associated with the adjoint, we have ν R,P ≤ ν.
Remark 2.19. Note that we could have reduced the index ν DAE to an explicit ODE after ν reduction steps, and subsequently, formed the adjoint. While this is formally equivalent to the above procedure by Remark 2.16, we prefer to keep the DAE in index 1 form. This is especially preferable from the viewpoint of numerics: if one reduces an index 1 DAE to an ODE and attempts to apply a numerical integrator, it is generically the case that the discrete flow drifts off the constraint manifold. For this reason, it is preferable to develop numerical integrators for the index 1 adjoint DAE system directly to prevent constraint violation.
where (q, u) ∈ R n × R m , f : R n × R m R n , g : R n R m , and ∂g ∂q ∂f ∂u is pointwise invertible. We reduce this to an index 1 DAE (2.12a)-(2.12b) as follows. Let M d = g −1 ({0}) be the dynamical configuration space which we will assume is a submanifold of R n . For example, this is true if g is a constant rank map. Furthermore, let M a = R m be the algebraic configuration space. To reduce the index, we differentiate the constraint g(q) = 0 with respect to time. This is equivalent to enforcing that the dynamics are tangent to M d . This gives u). Hence, we can form the semi-explicit index 1 system on M d × M a given bẏ q = f (q, u), 0 = φ(q, u).
The above system is an index 1 DAE since ∂φ ∂u = ∂g ∂q ∂f ∂u is pointwise invertible. We now form the adjoint DAE system associated with this index 1 DAE, (2.14a)-(2.14d). Expressing the constraint in terms of g and f , instead of φ, giveṡ We can then apply one iteration of the presymplectic constraint algorithm, as discussed above in the index ν = 1 case, to obtain the underlying nondegenerate Hamiltonian dynamics. Restricting to the primary constraint manifold, using the first constraint equation to solve for u = u(q) by the implicit function theorem and subsequently, using the second constraint equation to solve for λ = λ(q, p) by inverting ∂g ∂q ∂f ∂u T , gives the Hamiltonian systeṁ

2.3.2.
Adjoint Systems for DAEs with Augmented Hamiltonians. In Section 2.2.1, we augmented the adjoint ODE Hamiltonian by some function L. In this section, we do analogously for the adjoint DAE system.
To begin, let H(q, u, p, λ) = p, f (q, u) + λ, φ(q, u) be the Hamiltonian on the generalized phase space bundle corresponding to the DAEq = f (q, u), 0 = φ(q, u), and let L : M d × M a R be the function that we would like to augment. We identify L with its pullback through T * M d ⊕ Φ * M d × M a . Then, we define the augmented Hamiltonian (q, u, p, λ) H(q, u, p, λ) + L(q, u).
We define the augmented adjoint DAE system as the presymplectic system A direct calculation yields the coordinate expression, along an integral curve of such a (generally, partially-defined) vector field X H L ,q Remark 2.20. Observe that if the base DAE (2.12a)-(2.12b) has index 1, then the above system has index 1 by the exact same argument given in the nonaugmented case. After reduction by applying the presymplectic constraint algorithm and solving for u as a function of q and λ as a function of (q, p), the underlying nondegenerate Hamiltonian system on the primary (final) constraint manifold corresponds to the Hamiltonian which is the adjoint Hamiltonian for the ODEq ′ = f (q ′ , u(q ′ )), augmented by L(q ′ , u(q ′ )).
However, as we will discuss in Section 3.3, it is not uncommon in optimal control problems for ∂φ/∂u to be singular, but the presence of L dt in the minimization objective may uniquely specify the singular degrees of freedom.
We now prove an analogous proposition to Proposition 2.11, modified by the presence of L in the Hamiltonian. We again consider the variational equations (2.15a)-(2.15d) associated with the base DAE (2.12a)-(2.12b), which for simplicity we express in matrix derivative notation aṡ Proposition 2.12. For a solution (q, u, p, λ) of the augmented adjoint DAE system (2.17a)-(2.17d) and a solution (q, u, δq, δu) of the variational equations (2.18a)-(2.18d), covering the same solution (q, u) of the base DAE (2.12a)-(2.12b), Proof. This follows from a direct computation: where in the fourth equality above we used (2.18d) and in the sixth equality above we used (2.17d). Note that in our calculations below, the top row (the ODE case) can be formally obtained from the bottom row (the DAE case) simply by ignoring the algebraic variables (u, λ) and letting the constraint function φ be identically zero. Thus, we will focus on the bottom row, i.e., computing sensitivities of a terminal cost function and of a running cost function, subject to a DAE constraint.
In both cases, we will first show how the adjoint sensitivity can be derived using a traditional variational argument. Subsequently, we will show how the adjoint sensitivity can be derived more simply by using Propositions 2.11 and 2.12.

Adjoint Sensitivity of a Terminal Cost.
Consider the DAEq = f (q, u), 0 = φ(q, u) as in Section 2.3. We will assume that M d is a vector space and additionally, that the DAE has index 1. We would like to extract the gradient of a terminal cost function C(q(t f )) with respect to the initial condition q(0) = α, i.e., we want to extract the sensitivity of C(q(t f )) with respect to an infinitesimal perturbation in the initial condition, given by ∇ α C(q(t f )). Consider the functional J defined by Observe that for (q, u) satisfying the given DAE with initial condition q(0) = α, J coincides with C(q(t f )). We think of p 0 as a free parameter. For simplicity, we will use matrix derivative notation instead of indices. Computing the variation of J yields Integrating by parts in the term containing d dt δq and restricting to a solution (q, u, p, λ) of the adjoint DAE system (2.14a)-(2.14d) yields We enforce the endpoint condition p(t f ) = ∇ q C(q(t f )) and choose p 0 = p(0), which yields δJ = p(0), δα .
Hence, the sensitivity of C(q(t f )) is given by with initial condition q(0) = α and terminal condition p(t f ) = ∇ q C(q(t f )). Thus, the adjoint sensitivity can be computed by setting the terminal condition on p(t f ) above and subsequently, solving for the momenta p at time 0. In order for this to be well-defined, we have to verify that the given initial and terminal conditions lie on the primary constraint manifold P 1 . However, as discussed in Section 2.3.1, since the DAE has index 1, we can always solve for the algebraic variables u = u(q) and λ = λ(q, p) and thus, we are free to choose the initial and terminal values of q and p, respectively. For higher index DAEs, one has to ensure that these conditions are compatible with the final constraint manifold. For example, this is done in [11] in the case of Hessenberg index 2 DAEs. Alternatively, at least theoretically, for higher index DAEs, one can reduce the DAE to an index 1 DAE and then the above discussion applies, however, this reduction may fail in practice due to numerical cancellation.
Note that the above adjoint sensitivity result is also a consequence of the preservation of the quadratic invariant p, v as in Proposition 2.11. From this proposition, one has that where δq satisfies the variational equations. Setting p(t f ) = ∇ q C(q(t f )) and δq(0) = δα gives the same result. As mentioned in Remark 2.14, this quadratic invariant arises from the presymplecticity of the adjoint DAE system. Thus, a numerical integrator which preserves the presymplectic structure is desirable for computing adjoint sensitivities, as it exactly preserves the quadratic invariant that allows the adjoint sensitivities to be accurately and efficiently computed. We will discuss this in more detail in Section 3.2.
Adjoint Sensitivity of a Running Cost. Again, consider an index 1 DAEq = f (q, u), 0 = φ(q, u). We would like to extract the sensitivity of a running cost function where L : M d × M a R, with respect to an infinitesimal perturbation in the initial condition q(0) = α. Consider the functional J defined by Observe that when the DAE is satisfied with initial condition q(0) = α, J = t f 0 L dt. Now, we would to compute the implicit change in t f 0 L dt with respect to a perturbation δα in the initial condition. Taking the variation in J yields Restricting to a solution (q, u, p, λ) of the augmented adjoint DAE system (2.17a)-(2.17d), setting the terminal condition p(t f ) = 0, and choosing p 0 = p(0) gives δJ = p(0), δα . Hence, the implicit sensitivity of t f 0 L dt with respect to a change δα in the initial condition is given by Thus, the adjoint sensitivity of a running cost functional with respect to a perturbation in the initial condition can be computed by using the augmented adjoint DAE system (2.17a)-(2.17d) with terminal condition p(t f ) = 0 to solve for the momenta p at time 0.
Note that the above adjoint sensitivity result can be obtained from Proposition 2.12 as follows. We write equation ( to highlight that the right hand side measures the total induced variation of L. Now, we integrate this equation from 0 to t f , which gives Since we want to determine the change in the running cost functional with respect to a perturbation in the initial condition, we set p(t f ) = 0 which yields The right hand side is the total change induced on the running cost functional, whereas the left hand side tells us how this change is implicitly induced from a perturbation δq(0) in the initial condition. Note that a perturbation in the initial condition δq(0) will generally induce perturbations in both q and u, according to the variational equations. Such a curve (δq, δu) satisfying the variational equations exists in the index 1 case as noted in Remark 2.15. Thus, we arrive at the same conclusion as the variational argument: p(0) is the desired adjoint sensitivity.
To summarize, adjoint sensitivities for terminal and running costs can be computed using the properties of adjoint systems, such as the various aforementioned propositions regarding d dt p, δq , which is zero in the nonaugmented case and measures the variation of L in the augmented case. In the case of a terminal cost, one sets an inhomogeneous terminal condition p(t f ) = ∇ q C(q(t f )) and backpropagates the momenta through the nonaugmented adjoint DAE system (2.14a)-(2.14d) to obtain the sensitivity p(0). On the other hand, in the case of a running cost, one sets a homogeneous terminal condition p(t f ) = 0 and backpropagates the momenta through the augmented adjoint DAE system (2.17a)-(2.17d) to obtain the sensitivity p(0).
The various propositions used to derive the above adjoint sensitivity results are summarized below. We also include the ODE case, since it follows similarly. In Section 3.2, we will construct integrators that admit discrete analogues of the above propositions, and hence, are suitable for computing discrete adjoint sensitivities. [23] to construct structurepreserving integrators which admit discrete analogues of Propositions 2.3, 2.7, 2.11, and 2.12, and are therefore suitable for numerical adjoint sensitivity analysis. For brevity, the proofs of these discrete analogues can be found in Appendix A.

Structure-Preserving Discretizations of Adjoint Systems. In this section, we utilize the Galerkin Hamiltonian variational integrators of Leok and Zhang
We start by recalling the construction of Galerkin Hamiltonian variational integrators as introduced in Leok and Zhang [23]. We assume that the base manifold Q is a vector space and thus, we have the identification T * Q ∼ = Q × Q * . To construct a variational integrator for a Hamiltonian system on T * Q, one starts with the exact Type II generating function where one extremizes over C 2 curves on the cotangent bundle satisfying q(0) = q 0 , p(∆t) = p 1 . This is a Type II generating function in the sense that it defines a symplectic map (q 0 , p 1 ) (q 1 , p 0 ) by q 1 = D 2 H + d,exact (q 0 , p 1 ), p 0 = D 1 H + d,exact (q 0 , p 1 ). To approximate this generating function, one approximates the integral above using a quadrature rule and extremizes the resulting expression over a finite-dimensional subspace satisfying the prescribed boundary conditions. This yields the Galerkin discrete Hamiltonian where ∆t > 0 is the timestep, q 0 , q 1 , p 0 , p 1 are numerical approximations to q(0), q(∆t), p(0), p(∆t), respectively, b i > 0 are quadrature weights corresponding to quadrature nodes c i ∈ [0, 1], Q i and P i are internal stages representing q(c i ∆t), p(c i ∆t), respectively, and V is related to Q by Q i = q 0 + ∆t j a ij V j , where the coefficients a ij arise from the choice of function space. The expression above is extremized over the internal stages Q i , P i and subsequently, one applies the discrete right Hamilton's equations , to obtain a Galerkin Hamiltonian variational integrator. The extremization conditions and the discrete right Hamilton's equations can be expressed as where we interpret a ij as Runge-Kutta coefficients andã ij = (b i b j − b j a ji )/b i as the symplectic adjoint of the a ij coefficients. Thus, (3.1a)-(3.1d) can be viewed as a symplectic partitioned Runge-Kutta method.
We will consider such methods in four cases: adjoint systems corresponding to a base ODE or DAE, and whether or not the corresponding system is augmented. Note that in the DAE case, we will have to modify the above construction because the system is presymplectic. Furthermore, we will assume that all of the relevant configuration spaces are vector spaces.
Nonaugmented Adjoint ODE System. The simplest case to consider is the nonaugmented adjoint ODE system (2.6a)-(2.6b). Since the quadratic conservation law in Proposition 2.3, d dt p, δq = 0, arises from symplecticity, a structure-preserving discretization can be obtained by applying a symplectic integrator. This case is already discussed in Sanz-Serna [33], so we will only outline it briefly.
Applying the Galerkin Hamiltonian variational integrator (3.1a)-(3.1d) to the Hamiltonian for the adjoint ODE system, H(q, p) = p, f (q) , yields In the setting of adjoint sensitivity analysis of a terminal cost function, the appropriate boundary condition to prescribe on the momenta is p 1 = ∇ q C(q(t f )), as discussed in Section 3.1.
Since the above integrator is symplectic, we have the symplectic conservation law, dq 1 ∧ dp 1 = dq 0 ∧ dp 0 , when evaluated on discrete first variations of (3.2a)-(3.2d). In this setting, a discrete first variation can be identified with solutions of the linearization of (3.2a)-(3.2d). For the linearization of the equations in the position variables, (3.2a)-(3.2b), we have As observed in Sanz-Serna [33], while we obtained this by linearizing the discrete equations, one could also obtain this by first linearizing (2.1) and subsequently, applying the Runge-Kutta scheme to the linearization. For the linearization of the equations for the adjoint variables, (3.2c)-(3.2d), observe that they are already linear in the adjoint variables, so we can identify the linearization with itself. Thus, we can choose for first variations vector fields V as the first variation corresponding to the solution of the linearized position equation and W as the first variation corresponding to the solution of the adjoint equation itself. With these choices, the above symplectic conservation law yields 0 = dq 1 ∧ dp 1 (V, W )| (q 1 ,p 1 ) − dq 0 ∧ dp 0 (V, W )| (q 0 ,p 0 ) = p 1 , δq 1 − p 0 , δq 0 . This is of course a discrete analogue of Proposition 2.3. Note that one can derive the conservation law p 1 , δq 1 = p 0 , δq 0 directly by starting with the expression p 1 , δq 1 and substituting the discrete equations where appropriate. We will do this in the more general augmented case below.
Augmented Adjoint ODE System. We now consider the case of the augmented adjoint ODE system (2.11a)-(2.11b). In the continuous setting, we have from Proposition 2.7, d dt p, δq = − dL, δq .
We would like to construct an integrator which admits a discrete analogue of this equation. To do this, we apply the Galerkin Hamiltonian variational integrator, equations (3.1a)-(3.1d), to the augmented Hamiltonian H L (q, p) = p, f (q) + L(q). This gives We now prove a discrete analogue of Proposition 2.7. To do this, we again consider the discrete variational equations for the position variables, (3.3a)-(3.3b).
Proposition 3.1. With the above notation, the above integrator satisfies

Proof. See Appendix A.
Remark 3.1. To see that this is a discrete analogue of d dt p, δq = − dL, δq , we write it in integral form as Then, applying the quadrature rule on [0, ∆t] given by quadrature weights b i ∆t and quadrature nodes c i ∆t, the above integral is approximated by which yields equation (3.5). The discrete analogue is natural in the sense that the quadrature rule for which the discrete equation (3.5) approximates the continuous equation is the same as the quadrature rule used to approximate the exact discrete generating function. This occurs more generally for such Hamiltonian variational integrators, as noted in Tran and Leok [37] for the more general setting of multisymplectic Hamiltonian variational integrators.
For adjoint sensitivity analysis of a running cost L dt, the appropriate boundary condition to prescribe on the momenta is p 1 = 0, as discussed in Section 3.1. With such a boundary condition, equation (3.5) reduces to p 0 , δq 0 = ∆t Thus, p 0 gives the discrete sensitivity, i.e., the change in the quadrature approximation of L dt induced by a change in the initial condition along a discrete solution trajectory. One can compute this quantity directly via the direct method, where one needs to integrate the discrete variational equations for every desired search direction δq 0 . On the other hand, by the above proposition, one can compute this quantity using the adjoint method: one integrates the adjoint equation with p 1 = 0 once to compute p 0 and subsequently, pair p 0 with any search direction δq 0 to obtain the sensitivity in that direction. By the above proposition, both methods give the same sensitivities. However, assuming the search space has dimension n > 1, the adjoint method is more efficient since it only requires O(1) integrations and O(n) vector-vector products, whereas the direct method requires O(n) integrations and O(ns) vector-vector products where s ≥ 1 is the number of Runge-Kutta stages, since, in the direct method, one has to compute dL(Q i ), δQ i for each i and for each choice of δq 0 .
Nonaugmented Adjoint DAE System. We will now construct discrete Hamiltonian variational integrators for the adjoint DAE system (2.14a)-(2.14d), where we assume that the base DAE has index 1. To construct such a method, we have to modify the Galerkin Hamiltonian variational integrator (3.1a)-(3.1d), so that it is applicable to the presymplectic adjoint DAE system.
First, consider a general presymplectic system i X Ω ′ = dH. Note that, locally, any presymplectic system can be transformed to the canonical form (see, Cariñena et al. [12]), where, in these coordinates, Ω ′ = dq∧dp, so that ker(Ω ′ ) = span{∂/∂r}. The action for this system is given by ∆t 0 ( p,q − H(q, p, r))dt. We approximate this integral by quadrature, introduce internal stages for q, p as before, and additionally introduce internal stages R i = r(c i h). This gives the discrete generating function where again V is related to the internal stages of Q by Q i = q 0 + ∆t j a ij V j and the above expression is extremized over the internal stages Q i , P i , R i . The discrete right Hamilton's equations are again given by , which we interpret as the evolution equations of the system. There are no evolution equations for r due to the presymplectic structure and the absence of derivatives of r in the action. This gives the integrator where (3.6b), (3.6d), (3.6e) arise from extremizing with respect to P i , Q i , R i , respectively, while (3.6a) and (3.6c) arise from the discrete right Hamilton's equations. This integrator is presymplectic, in the sense that dq 1 ∧ dp 1 = dq 0 ∧ dp 0 , when evaluated on discrete first variations. The proof is formally identical to the symplectic case. For this reason, we refer to (3.6a)-(3.6e) as a presymplectic Galerkin Hamiltonian variational integrator.
Remark 3.2. In general, the system (3.6a)-(3.6e) evolves on the primary constraint manifold given implicitly by the zero level set of D r H, however, it may not evolve on the final constraint manifold. This is not an issue for us since we are dealing with adjoint DAE systems for index 1 DAEs, for which we know the primary constraint manifold and the final constraint manifold coincide. For the general case, one may need to additionally differentiate the constraint equation D r H = 0 to obtain hidden constraints.
Thus, the method (3.6a)-(3.6e) is generally only applicable to index 1 presymplectic systems, unless we add in further hidden constraints. In order for the continuous presymplectic system to have index 1, it is sufficient that the Hessian of H with respect to the algebraic variables, D 2 r H, is (pointwise) invertible on the primary constraint manifold. This is the case for the adjoint DAE system corresponding to an index 1 DAE.
We now specialize to the adjoint DAE system (2.14a)-(2.14d), corresponding to an index 1 DAE, which is already in the above canonical form with r = (u, λ) and H(q, u, p, λ) = p, f (q, u) + λ, φ(q, u) . Note that we reordered the argument of H, (q, p, r) = (q, p, u, λ) (q, u, p, λ), in order to be consistent with the previous notation used throughout. We label the internal stages for the algebraic variables as R i = (U i , Λ i ). Applying the presymplectic Galerkin Hamiltonian variational integrator to this particular system yields where (3.7b), (3.7d), (3.7e), (3.7f) arise from extremizing over P i , Q i , Λ i , U i , respectively, while (3.7a), (3.7c) arise from the discrete right Hamilton's equations. Remark 3.3. In order for q 1 to appropriately satisfy the constraint, we should take the final quadrature point to be c s = 1 (for an s-stage method), so that φ(q 1 , U s ) = φ(Q s , U s ) = 0. In this case, equation (3.7a) and equation (3.7b) with i = s are redundant. Note that with the choice c s = 1, they are still consistent (i.e., are the same equation), since in the Galerkin construction, the coefficients a ij and b i are defined as where φ j are functions on [0, 1] which interpolate the nodes c j (see, Leok and Zhang [23]). Hence, a sj = b j , so that the two equations are consistent. However, we will write the system as above for conceptual clarity. Furthermore, even in the case where one does not take c s = 1, the proposition that we prove below still holds, despite the possibility of constraint violations.
A similar remark holds for the adjoint variable p and the associated constraint (3.7f), except we think of p 0 as the unknown, instead of p 1 .
Proposition 3.2. With the above notation, the above integrator satisfies Thus, the above integrator admits a discrete analogue of Proposition 2.11 for the nonaugmented adjoint DAE system. By setting p 1 = ∇ q C(q(t f )), one can use this integrator to compute the sensitivity p 0 of a terminal cost function with respect to a perturbation in the initial condition. As discussed before, this only requires O(1) integrations instead of O(n) integrations via the direct method (for a dimension n search space). Furthermore, the adjoint method requires only O(1) numerical solves of the constraints, while the direct method requires O(n) numerical solves.
Remark 3.4. Since we are assuming the DAE has index 1, it is always possible to prescribe an arbitrary initial condition q 0 (and δq 0 ) and terminal condition p 1 , since the corresponding algebraic variables can always formally be solved for using the corresponding constraints. In practice, one generally has to solve the constraints to some tolerance, e.g., through an iterative scheme. If the constraints are only satisfied to a tolerance O(ǫ), then the above proposition holds to O(sǫ), where s is the number of Runge-Kutta stages.
Remark 3.5. The above method (3.7a)-(3.7f) is presymplectic, since it is a special case of the more general presymplectic Galerkin Hamiltonian variational integrator (3.6a)-(3.6e). Although we proved it directly, the above proposition could also have been proven from presymplecticity, with the appropriate choices of first variations.
The presymplectic integrator is then The associated variational equations are again (3.8a)-(3.8c). Remarks analogous to the nonaugmented case regarding setting the quadrature node c s = 1 and solvability of these systems under the index 1 assumption can be made. Proposition 3.3. With the above notation, the above integrator satisfies

Proof. See Appendix A.
Remark 3.6. Analogous to the remark in the augmented adjoint ODE case, the above proposition is a discrete analogue of Proposition 2.12, in integral form, The discrete analogue is natural in the sense that it is just quadrature applied to the right hand side of this equation, with the same quadrature rule used to discretize the generating function. Remark 3.7. As with the augmented adjoint ODE case, the above proposition allows one to compute numerical sensitivities of a running cost function by solving for p 0 with p 1 = 0, which is more efficient than the direct method.
To summarize, we have utilized Galerkin Hamiltonian variational integrators to construct methods which admit natural discrete analogues of the various propositions used for sensitivity analysis. We summarize the results below.

Terminal Cost
Running Cost 3.2.1. Naturality of the Adjoint DAE System Discretization. To conclude our discussion of discretizing adjoint systems, we prove a discrete extension of the fact that, for an index 1 DAE, the process of index reduction and forming the adjoint system commute, as discussed in Section 2.3.1. Namely, we will show that, starting from an index 1 DAE (2.12a)-(2.12b), the processes of reduction, forming the adjoint system, and discretization all commute, for particular choices of these processes which we will define and choose below. This can be summarized in the following commutative diagram. In the above diagram, we will use the convention that the "Discretize" arrows point forward, the "Adjoint" arrows point downward, and the "Reduce" arrows point to the right. For the "Discretize" arrows on the top face, we take the discretization to be a Runge-Kutta discretization (of a DAE on the left and of an ODE on the right, with the same Runge-Kutta coefficients in both cases). For the "Discretize" arrows on the bottom face, we take the discretization to be the symplectic partitioned Runge-Kutta discretization induced by the discretization of the base DAE or ODE, i.e., the momenta expansion coefficientsã ij are the symplectic adjoint of the coefficients a ij used on the top face. We have already defined the "Adjoint" arrows on the back face, as discussed in Section 2.
For the "Adjoint" arrows on the front face, we define them as forming the discrete adjoint system corresponding to a discrete (and generally nonlinear) system of equations and we will review this notion where needed in the proof. We have already defined the "Reduce" arrows on the back face, as discussed in Section 2.3.1. For the "Reduce" arrows on the front face, we define this as solving for the discrete algebraic variables in terms of the discrete kinematic variables through the discrete constraint equations. With these choices, the above diagram commutes, as we will show. To prove this, it suffices to prove that the diagrams on each of the six faces commutes. To keep the exposition concise, we provide the proof in Appendix B and move on to discuss the implications of this result.
The previous discussion shows that the presymplectic Galerkin Hamiltonian variational integrator construction is natural for discretizing adjoint (index 1) DAE systems, in the sense that the integrator is equivalent to the integrator produced from applying a symplectic Galerkin Hamiltonian variational integrator to the underlying nondegenerate Hamiltonian system. Of course, in practice, one cannot generally determine the function u = u(q) needed to reduce the DAE to an ODE. Therefore, one generally works with the presymplectic Galerkin Hamiltonian variational integrator instead, where one iteratively solves the constraint equations. However, although reduction then symplectic integration is often impractical, one can utilize this naturality to derive properties of the presymplectic integrator. For example, we will use this naturality to prove a variational error analysis result.
The basic idea for the variational error analysis result goes as follows: one utilizes the naturality to relate the presymplectic variational integrator to a symplectic variational integrator of the underlying nondegenerate Hamiltonian system and subsequently, applies the variational error analysis Consider the following optimal control problem in Bolza form, subject to a DAE constraint, which we refer to as (OCP-DAE), where the DAE systemq = f (q, u), 0 = φ(q, u) is over M d × M a as described in Section 2.3, R is the running cost, the initial condition q(0) = q 0 is prescribed, and for generality, a terminal constraint φ f (q(t f )) = 0 is also imposed, where φ f is a map from M d into some vector space V .
We assume a local optimum to (OCP-DAE). We then adjoin the constraints to J using adjoint variables, which gives the adjoined functional The optimality conditions are given by the condition that J is stationary about the local optimum, δJ = 0 (Biegler [7]). For simplicity in the notation, we will use matrix derivative instead of indices. Note also that we will implicitly leave out the variation of the adjoint variables, since those terms pair with the DAE constraints, which vanish at the local optimum. The optimality condition δJ = 0 is then where we integrated by parts on the term p, d dt δq and used δq(0) = 0 since the initial condition is fixed. Enforcing stationarity for all such variations gives the optimality conditions, The first four optimality conditions (3.10a)-(3.10d) are precisely the augmented adjoint DAE equations, (2.17a)-(2.17d). The last two optimality conditions (3.10e), (3.10f) are the terminal constraint and the associated transversality condition, respectively. Note that these conditions are only sufficient for a trajectory (q, u, p, λ) to be an extremum of the optimal control problem; whether or not the trajectory is optimal depends on the properties of the DAE constraint and cost function, e.g., convexity of L.
Regular Index 1 Optimal Control. In the literature, the problem (OCP-DAE) is usually formulated by making a distinction between algebraic variables and control variables, (q, y, u), instead of (q, u) (see, for example, Biegler [7] and Aguiar et al. [2]). This does not change any of the previous discussion of the optimality conditions, except that (3.10d) splits into two equations for y and u. That is, the distinction is not formally important for the previous discussion. It is of course important when actually solving such an optimal control problem. For example, the constraint function φ(q, y, u) may have a singular matrix derivative with respect to (y, u) but may have a nonsingular matrix derivative with respect to y. In such a case, one interprets y as the algebraic variable, in that it can locally be solved in terms of (q, u) via the constraint, and the control variable u as "free" to optimize over. We now briefly elaborate on this case.
We take the configuration manifold for the algebraic variables to be M a = Y a × U ∋ (y, u), where y is interpreted as the algebraic constraint variable and u is interpreted as the control variable. We will assume that the control space U is compact. The constraint has the form φ(q, y, u) = 0, and we assume that ∂φ/∂y is pointwise invertible. We consider the following optimal control problem, We perform an analogous argument to before, except that, in this case, since U may have a boundary, the optimality for the control variable u will either require u to lie on ∂U or will require the stationarity of the adjoined functional with respect to variations in u. In any case, the necessary conditions for optimality can be expressed aṡ where H L is the augmented Hamiltonian H L (q, y, u) = L(q, y, u) + p, f (q, y, u) + λ, φ(q, y, u) . Assuming that u lies in the interior of U , (3.11e) can be expressed as 0 = ∇ u L(q, y, u) + [D u f (q, y, u)] * p + [D u φ(q, y, u)] * λ, or D u H L (q, y, u) = 0. We say that an optimal control problem with a DAE constraint forms a regular index 1 system if both ∂φ/∂y and the Hessian D 2 u H L are pointwise invertible. In this case, whenever u lies on the interior of U , (y, u, λ) can be locally solved as functions of (q, p). Thus, in principle, the resulting Hamiltonian ODE for (q, p) can be integrated to yield extremal trajectories for the optimal control problem. As mentioned before, without additional assumptions on the DAE and cost function, such a trajectory will only generally be an extremum but not necessarily optimal.
Of course, in practice, one cannot generally analytically integrate the resulting ODE nor determine the functions which give (y, u, λ) in terms of (q, p). Thus, the only practical option is to discretize the presymplectic system above to compute approximate extremal trajectories. To integrate such a presymplectic system, one can again use the presymplectic Galerkin Hamiltonian variational integrator construction discussed in Section 3.2. Such an integrator would be natural in the following sense. First, as discussed in Section 3.2, a presymplectic Galerkin Hamiltonian variational integrator applied to the augmented adjoint DAE system is equivalent to applying a symplectic Galerkin Hamiltonian variational integrator to the underlying Hamiltonian ODE, with the same Runge-Kutta expansions for q 1 , Q i in both methods. Furthermore, as shown in Sanz-Serna [33], utilizing a symplectic integrator to discretize the extremality conditions is equivalent to first discretizing the ODE constraint by a Runge-Kutta method and then enforcing the associated discrete extremality conditions. This also holds in the DAE case.
More precisely, beginning with a regular index 1 optimal control problem, the processes of reduction, extremization, and discretization commute, for suitable choices of these processes, analogous to those used in the naturality result discussed in Section 3.2.1. The proof is similar to the naturality result discussed in Section 3.2.1, where the arrow given by forming the adjoint is replaced by extremization. In essence, these are the same, since the extremization condition is given by the adjoint system, so we will just elaborate briefly. We already know how to extremize the continuous optimal control problem, with either a DAE constraint or an ODE constraint after reduction, which results in an adjoint system. We also already know how to discretize the resulting adjoint system after discretization, using a (pre)symplectic partitioned Runge-Kutta method. Furthermore, at any step, reduction is just defined to be solving the continuous or discrete constraints for y in terms of (q, u). Thus, the only major difference compared to the previous naturality result is defining the discretization of the optimal control problem and subsequently, how to extremize the discrete optimal control problem. For the regular index 1 optimal control problem, min t f 0 L(q, y, u)dt subject tȯ q = f (q, y, u), 0 = φ(q, y, u), q 0 = q(0), its discretization is obtained by replacing the constraints with a Runge-Kutta discretization and replacing the cost function with its quadrature approximation, using the same quadrature weights as those in the Runge-Kutta discretization. This can be written as where Q i = q 0 + ∆t j a ij V j , which implicitly encodes q(0) = q 0 . One can then extremize this discrete system, which is given by the discrete Euler-Lagrange equations for the discrete action That is, we enforce the discrete constraints by adding to the discrete Lagrangian the appropriate Lagrange multiplier terms paired with the constraints, where we weighted the Lagrange multipliers P i , Λ i by ∆tb i just as convention, in order to interpret them as the appropriate variables, as discussed in Appendix B. Enforcing extremality of this action recovers a partitioned Runge-Kutta method applied to the adjoint system corresponding to extremizing the continuous optimal control problem, as discussed in Appendix B, where the Runge-Kutta coefficients for the momenta are the symplectic adjoint of the original Runge-Kutta coefficients. Alternatively, starting from the original continuous optimal control problem, one could first reduce the DAE constraint to an ODE constraint using the invertibility of D y φ to give min t f 0 L(q, y(q, u), u)dt subject tȯ q = f (q, y(q, u), u), One can then discretize this using the same Runge-Kutta method as before, where the cost function is replaced with a quadrature approximation, and then extremize using Lagrange multipliers. Alternatively, one can extremize the continuous problem to yield an adjoint system and then apply a partitioned Runge-Kutta method to that system, where the momenta Runge-Kutta coefficients are again the symplectic adjoint of the original Runge-Kutta coefficients. Having defined all of these processes, a direct computation yields that all of the processes commute, analogous to the computation in Appendix B.

Conclusion and Future Research Directions
In this paper, we utilized symplectic and presymplectic geometry to study the properties of adjoint systems associated with ODEs and DAEs, respectively. The (pre)symplectic structure of these adjoint systems led us to a geometric characterization of the adjoint variational quadratic conservation law used in adjoint sensitivity analysis. As an application of this geometric characterization, we constructed structure-preserving discretizations of adjoint systems by utilizing (pre)symplectic integrators, which led to natural discrete analogues of the quadratic conservation laws.
A natural research direction is to extend the current framework to adjoint systems for differential equations with nonholonomic constraints, in order to more generally allow for constraints between configuration variables and their derivatives. In this setting, it is reasonable to expect that the geometry of the associated adjoint systems can be described using Dirac structures (see, for example, Yoshimura and Marsden [40,41]), which generalize the symplectic and presymplectic structures of adjoint ODE and DAE systems, respectively. Structure-preserving discretizations of such systems could then be studied through the lens of discrete Dirac structures (Leok and Ohsawa [22]). These discrete Dirac structures make use of the notion of a retraction (Absil et al. [1]). The tangent and cotangent lifts of a retraction also provide a useful framework for constructing geometric integrators (Barbero-Liñán and Martín de Diego [3]). It would be interesting to synthesize the notion of tangent and cotangent lifts of retraction maps with discrete Dirac structures in order to construct discrete Dirac integrators for adjoint systems with nonholonomic constraints which generalize the presymplectic integrators constructed in Barbero-Liñán and Martín de Diego [4].
Another natural research direction is to extend the current framework to evolutionary partial differential equations (PDEs). There are two possible approaches in this direction. The first is to consider evolutionary PDEs as ODEs evolving on infinite-dimensional spaces, such as Banach or Hilbert manifolds. One can then investigate the geometry of the infinite-dimensional symplectic structure associated with the corresponding adjoint system. In practice, adjoint systems for evolutionary PDEs are often formed after semi-discretization, leading to an ODE on a finite-dimensional space. Understanding the reduction of the infinite-dimensional symplectic structure of the adjoint system to a finite-dimensional symplectic structure under semi-discretization could provide useful insights into structure-preservation. The second approach would be to explore the multisymplectic structure of the adjoint system associated with a PDE. This approach would be insightful for several reasons. First, an adjoint variational quadratic conservation law arising from multisymplecticity would be adapted to spacetime instead of just time. With appropriate spacetime splitting and boundary conditions, such a quadratic conservation law would induce either a temporal or spatial conservation law. As such, one could use the multisymplectic conservation law to determine adjoint sensitivities for a PDE with respect to spatial or temporal directions, which could be useful in practice [24]. Furthermore, the multisymplectic framework would apply equally as well to nonevolutionary (elliptic) PDEs, where there is no interpretation of a PDE as an infinite-dimensional evolutionary ODE. Additionally, adjoint systems for PDEs with constraints could be investigated with multi-Dirac structures (Vankerschaver et al. [38]). In future work, we aim to explore both approaches, relate them once a spacetime splitting has been chosen, and investigate structurepreserving discretizations of such systems by utilizing the multisymplectic variational integrators constructed in Tran and Leok [37].
where, in the last equality, we substituted (3.4d) and (3.3b). We now group and simplify the above expression, where, in the last equality, we used (3.3b).

Proof of Proposition 3.2.
Proof. For brevity, we denote Starting from p 1 , δq 1 , we substitute the evolution equations (3.7c), (3.7d), (3.8a), (3.8b), where in the third to last equality, we used the constraint equation (3.7f) and in the last equality, we used the constraint equation (3.8c).

Proof of Proposition 3.3.
Proof. The proof uses computations analogous to those used in the proofs of Back Face. We have already proved that the back face commutes (i.e., that reduction and forming the adjoint commute when starting with an index 1 DAE), as discussed in Section 2.3.1. One can then interpret the above diagram as an extension of this result with an extra dimension corresponding to discretization.
Right Face. This was proven in Sanz-Serna [33]. One can then interpret the above diagram as an extension of the result in Sanz-Serna [33] by adding the reduction operation.
Bottom Face. Consider the augmented adjoint DAE system corresponding to the DAE (2.12a)-(2.12b), which we take to have index 1, i.e., ∂φ/∂u is pointwise invertible. We consider the augmented case because the nonaugmented case can be obtained by taking L ≡ 0. We show that reducing the system first and then applying a symplectic Galerkin Hamiltonian variational integrator is equivalent to applying a presymplectic Galerkin Hamiltonian variational integrator, with the same partitioned Runge-Kutta coefficients, and then reducing.
We start with the former approach. The symplectic adjoint ODE system given by reduction, as discussed in Section 2.3.1, is the Hamiltonian system corresponding to the Hamiltonian where we have solved u = u(q ′ ) and defined f ′ (q ′ ) ≡ f (q ′ , u(q ′ )), L ′ (q ′ ) ≡ L(q ′ , u(q ′ )). Applying the symplectic Galerkin Hamiltonian variational integrator construction yields the integrator Note that the derivative Df ′ can be equivalently expressed as where D i denotes differentiation with respect to the i th argument. We switch to indexing the derivative operator here, so we do not have to make the distinction between total derivatives D q and partial derivatives ∂ q . Similarly, we can express dL ′ as follows. First, note that we have been implicitly identifying the row vector dL ′ with the column vector given by its transpose ∇L ′ . Thus, dL ′ in equations (B.1c)-(B.1d) should really be written as ∇L ′ . Thus, Now, we show that the second approach is equivalent to the above system. The starting point is the presymplectic Galerkin Hamiltonian variational integrator, equations (3.9a)-(3.9f). From (3.9e), we can solve for U i in terms of Q i as U i = u(Q i ). Plugging this into (3.9a)-(3.9b) gives precisely (B.1a)-(B.1b). Thus, we just need to see that, after solving the constraint (3.9f) for Λ i , the two momenta equations (3.9c)-(3.9d) are equivalent to (B.1c)-(B.1d). Solving (3.9f) for Λ i gives

Multiplying both sides by
where in the second equality, we used D 1 φ(Q i , u(Q i )) = −D 2 φ(Q i , u(Q i ))Du(Q i ) from the implicit function theorem. Plugging this expression and U i = u(Q i ) into (3.9c)-(3.9d) yields (B.1c)-(B.1d), noting the above expressions for Df ′ , dL ′ .
Remark B.1. Note that, in the above, we used the implicit function theorem to obtain the local function u = u(q). This is sufficient to prove that the two processes are the same for a single integration step, assuming that the timestep ∆t is sufficiently small and the vector field f and constraint φ are sufficiently regular, so that q 0 , q 1 , and all of the internal stages Q i are in the neighborhood where the local function is defined. For each subsequent time step, one generally needs a different local function. This does not matter in practice since one works directly with the presymplectic integrator and solves the constraints iteratively.
Top Face. We want to prove that, starting from an index 1 DAE, the processes of discretization and reduction commute, where the discretization of the ODE and DAE have the same Runge-Kutta coefficients.
We start first with reduction then discretization. Starting from the index 1 DAEq = f (q, u), φ(q, u) = 0, we apply the reduction operation, which gives the ODEq = f (q, u(q)). Applying a Runge-Kutta discretization gives On the other hand, we can discretize the DAE and then reduce. We discretize the DAEq = f (q, u), φ(q, u) = 0 by applying a Runge-Kutta discretization with the same coefficients as before, To reduce this system, we solve the constraint equations U i = u(Q i ) and substitute these into the two evolution equations, which yields the same system obtained from first reducing and then discretizing.
Front Face. The starting point for this loop is a discrete DAE system, which arises as a Runge-Kutta discretization of an index 1 DAE, i.e., it is given by the discrete system From here, we wish to show that reducing and forming the discrete adjoint system commute.
First, we recall the notion of a discrete adjoint system. Suppose we are given a generally nonlinear system of equations, F (x 1 ) = x 0 , where x 1 ∈ V is unknown, x 0 ∈ W is given, and F : V W (where V and W are vector spaces). To define the adjoint system, we first consider the variational equations associated with this nonlinear system given by its linearization, where DF (x 1 ) is a linear map V W and δx 0 ∈ W is given. Suppose that we are interested in computing the quantity s 1 , δx 1 for a given vector s 1 ∈ V * . In the setting of adjoint sensitivity analysis, the quantity s 1 , δx 1 is the sensitivity of the terminal cost function. We define the associated adjoint equation as [DF (x 1 )] * s 0 = s 1 .
Thus, to compute s 1 , δx 1 , one could solve the variational equation for δx 1 and pair it with s 1 which is given, or, alternatively, solve the adjoint equation for s 0 and pair it with δx 0 which is given, since these linear systems are solvable by assumption. We define the adjoint system associated with the equation F (x 1 ) = x 0 as this equation combined with the associated adjoint equation, i.e., as the combined system [Df (x 1 )] * s 0 = s 1 .
Following Ibragimov [20], we will utilize an alternative characterization of the adjoint system. We define the discrete adjoint action S(x 1 , s 0 ) ≡ s 0 , F (x 1 ) .
Then, observe that S is a generating function for the adjoint system (x 1 , s 0 ) (x 0 , s 1 ), in the sense that x 0 = δ δs 0 S(x 1 , s 0 ) = F (x 1 ), This characterization serves two purposes. First, it will simplify the calculation of the adjoint system for the case at hand. Furthermore, it resembles the process of forming the adjoint at the continuous level: starting from the (discrete or continuous) differential(-algebraic) equation at hand, one forms the (discrete or continuous) adjoint action and applies the variational principle to obtain the adjoint system. To obtain the augmented adjoint system, we add a discrete Lagrangian L : V R to the action (as a convention, we subtract the discrete Lagrangian). We define the augmented discrete adjoint action to be S L (x 1 , s 0 ) ≡ s 0 , F (x 1 ) − L(x 1 ).
The map that this generates defines the augmented discrete adjoint system, Observe that this definition of an augmented discrete adjoint system is natural in the sense that, s 1 , δx 1 = [Df (x 1 )] * s 0 + dL(x 1 ), δx 1 = s 0 , x 0 − dL(x 1 ), δx 1 , which resembles the continuous analogue of the adjoint sensitivity result for a running cost function. Now, we use this notion of a discrete adjoint system for the problem at hand. We begin first with reduction and then forming the adjoint system. Applying the reduction operation to the discrete DAE system (B.2a)-(B.2c), given by solving φ(Q i , U i ) = 0 for U i = u(Q i ), Let us define Q i = q 0 + ∆t j a ij V j . We think of the internal stages Q as functions of the internal stages V , which are the internal stage proxies forq. Our discrete system (B.3a)-(B.3b) can then be defined by where s is the number of internal stages, and Observe that F = 0 only gives the internal stage equations (B.3b). We do this for simplicity, since we will assume c s = 1 as is typical for a Runge-Kutta discretization of a DAE as previously discussed and hence, equation (B.3a) is redundant, since a sj = b j .
We define F and x 1 in terms of V instead of Q because when we form the adjoint action, we pair the components of F with the dual variable s 0 . In order to interpret s 0 as representing the momenta internal stages P i , it should be paired with the proxy for the tangent vector V , instead of Q. We now form the discrete adjoint action. We define the dual variable for the adjoint system to be s 0 = {∆tb i P i } s i=1 . The normalization factor ∆tb i is used so that the discrete action is the quadrature approximation of the continuous action. This is just a convention, but we would have to reinterpret the components of s 0 if we did not choose this convention. Finally, we define the discrete Lagrangian to be the quadrature approximation of the continuous Lagrangian L ′ (q) ≡ L(q, u(q)), i.e., L(x 1 ) = ∆t i b i L ′ (Q i (V )). This is the natural choice because the discrete sensitivity of a running cost function is ∆t i b i dL ′ (Q i (V )), δQ i (V ) , which equals dL(x 1 ), δx 1 with the above choice of L. The augmented discrete adjoint action is then To define the discrete adjoint system, we have to give s 1 , which we take to be s 1 = {∆tb i p 1 } s i=1 , where p 1 is given. Thus, the augmented discrete adjoint system is given by , u(Q i (V )))Du(Q i (V ))] * P i + dL ′ (Q i (V )) .
The first set of equations above, combined with the definition of Q in terms of V , gives (B.3b). For the second set of equations, we first divide through by ∆tb k and rearrange to obtain Note that this is the usual symplectic partitioned Runge-Kutta expansion for the internal stages P i , expressed in terms of p 1 instead of p 0 . Thus, the full adjoint system, combined with the redundant k = s stages, yields a symplectic partitioned Runge-Kutta method.
Now, in the other direction, we first form the adjoint system corresponding to the discrete DAE system and subsequently reduce. We begin by forming the adjoint system. We form the discrete action analogously to before, but now the discrete system (B.2a)-(B.2c) also has constraints which we must incorporate into F , since we have not yet reduced the system. We take We define F as Note again that Q is a function of V as Q i = q 0 + ∆t j a ij V j . It is not a priori a function of U because the condition V i = f (Q i (V ), U i ) has not yet been enforced. Rather, it is a consequence of the variational principle, which formally matters when one computes the variation of the discrete action. Define the discrete Lagrangian L(x 1 ) = i b i L(Q i (V ), U i ). We form the augmented discrete adjoint action We use this as a generating function to compute the adjoint system as before. The computation is analogous so we will just state the result, Finally, we reduce by solving the last two equations for U i , Λ i as functions of Q i (V i ), P i . Finally, an implicit function theorem computation analogous to the proof of the bottom face shows that this is the same as the system obtained by first reducing and then forming the discrete adjoint.
Left Face. The proof for the left face is formally similar to the right face, but since we have already computed both directions, we will include it for completeness. Starting from an index 1 DAE, forming the adjoint and then discretizing just give the presymplectic Galerkin Hamiltonian variational integrator (3.9a)-(3.9f). In the other direction, we first discretize the DAE and then take the adjoint which we did in the proof of the front face. Expressed in terms of Q, instead of V , this is Q k = q 0 + ∆t j a ij f (Q j , U j ), Returning to the system given by first forming the adjoint and then discretizing, (3.9a)-(3.9f), one substitutes (3.9c) into (3.9d) to write the internal stages for P i in terms of p 1 , and this gives the above system.