Dynamical reduced basis methods for Hamiltonian systems

We consider model order reduction of parameterized Hamiltonian systems describing nondissipative phenomena, like wave-type and transport dominated problems. The development of reduced basis methods for such models is challenged by two main factors: the rich geometric structure encoding the physical and stability properties of the dynamics and its local low-rank nature. To address these aspects, we propose a nonlinear structure-preserving model reduction where the reduced phase space evolves in time. In the spirit of dynamical low-rank approximation, the reduced dynamics is obtained by a symplectic projection of the Hamiltonian vector field onto the tangent space of the approximation manifold at each reduced state. A priori error estimates are established in terms of the projection error of the full model solution onto the reduced manifold. For the temporal discretization of the reduced dynamics we employ splitting techniques. The reduced basis satisfies an evolution equation on the manifold of symplectic and orthogonal rectangular matrices having one dimension equal to the size of the full model. We recast the problem on the tangent space of the matrix manifold and develop intrinsic temporal integrators based on Lie group techniques together with explicit Runge–Kutta (RK) schemes. The resulting methods are shown to converge with the order of the RK integrator and their computational complexity depends only linearly on the dimension of the full model, provided the evaluation of the reduced flow velocity has a comparable cost.

Hamiltonian systems can be viewed as dynamical extension of the first law of thermodynamics. In this work, we consider parameterized finite-dimensional canonical Hamiltonian systems: these can model energy-conserving nondissipative flows or can ensue from the numerical discretization of partial differential equations derived from action principles. Many relevant models in mathematical physics can be written as Hamiltonian systems, and find application in, for example, classical mechanics, quantum dynamics, population and epidemics dynamics. Furthermore, partial differential equations that can be derived from action principles include Maxwell's equations, Schrödinger's equation, Korteweg-de Vries and the wave equation, compressible and incompressible Euler equations, Vlasov-Poisson and Vlasov-Maxwell equations.
Our target problem is as follows. Let T := (t 0 , T ] be a temporal interval and let V 2N be a 2N -dimensional vector space. Let ⊂ R d , with d ≥ 1, be a compact set of parameters. For each η ∈ , we consider the initial value problem: For u 0 (η) ∈ V 2N , find u(·, η) ∈ C 1 (T , V 2N ) such that for t ∈ T , u(t 0 , η) = u 0 (η), (1.1) where X H (u, η) ∈ V 2N is the Hamiltonian vector field at time t ∈ T , and C 1 (T , V 2N ) denotes continuous differentiable functions in time taking values in V 2N . Numerical simulations of systems like (1.1) can become prohibitively expensive, in terms of computational cost, if the number 2N of degrees of freedom is large. In the context of long-time and many-query simulations, this often leads to unmanageable demands on computational resources. Model order reduction aims at alleviating this computational burden by replacing the original high-dimensional problem with a low-dimensional, efficient model that is fast to solve but that approximates well the underlying fullorder dynamics. When dealing with Hamiltonian systems additional difficulties are encountered to ensure that the geometric structure of the phase space, the stability and the conservation properties of the original system are not hindered during the reduction. The main goal of this work is to develop and analyze structure-preserving model order reduction methods for the efficient, accurate, and physically consistent approximation of high-dimensional parametric Hamiltonian systems. Within model order reduction techniques, projection-based reduced basis methods (RBM) consist in building, during a computationally intensive offline phase, a reduced basis from a proper orthogonal decomposition of a set of high-fidelity simulations (referred to as snapshots) at sampled values of time and parameters. A reduced dynamics is then obtained via projection of the full model onto the lower dimension space spanned by the reduced basis. Projection-based RBM for Hamiltonian systems tailored to preserve the geometric structure of the dynamics were developed in [21] and [6] using a variational Lagrangian formulation of the problem, in [2,4,30] for canonically symplectic dynamical systems, and in [16] to deal with Hamiltonian problems whose phase space is endowed with a state-dependent Poisson manifold structure. Although the aforementioned approaches can provide robust and efficient reduced models, they might require a sufficiently large approximation space to achieve even moderate accuracy. This can be ascribed to the fact that nondissipative phenomena, like advection and wave-type problems, do not possess a global low-rank structure, and are therefore characterized by slowly decaying Kolmogorov widths, as highlighted in [11]. Hence, local reduced spaces seem to provide a more effective instrument to deal with this kind of dynamical systems.
In this work we propose a nonlinear projection-based model order reduction of parameterized Hamiltonian systems where the reduced basis is dynamically evolving in time. The idea is to consider a modal decomposition of the approximate solution to (1.1) of the form where the reduced basis {U i } 2n i=1 ⊂ R 2N , and the expansion coefficients {Z i } 2n i=1 ⊂ R can both change in time. The approximate reduced flow is then generated by the velocity field resulting from the projection of the vector field X H in (1.1) into the tangent space of the reduced space at the current state. By imposing that the evolving reduced space spanned by {U i } 2n i=1 is a symplectic manifold at every time the continuous reduced dynamics preserves the geometric structure of the full model.
Low-rank approximations based on a modal decomposition of the approximate solution with dynamically evolving modes similar to (1.2), have been widely studied in quantum mechanics in the multiconfiguration time-dependent Hartree (MCTDH) method, see e.g. [23]. In the finite dimensional setting, a similar approach, known as dynamical low-rank approximation [20], provides a low-rank factorization updating technique to efficiently compute approximations of time-dependent large data matrices, by projecting the matrix time derivative onto the tangent space of the low-rank matrix manifold. For the discretization of time-dependent stochastic PDEs, Sapsis and Lermusiaux proposed in [31] the so-called dynamically orthogonal (DO) scheme, where the deterministic approximation space adapts over time by evolving according to the differential operator describing the stochastic problem. A connection between dynamical low-rank approximations and DO methods was established in [29]. Further, a geometric perspective on the relation between dynamical low-rank approximation, DO field equations and model order reduction in the context of time-dependent matrices has been investigated in [14]. To the best of our knowledge, the only work to address structure-preserving dynamical low-rank approximations is [28], where the authors develop a DO discretization of stochastic PDEs possessing a symplectic Hamiltonian structure. The method proposed in [28] consists in recasting the continuous PDE into the complex setting and then applying a dynamical low-rank strategy to derive field equations for the evolution of the stochastic modal decomposition of the approximate solution. The approach we propose for the nonlinear model order reduction of problem (1.1) adopts a geometric perspective similar to [14] and yields an evolution equation for the reduced solution analogous to [28], although we do not resort to a reformulation of the evolution problem in a complex framework.
Concerning the temporal discretization of the reduced dynamics describing the evolution of the approximate solution (1.2), the low-dimensional system for the expansion coefficients {Z i } 2n i=1 is Hamiltonian and can be approximated using standard symplectic integrators. On the other hand, the development of numerical schemes for the evolution of the reduced basis is more involved as two major challenges need to be addressed: (i) a structure-preserving approximation requires that the discrete evolution remains on the manifold of symplectic and (semi-)orthogonal rectangular matrices; (ii) since the reduced basis forms a matrix with one dimension equal to the size of the full model, the effectiveness of the model reduction might be thwarted by the computational cost associated with the numerical solution of the corresponding evolution equation. Various methods have been proposed in the literature to solve differential equations on manifolds, see e.g. [15,Chapter IV]. Most notably projection methods apply a conventional discretization scheme and, after each time step, a "correction" is made by projecting the updated approximate solution to the constrained manifold. Alternatively, methods based on the use of local parameterizations of the manifold, so-called intrinsic, are well-developed in the context of differential equations on Lie groups, cf. [15,Sect. IV.8]. The idea is to recast the evolution equation in the corresponding Lie algebra, which is a linear space, and to then recover an approximate solution in the Lie group via local coordinate maps. Instrinsic methods possess excellent structurepreserving properties provided the local coordinate map can be computed exactly. However, they usually require a considerable computational cost associated with the evaluation of the coordinate map and its inverse at every time step (possibly at every stage within each step).
We propose and analyze two structure-preserving temporal approximations and show that their computational complexity scales linearly with the dimension of the full model, under the assumption that the velocity field of the reduced flow can be evaluated at a comparable cost. The first algorithm we propose is a Runge-Kutta Munthe-Kaas (RK-MK) method [24], and we rely on the action on the orthosymplectic matrix manifold by the quadratic Lie group of unitary matrices. By exploiting the structure of our dynamical low-rank approximation and the properties of the local coordinate map supplied by the Cayley transform, we prove the computational efficiency of this algorithm with respect to the dimension of the high-fidelity model. However, a polynomial dependence on the number of stages of the RK temporal integrator might yield high computational costs in the presence of full models of moderate dimension. To overcome this issue, we propose a discretization scheme based on the use of retraction maps to recast the local evolution of the reduced basis on the tangent space of the matrix manifold at the current state, inspired by the works [9,10] on intrinsic temporal integrators for orthogonal flows.
The remainder of the paper is organized as follows. In Sect. 2 the geometric structure underlying the dynamics of Hamiltonian systems is presented, and the concept of orthosymplectic basis spanning the approximate phase space is introduced. In Sect. 3 we describe the properties of linear symplectic maps needed to guarantee that the geometric structure of the full dynamics is inherited by the reduced problem. Subsequently, in Sect. 4 we develop and analyze a dynamical low-rank approximation strategy resulting in dynamical systems for the reduced orthosymplectic basis and the corresponding expansion coefficients in (1.2). In Sect. 5 efficient and structurepreserving temporal integrators for the reduced basis evolution problem are derived. Section 6 concerns a numerical test where the proposed method is compared to a global reduced basis approach. We present some concluding remarks and open questions in Sect. 7.

Hamiltonian dynamics on symplectic manifolds
The phase space of Hamiltonian dynamical systems is endowed with a differential Poisson manifold structure which underpins the physical properties of the system. Most prominently, Poisson structures encode a family of conserved quantities that, by Noether's theorem, are related to symmetries of the Hamiltonian. Here we focus on dynamical systems whose phase space has a global Poisson structure that is canonical and nondegenerate, namely symplectic.
On a finite 2N -dimensional smooth manifold V 2N , let ω be a 2-form, that is, for any p ∈ V 2N , the map ω p : T p V 2N × T p V 2N → R is skew-symmetric and bilinear on the tangent space to V 2N at p, and it varies smoothly in p. The 2-form ω is a symplectic structure if it is closed and ω p is symplectic for all p ∈ V 2N , in the sense of Definition 2.1. A manifold V 2N endowed with a symplectic structure ω is called a symplectic manifold and denoted by (V 2N , ω). The algebraic structure of a symplectic manifold (V 2N , ω) can be characterized through the definition of a bracket: Let dF be the 1-form given by the exterior derivative of a given smooth function F. Then, for where T * V 2N ·, · T V 2N denotes the duality pairing between the cotangent and the tangent bundle. The function J 2N : To any function H ∈ C ∞ (V 2N ), the symplectic form ω allows to associate a vector field X H ∈ T V 2N , called Hamiltonian vector field, via the relation where i denotes the contraction operator. Since ω is nondegenerate, X H ∈ T V 2N is unique. Any vector field X H on a manifold V 2N determines a phase flow, namely a oneparameter group of diffeomorphisms t The flow of a Hamiltonian vector field satisfies ( t X H ) * ω = ω, for each t ∈ T , that is t X H is a symplectic diffeomorphism (symplectomorphism) on its domain. : In addition to possessing a symplectic phase flow, Hamiltonian dynamics is characterized by the existence of differential invariants, and symmetry-related conservation laws.

Definition 2.3 (Invariants of motion)
The Hamiltonian, if time-independent, is an invariant of motion. A particular subset of the invariants of motion of a dynamical system is given by the Casimir invariants, smooth functions C on V 2N that {·, ·} 2N -commute with every other functions, i.e. {C, F} 2N = 0 for all F ∈ C ∞ (V 2N ). Since Casimir invariants are associated with the center of the Lie algebra (C ∞ (V 2N ), {·, ·} 2N ), symplectic manifolds only possess trivial Casimir invariants.
Resorting to a coordinate system, the canonical structure on a symplectic manifold can be characterized by canonical charts whose existence is postulated in [1, In the local canonical coordinates introduced in Definition 2.4, the vector bundle map J 2N , defined in (2.1), takes the canonical symplectic form where Id and 0 denote the identity and zero map, respectively. Symplectic canonical charts on a symplectic vector space allow to identify a Kälher structure, namely a compatible combination of a scalar product and symplectic form, as follows. On a symplectic vector space (V 2N , ω), the operator J 2N is an almost complex structure, that is a linear map on V 2N such that J 2N • J 2N = − Id . Furthermore, J 2N is compatible with the symplectic structure ω, namely, for any u, v ∈ V 2N , u = 0, it holds A symplectic form ω on a vector space V 2N together with a compatible positive almost complex structure J 2N determines an inner product on V 2N , given by A symplectic basis on (V 2N , ω) is an orthonormal basis for the compatible inner product (2.3), and we refer to it as orthosymplectic. A subspace U of a symplectic vector space (V 2N , ω) is called Lagrangian if it coincides with its symplectic complement in V 2N , namely the set of u ∈ V 2N such that ω(u, v) = 0 for all v ∈ U. As a consequence of the fact that any basis of a Lagrangian subspace of a symplectic vector space can be extended to a symplectic basis, every symplectic vector space admits an orthosymplectic basis, cf. for example [5,Sect. 1.2]. With the definitions introduced hitherto, we can recast the dynamical system (1.1) on a symplectic vector space (V 2N , ω) as a Hamiltonian initial value problem. For each η ∈ , and for u 0 (η) where H(·, η) ∈ C ∞ (V 2N ) is the Hamiltonian function, and ∇ u denotes the gradient with respect to the variable u. The well-posedness of (2.4) is guaranteed by assuming that, for any fixed η ∈ , the operator X H : is Lipschitz continuous in u uniformly in t ∈ T in a suitable norm.

Orthosymplectic matrices
In order to construct surrogate models preserving the physical and geometric properties of the original Hamiltonian dynamics we build approximation spaces of reduced dimension endowed with the same geometric structure of the full model. To this aim, the reduced space is constructed as the span of suitable symplectic and orthonormal time-dependent bases, so that the reduced space inherits the geometric structure of the original dynamical system. In this Section we describe the properties of linear symplectic maps between finite dimensional symplectic vector spaces. Analogously to [1, p. 168], we can easily extend the characterization of symplectic linear maps to the case of vector spaces of different dimension as in the following result. We define symplectic right inverse of the symplectic matrix M + ∈ R 2n×2N the matrix M = J 2N M + J 2n ∈ R 2N ×2n . It can be easily verified that M + M = I 2n , and that M : is the adjoint operator with respect to the symplectic form ω, i.e. ω(M + u, y) = ω(u, My) for any u ∈ V 2N , and y ∈ V 2n . Furthermore, the symplectic condition M + J 2N M + = J 2n is equivalent to M J 2N M = J 2n . Owing to this equivalence, with a small abuse of notation, we will say that M ∈ R 2N ×2n is symplectic if it belongs to the space Sp 2n, R 2N := L ∈ R 2N ×2n : L J 2N L = J 2n .
Orthosymplectic rectangular matrices can be characterized as follows.
Conversely, the symplecticity of In order to design numerical methods for evolution problems on the manifold U(2n, R 2N ) of orthosymplectic rectangular matrices, we will need to characterize its tangent space. To this aim we introduce the vector space so(2n) of skew-symmetric 2n × 2n real matrices so(2n) := {M ∈ R 2n×2n : M + M = 0 2n }, and the vector space sp(2n) of Hamiltonian 2n ×2n real matrices, namely sp(2n) := {M ∈ R 2n×2n : M J 2n + J 2n M = 0 2n }. Throughout, if not otherwise specified, we will denote with G 2n := U(2n) the Lie group of orthosymplectic 2n × 2n matrices and with g 2n the corresponding Lie algebra g 2n := so(2n) ∩ sp(2n), with bracket given by the matrix

Orthosymplectic dynamical reduced basis method
Assume we want to solve the parameterized Hamiltonian problem (2.4) at p ∈ N samples of the parameter {η j } p j=1 =: h ⊂ R pd . To simplify the notation we take d = 1, namely we assume that the parameter η is a scalar quantity, for vector-valued η the derivation henceforth applies mutatis mutandis. Then, the Hamiltonian system (2.4) can be recast as a set of ordinary differential equations in a 2N × p matrix unknown. Let η h ∈ R p denote the vector of sampled parameters, the evolution problem reads: (4.1) Let n N , to characterize the reduced solution manifold we consider an approximation of the solution of (4.1) of the form . . , 2n, and j = 1, . . . , p. Since we aim at a structure-preserving model order reduction of (4.1), we impose that the basis U (t) is orthosymplectic at all t ∈ T , in analogy with the symplectic reduction techniques employing globally defined reduced spaces. Here, since U is changing in time, this means that we constrain its evolution to the manifold U(2n, R 2N ) from Definition 3.2. With this in mind, the reduced solution is sought in the reduced space defined as If the full-rank condition (4.5) is satisfied, then the number p of samples of the parameter η ∈ satisfies p ≥ n. This means that, for a fixed p, a too large reduced basis might lead to a violation of the full rank condition, which would entail a rank-deficient evolution problem for the coefficient matrix Z ∈ R p×2n . This is related to the problem of overapproximation in dynamical low-rank techniques, see [20,Sect. 5.3]. Observe also that if p ≥ 2n and rank(Z ) = 2n then the full rank condition (4.5) is always satisfied. In general, the elements of M spl 2n might not have full rank 2n: for any R ∈ M spl 2n it holds rank(Z ) ≤ rank(R) ≤ min{2n, p}.
is a fiber bundle with fibers given by the group of unitary matrices U(2n), and M spl 2n is isomorphic to (M/ U (2n)) × V p×2n . Indeed, let U 1 ∈ M and Z 1 ∈ V p×2n , then, for any arbitrary A ∈ U(2n), it holds In dynamically orthogonal approximations [31] a characterization of the reduced solution is obtained by fixing a gauge constraint in the tangent space of the reduced solution manifold. For the manifold M spl 2n the tangent space at R ∈ M spl 2n is defined as the set of X ∈ R 2N × p such that there exists a differentiable path γ : 2n is of the form X =U Z + UŻ , whereU andŻ denote the time derivatives of U (t) and Z (t), respectively. Taking the derivative of the orthogonality constraint on U yieldsU U + U U = 0. Analogously, the symplecticity constraint giveṡ (4.6) However, this parameterization is not unique. Indeed, let S ∈ g 2n be arbitrary: if X U U ∈ g 2n then the matrix (X U + U S) U belongs to g 2n , and the pairs (X U , X Z ) and (X U + U S, X Z + Z S) identify the same tangent vector X := X U Z + U X Z . We fix the parameterization of the tangent space as follows.
This means that the map . By the definition of the tangent space (4.6), the zero vector admits the representation 0 which implies X U = 0 in view of the full-rank condition (4.5).
For the surjectivity of we show that Any X ∈ T U Z M spl 2n can be written as X =U Z + UŻ whereŻ ∈ R p×2n anḋ U ∈ R 2N ×2n satisfiesU U ∈ g 2n . Hence, the tangent vector X can be recast as We need to show that the pair (X U , X Z ), defined as X U := (I 2N − UU )U and X Z :=Ż + ZU U , belongs to the space H (U ,Z ) . From the orthogonality of U it easily follows that , and use the symplectic constraint on U and its temporal derivative to get Then, using the commutativity of the symplectic unit J 2N and the projection onto the orthogonal complement to the space spanned by U , i.e.  .7), is a horizontal space. This decomposition into the subset of directions tangent to the fiber and its complementary space provides a unique parameterization of the tangent space. We refer the reader to e.g. [13] and [19,Chapter 2], for further details on the topic.
Owing to Proposition 4.1, the tangent space of M spl 2n can be characterized as Henceforth, we consider M endowed with the metric induced by the ambient space C 2N ×2n , namely the Frobenius inner product A, B := tr(A H B), where A H denotes the conjugate transpose of the complex matrix A, and we will denote with · the Frobenius norm. Note that, on simple Lie algebras, the Frobenius inner product is a multiple of the Killing form.

Dynamical low-rank symplectic variational principle
For any fixed η ∈ , the vector field X H in (2.4) at time t belongs to T u(t) V 2N . Taking the cue from dynamical low-rank approximations [20], we derive a dynamical system on the reduced space M spl 2n via projection of the velocity field X H of the full dynamical system (4.1) onto the tangent space of M spl 2n at the current state. The reduced dynamical system is therefore optimal in the sense that the resulting vector field is the best dynamic approximation of X H , in the Frobenius norm, at every point on the manifold V 2N . To preserve the geometric structure of the full dynamics we construct a projection which is symplectic for each value of the parameter η j ∈ h , with 1 ≤ j ≤ p. To this aim, let us introduce on the symplectic vector space (V 2N , ω) the family of skew-symmetric bilinear forms ω j : where a j ∈ R 2N denotes the j-th column of the matrix a ∈ R 2N × p , and similarly for b j ∈ R 2N .
where ω j is defined in (4.8).
is a projection.
It can be easily verified that where X U = X U (v) and X Z = X Z (v), but henceforth we omit the dependence on v.
Using the definition of X Z and the symplecticity of the basis U the last term becomes Moreover, it can be easily checked that ω(X U Z , U Y Z ) = 0 by definition of X U and by the orthosymplecticity of U . Hence, the only non-trivial terms are We need to prove that T 1 and T 2 coincide. Let M i ∈ R 2N denote the i-th column vector of a given matrix M ∈ R 2N ×2n . The properties of the symplectic canonical form ω yield To deal with the term Moreover, using once more the fact that Y U ∈ H U results in The result follows by definition of X U (v).

Remark 4.4
Owing to the inner product structure (2.3), the projection operator from Proposition 4.3 is orthogonal in the Frobenius norm since This means that the projection gives the best low-rank approximation of the velocity vector, and hence the reduced dynamics is associated with the flow field ensuing from the best approximation in the tangent space to the reduced manifold.
To compute the initial condition of the reduced problem, we perform the complex SVD of R 0 (η h ) ∈ R 2N × p truncated at the n-th mode. Then the initial value U 0 ∈ M is obtained from the resulting unitary matrix of left singular vectors of R 0 (η h ) by exploiting the isomorphism between M and St(n, C N ), cf. Lemma 4.8. The expansion coefficients matrix is initialized as Z 0 = R 0 (η h ) U 0 . Therefore, the dynamical system for the approximate reduced solution (4.2) reads: for t ∈ T , For any 1 ≤ j ≤ p and t ∈ T , let Z j (t) ∈ R 1×2n be the j-th row of the matrix , and ∇ U Z j denotes the gradient with respect to U Z j . Using the decomposition R = U Z in (4.3), we can now derive from (4.9) evolution equations for U and Z : The reduced problem (4.10) is analogous to the system derived in [28,Proposition 6.9]. The evolution equations for the coefficients Z form a system of p equations in 2n unknowns and correspond to the Galerkin projection onto the space spanned by the columns of U , as obtained with a standard reduced basis method. Here, however, the projection is changing over time as the reduced basis U is evolving. For U fixed, the flow map characterizing the evolution of each Z j , for 1 ≤ j ≤ p, is a symplectomorphism (cf. Definition 2.2), i.e. the dynamics is canonically Hamiltonian. The evolution problem satisfied by the basis U is a matrix equation in 2N × 2n unknowns on the manifold of orthosymplectic rectangular matrices introduced in Definition 3.2, as shown in the following result.
Proof We first show that, for any matrix by the assumption on the initial condition. Moreover, the conditionẆ = J 2NẆ J 2n together with the dynamical orthogo- Owing to the reasoning above, we only need to verify that the solution of (4.10) satisfiesU ∈ H U . The dynamical orthogonal conditionU U = 0 is trivially satisfied. Moreover, since S J 2n = J 2n S, the constraintU = J 2NU J 2n is satisfied ifU S J 2n = J 2NU S. One can easily show that

Remark 4.6
Observe that the dynamical reduced basis technique proposed in the previous Section can be extended to more general Hamiltonian systems endowed with a degenerate constant Poisson structure. The idea is to proceed as in [16,Sect. 3] by splitting the dynamics into the evolution on a symplectic submanifold of the phase space and the trivial evolution of the Casimir invariants. The symplectic dynamical model order reduction developed in Sect. 4 can then be performed on the symplectic component of the dynamics.

Conservation properties of the reduced dynamics
The velocity field of the reduced flow (4.9) is the symplectic projection of the full model velocity onto the tangent space of the reduced manifold. For any fixed parameter η j ∈ h , let H j := H(·, η j ). In view of Proposition 4.3, the reduced solution R ∈ This implies that the Hamiltonian H is a conserved quantity of the continuous reduced problem (4.10). Indeed, To deal with the other invariants of motion, let us assume for simplicity that p = 1.
Since the linear map R 2N → span{U (t)} associated with the reduced basis at any time t ∈ T cannot be symplectic, the invariants of motion of the full and reduced model cannot be in one-to-one correspondence. Nevertheless, a result analogous to [16,Lemma 3.9] holds.

Lemma 4.7
Let π * +,t be the pullback of the linear map associated with the reduced basis U (t) at time t ∈ T . Assume that H ∈ Im(π * +,t ) for any t ∈ T . Then,

Convergence estimates with respect to the best low-rank approximation
In order to derive error estimates for the reduced solution of problem (4.9), we extend to our setting the error analysis of [14,Sect. 5] which shows that the error committed by the dynamical approximation with respect to the best low-rank approximation is bounded by the projection error of the full model solution onto the reduced manifold of low-rank matrices. To this aim, we resort to the isomorphism between the reduced symplectic manifold M spl 2n defined in (4.3) and the manifold M n of rank-n complex matrices, already established in [28, Lemma 6.1]. Then, we derive the dynamical orthogonal approximation of the resulting problem in the complex setting and prove that it is isomorphic to the solution of the reduced Hamiltonian system (4.9). The differentiability properties of orthogonal projections onto smooth embedded manifolds and the trivial extension to complex matrices of the curvature bounds in [14] allows to derive an error estimate.
Let L( ) denote the set of functions with values in the vector space , and let F : Then, problem (4.1) can be recast in the complex setting as: For C(t 0 ) ∈ M n associated with R(t 0 ) ∈ M spl 2n via the map (4.13), we can therefore derive the DO dynamical system: find C ∈ C 1 (T , M n ) such thaṫ for t ∈ T , (4.14) where T C M n is the projection onto the tangent space of M n at C = W Y , defined as The so-called dynamically orthogonal condition X H W W = 0, allows to uniquely parameterize the tangent space T C M n by imposing that the complex reduced basis evolves orthogonally to itself.
Let M * indicate the complex conjugate of a given matrix M. The projection onto the tangent space of M n can be characterized as in the following result. Lemma 4.9 At every C = W Y ∈ M n , the map is the · -orthogonal projection onto the tangent space of M n at C.
Proof The result can be derived similarly to the proof of [14,Proposition 7] by minimizing the convex functional J(X W , Using the expression (4.15) for the projection onto the tangent space of M n , we can derive from (4.14) evolution equations for the terms W and Y : Given (4.16)

Proposition 4.10
Under the assumption of well-posedness, problem (4.9) is equivalent to problem (4.14).

Proof
The proof easily follows from algebraic manipulations of the field equations (4.10) and (4.16) and from the definition of the isomorphism (4.13).
In view of Proposition 4.10, we can revert to the error estimate established in [14].
Theorem 4.11 ([14,Theorem 32]). Let C ∈ C 1 (T , C N × p ) denote the exact solution of (4.12) and let C ∈ C 1 (T , M n ) be the solution of (4.14) at time t ∈ T . Assume that no crossing of the singular values of C occurs, namely Let M n be the · -orthogonal projection onto M n . Then, at any time t ∈ T , the error between the approximate solution C(t) and the best rank-n approximation of C(t) can be bounded as where μ ∈ R and ν ∈ R are defined as and L X ∈ R denotes the Lipschitz continuity constant of X H .
The remainder of this work pertains to numerical methods for the temporal discretization of the reduced dynamics (4.10). Since we consider splitting techniques, see e.g. [15,Sect. II.5], the evolution problems for the expansion coefficients and for the reduced basis are examined separately. The coefficients Z (t) ∈ V p×2n of the expansion (4.2) satisfy a Hamiltonian dynamical system (4.10) in the reduced symplectic manifold of dimension 2n spanned by the evolving orthosymplectic basis U (t) ∈ M. The numerical approximation of the evolution equation for Z (t) can, thus, be performed using symplectic integrators, cf. [15,Sect. VI]. Observe that the use of standard splitting techniques might require the approximate reduced solution, at a given time step, to be projected into the space spanned by the updated basis. This might cause an error in the conservation of the invariants due to the projection step, that, however, can be controlled under sufficiently small time steps. In principle, exact conservation can be guaranteed if the evolution of the reduced basis evolves smoothly at the interface of temporal interval (or temporal subintervals associated with the splitting), or, in other words, if the splitting is synchronous and the two systems are concurrently advanced in time. We postpone to future work the investigation and the numerical study of splitting methods that exactly preserve the Hamiltonian.

Numerical methods for the evolution of the reduced basis
Contrary to global projection-based model order reduction, dynamical reduced basis methods eschew the standard online-offline paradigm. The construction and evolution of the local reduced basis (4.10) does not require queries of the high-fidelity model so that the method does not incur a computationally expensive offline phase. However, the evolution of the reduced basis entails the solution of a matrix equation in which one dimension equals the size of the full model. Numerical methods for the solution of (4.10) will have arithmetic complexity min{C R , C F } where C F is the computational cost required to evaluate the velocity field of (4.10), and C R denotes the cost associated with all other operations. Assume that the cost to evaluate the Hamiltonian at the reduced solution has order O(α(N )). Then, a standard algorithm for the evaluation of the right hand side of (4.10) will have arithmetic complexity , where the last two terms are associated with the computation of Y Z, and the inversion of C + J 2n C J 2n , respectively. This Section focuses on the development of structure-preserving numerical methods for the solution of (4.10) such that C R is at most linear in N . The efficient treatment of the nonlinear terms is out of the scope of the present study and will be the subject of future investigations on structure-preserving hyper-reduction techniques.
To simplify the notation, we recast (4.10) as: For Q ∈ M, find U ∈ C 1 (T , R 2N ×2n ) such that where, for any fixed t ∈ T , where H U is defined as in (4.7), and T U M = {V ∈ R 2N ×2n : U V ∈ g 2n }. In a temporal splitting perspective, we assume that the matrix Z (t) ∈ V p×2n is given at each time instant t ∈ T . Owing to Proposition 4.5, if Q ∈ M, then U (t) ∈ M for all t ∈ T . Then, the goal is to develop an efficient numerical scheme such that the discretization of (5.1) yields an approximate flow map with trajectories belonging to M. We propose two intrinsic numerical methods for the solution of the differential equation

Cayley transform as coordinate map
Orthosymplectic square matrices form a subgroup U(2N ) of a quadratic Lie group. We can, therefore, use the Cayley transform to induce a local parameterization of the Lie group U(2N ) near the identity, with the corresponding Lie algebra as parameter space. The following results extend to orthosymplectic matrices the properties of the Cayley transform presented in e.g. [ Then, (i) cay maps the Lie algebra g 2N into the Lie group G 2N .
(ii) cay is a diffeomorphism in a neighborhood of the zero matrix 0 ∈ g 2N . The differential of cay at ∈ g 2N is the map dcay : and its inverse is  The factor 1/2 in the definition (5.3) of the Cayley transform is arbitrary and has been introduced to guarantee that dcay 0 = I 2N , which will be used in Sect. 5.3 for the construction of retraction maps.
To derive computationally efficient numerical schemes for the solution of the basis evolution equation (5.1) we exploit the properties of analytic functions evaluated at the product of rectangular matrices. Proof Since has rank k it admits the splitting = αβ for some α, β ∈ R 2N ×k . To evaluate the Cayley transform in a computationally efficient way we exploit the properties of analytic functions of low-rank matrices. More in details, let f (z) := z −1 (cay(z) − 1) for any z ∈ C. The function f has a removable pole at z = 0. Its analytic extension reads, The cost to compute A := β α ∈ R k×k is O (N k 2 ). Moreover, The approach suggested hitherto is clearly not unique. The invertibility of the matrix A is ensured under the condition that the low-rank factors α and β are full rank. Although a low-rank decomposition with full rank factors is achievable [8,Proposition 4], one could alternatively envision the use of Woodbury matrix identity [33] to compute the matrix inverse appearing in the definition (5.3) of the Cayley transform. This yields the formula which can also be evaluated in O(Nrk) + O(k 2 r ) + O(k 3 ) operations.

Numerical integrators based on Lie groups acting on manifolds
In this Section we propose a numerical scheme for the solution of (5.1) based on Lie group methods, cf. [18]. The idea is to consider M as a manifold acted upon by the Lie group G 2N = U(2N ) of square orthosymplectic matrices. Then, since the local structure in a neighbourhood of any point of G 2N can be described by the corresponding Lie algebra g 2N , a local coordinate map is employed to derive a differential equation on g 2N . Since Lie algebras are linear spaces, using Runge-Kutta methods to solve the equation on g 2N allows to derive discrete trajectories that remain on the Lie algebra. This approach falls within the class of numerical integration schemes based on canonical coordinates of the first kind, also known as Runge-Kutta Munthe-Kaas (RK-MK) methods [24][25][26][27].
for any U ∈ M, then, Proof Let us consider, at each time t ∈ T , an orthosymplectic extension and YẎ = UU + WẆ . Expressing A(Y ,Ẏ ) explicitly in terms of U and W , and using the evolution equation satisfied by U , yields Moreover, since A is skew-symmetric, it holdṡ Once we have recast (5.1) into the equivalent problem (5.6), the idea is to derive an evolution equation on the Lie algebra g 2N via a coordinate map. A coordinate map of the first kind is a smooth function ψ : g 2N → G 2N such that ψ(0) = Id ∈ G 2N and dψ 0 = Id , where dψ : g 2N × g 2N → g 2N is the right trivialized tangent of ψ defined as For sufficiently small t ≥ t 0 , the solution of (5.6) is given by for t ∈ T , (t 0 ) = 0.
(5.12) Problem (5.12) can be solved using traditional RK methods . Let (b i , a i, j ), for  i = 1, . . . , s and j = 1, . . . , s, be the coefficients of the Butcher tableau describing an s-stage explicit RK method. Then, the numerical approximation of (5.12) in the interval (t m , t m+1 ] is performed as in Algorithm 1. As anticipated in Sect. 5.1, we resort to the Cayley transform as coordinate map in Algorithm 1. The use of the Cayley transform in the solution of matrix differential equations on Lie groups was proposed in [12,17,22]. Analogously to [12,Theorem 5], it can be shown that the invertibility of cay and dcay is guaranteed if U (t) ∈ M solution of (5.6) satisfies −1 / ∈ t∈T σ (U (t)). Note that choosing a sufficiently small time step for the temporal integrator can prevent the numerical solution from having an eigenvalue close to −1, for some t ∈ T . Alternatively, restarting procedures of the Algorithm 1 can be implemented similarly to [12, pp. 323, 324].
The computational cost of Algorithm 1 with ψ = cay is assessed in the following result.
For i = 1, 1 m = 0 and, hence, 1 = L(U 1 m ) = γ 1 δ 1 owing to (5.13). Using definition (5.4), it holds where e i , f i ∈ R 2N ×4n are defined as Using Line 3 of Algorithm 1, the rank of i m can be bounded as In principle one can solve the evolution equation (5.12) on the Lie algebra g 2N using the matrix exponential as coordinate map instead of the Cayley transform, in the spirit of [26]. However, there is no significant gain in terms of computational cost, as shown in details in Appendix A.
Although the computational complexity of Algorithm 1 is linear in the full dimension, it presents a suboptimal dependence on the number s of stages of the RK scheme. However, in practical implementations, the computational complexity of Proposition 5.4 might prove to be pessimistic in s, and might be mitigated with techniques that exploit the structure of the operators involved.
In the following Section we improve the efficiency of the numerical approximation of (5.1) by developing a scheme which is structure-preserving and has a computational cost O(N n 2 s), namely only linear in the dimension N of the full model and in the number s of RK stages.

Tangent methods on the orthosymplectic matrix manifold
In this Section we derive a tangent method based on retraction maps for the numerical solution of the reduced basis evolution problem (5.1). The idea of tangent methods is presented in [9,Sect. 2] and consists in expressing any U (t) ∈ M in a neighborhood of a given Q ∈ M, via a smooth local map R Q : T Q M → M, as By the continuity of V and the local rigidity condition, the map dR Q is invertible for sufficiently small t (i.e., V (t) sufficiently close to 0 ∈ T Q M) and hencė This strategy allows to solve the ODE (5.15) on the tangent space TM, which is a linear space, with a standard temporal integrator and then recover the approximate solution on the manifold M via the retraction map as in (5.14). If the retraction map can be computed exactly, this approach yields, by construction, a structure-preserving discretization. The key issue here is to build a suitable smooth retraction R : TM → M such that its evaluation and the computation of the inverse of its tangent map can be performed exactly at a computational cost that depends only linearly on the dimension of the full model.
In order to locally solve the evolution problem (5.15) on the tangent space to the manifold M at a point Q ∈ M we follow a similar approach to the one proposed in [10] for the solution of differential equations on the Stiefel manifold. Observe that, for any Q ∈ M, the velocity field F(Q) in (5.2), which describes the flow of the reduced basis on the manifold M, belongs to the space H Q defined in (4.7). We thus construct a retraction R Q : H Q → M as composition of three functions: a linear map ϒ Q from the space H Q to the Lie algebra g 2N associated with the Lie group G 2N acting on the manifold M, the Cayley transform (5.3) as coordinate map from the Lie algebra to the Lie group and the group action : that we take to be the matrix multiplication. This is summarized in the diagram below, In more details, we take ϒ Q to be, for each Q ∈ M, the linear map ϒ Q : The space T Q M can be characterized as follows.

Proposition 5.5 Let Q ∈ M be arbitrary. Then, V ∈ T Q M if and only if
To prove that V ∈ T Q M, we verify that Q V ∈ g 2n . Using the orthogonality of Q, and the assumption Q ∈ sp(2n) results in with S ∈ Sym(2n) ∩ sp(2n) arbitrary. We first verify that Q ∈ sp(2n). Using the orthogonality of Q, the fact that V ∈ T Q M and S ∈ sp(2n) results in We then verify that, with the above definition of , the matrix ( Q − Q )Q = − Q Q coincides with V . Using the fact that S ∈ Sym(2n) and V ∈ T Q M yields We can therefore characterize the tangent space of the orthosymplectic matrix manifold as This suggests that the linear map ϒ Q can be defined as where the last equality follows by (5.16). Note that Q = d Then the map R Q : H Q → M defined for any V ∈ H Q as Note that the matrix S ∈ Sym(2n) ∩ sp(2n) in the definition of the retraction (5.18) is of the form Its choice affects the numerical performances of the algorithm for the computation of the retraction and its inverse tangent map, as pointed out in [10,Sect. 3].
In the following Subsections we propose a temporal discretization of (5.15) with an s-stage explicit Runge-Kutta method and show that the resulting algorithm has arithmetic complexity of order C F + O(N n 2 ) at every stage of the temporal solver.

Efficient computation of retraction and inverse tangent map
In the interval (t m , t m+1 ] the local evolution on the tangent space, corresponding to (5.15), readsV Let (b i , a i, j ) for i = 1, . . . , s and j = 1, . . . , i − 1 be the coefficients of the Butcher tableau describing the s-stage explicit Runge-Kutta method. Then the numerical approximation of (5.15)-(5.14) with U 0 := Q ∈ M and V 0 = 0 ∈ T Q M is given in Algorithm 2.
Other than the evaluation of the velocity field F at R U m (V ), the crucial points of Algorithm 2 in terms of computational cost, are the evaluation of the retraction and the computation of its inverse tangent map. If we assume that both operations can be performed with a computational cost of order O(N n 2 ), then Algorithm 2 has an overall arithmetic complexity of order O(N n 2 s) + C F s, where C F is the cost to compute F(U ) in (5.2) at any given U ∈ M.
Computation of the retraction. A standard algorithm to compute the retraction R Q (5.18) at the matrix V ∈ R 2N ×2n requires O(N 2 n) for the multiplication between cay(ϒ Q (V )) and Q, plus the computational cost to evaluate the Cayley transform at ϒ Q (V ) ∈ R 2N ×2N . However, for any V ∈ H Q , the matrix ϒ Q (V ) ∈ g 2N admits the low-rank splitting We can revert to the results of Proposition 5.2 (with k = 4n) so that the retraction (5.18) can be computed as with computational cost of order O(N n 2 ).
Computation of the inverse tangent map of the retraction. Let Q ∈ M and V ∈ H Q . Using the definition of retraction (5.18) we have Then, the tangent map dR Q reads Fixing the fiber on T H Q corresponding to V ∈ H Q results in where we have used the linearity of the map ϒ Q . Assume we know W ∈ H R Q (V ) . We want to compute V ∈ H Q such that It is possible to solve problem (5.20) with arithmetic complexity O(N n 2 ) by proceeding as in [10,Sect. 3.2.1]. Since, for our algorithm, the result of [10] can be extended to the case of arbitrary matrix S ∈ Sym(2n) ∩ sp(2n) in (5.18), we report the more general derivation in Appendix A. Note that, for S = 0 and explicit Euler scheme, the two Algorithms 1 and 2 are equivalent.

Convergence estimates for the tangent method
Since the retraction and its inverse tangent map in Algorithm 2 can be computed exactly, the smoothness properties of R allow to derive error estimates for the approximate reduced basis in terms of the numerical solution of the evolution problem (5.15) in the tangent space.
Using the definition of Cayley transform (5.3) we have, for ϒ Q (·) := ϒ Q (·)/2, Since ϒ Q is skew-symmetric I 2N − ϒ Q (V ) −1 is normal. Then, Hence, since Q and Y are (semi-)orthogonal matrices, it holds Using the definition of ϒ Q from (5.17) results in It follows that the solution of Algorithm 2 can be computed with the same order of accuracy of the RK temporal scheme.

Numerical experiment
To gauge the performances of the proposed method, we consider the numerical simulation of the finite-dimensional parametrized Hamiltonian system arising from the spatial approximation of the one-dimensional shallow water equations (SWE). The shallow water equations are used in oceanography to describe the kinematic behavior of thin inviscid single fluid layers flowing over a changing topography. Under the assumptions of irrotational flow and flat bottom topography, the fluid is described by the scalar potential φ and the height h of the free-surface, normalized by its mean value, via the nonlinear system of PDEs  (h 1 , . . . , h N , φ 1 , . . . , φ N ). The discrete set of parameters h is obtained by uniformly sampling with 10 samples per dimension, for a total of p = 100 different configurations. This implies that the full model variable R in (4.1) is the N × p matrix given by 2N , for any k ∈ {1, . . . , p}, where η k h denotes the k-th entry of the vector η h ∈ R p containing the samples of the parameters. We consider second order accurate centered finite difference schemes to discretize the first order spatial derivative in (6.1). The evolution problem (6.1) admits a canonical symplectic Hamiltonian. Spatial discretization with centered finite differences yields a Hamiltonian dynamical system where the Hamiltonian associated with the k-th parameter is given by Wave-type phenomena often exhibit a low-rank behavior only locally in time, and, hence, global (in time) model order reduction proves ineffective in these situations. We show this behavior by comparing the performances of our dynamical reduced basis method with the global symplectic reduced basis approach of [30] based on complex SVD. For the latter, a symplectic reduced space is obtained from the full model obtained by discretizing (6.1) with centred finite differences in space, with N = 1000, and the implicit midpoint rule in time, with t = 10 −3 . We consider snapshots every 10 time steps and 4 uniformly distributed samples of per dimension. Concerning the dynamical reduced model, we evaluate the initial condition (h h (0; η h ), φ h (0; η h )) at all values η h and compute the matrix R 0 ∈ R 2N × p having as columns each of the evaluations. As initial condition for the reduced system (4.10), we use U (0) = U 0 ∈ R 2N ×2n obtained via complex SVD of the matrix R 0 truncated at n, while Z (0) = U T 0 R 0 . Then, we solve system (4.10) with a 2-stage partitioned Runge-Kutta method obtained as follows: the evolution equation for the coefficients Z is discretized with the implicit midpoint rule; while the evolution equation (5.1) for the reduced basis is solved using the tangent method described in Algorithm 2 with the explicit midpoint scheme, i.e. s = 2, b 1 = 0, b 2 = 1, and a 1,1 = a 1,2 = a 2,2 = 0, a 2,1 = 1/2. Note that the resulting partitioned RK method has order of accuracy 2 and the numerical integrator for Z is symplectic [15,Sect. III.2]. Finally, the nonlinear quadratic operator in (6.1), is reduced by using tensorial techniques [32].
In Fig. 1 we report the error in the Frobenius norm, at final time, between the full model solution and the reduced solution obtained with the two different approaches and various dimensions of the reduced space. Note that the runtime includes also the offline phase for the global approach.
The results of Fig. 1 show that the dynamical reduced basis method outperforms the global approach by reaching comparable accuracy at a reduced computational cost. Moreover, as the dimension of the reduced space increases, the runtime of the global Fig. 2 Evolution of the error in the conservation of the Hamiltonian method becomes comparable to the one required to solve the high-fidelity problem, meaning that there is no gain in performing global model order reduction. Figure 2 shows the evolution of the error in the conservation of the discrete Hamiltonian, averaged over all p values of the parameter. Since the Hamiltonian is a cubic quantity, we do not expect exact conservation associated with the proposed partitioned RK scheme. In addition, as pointed out at the end of Sect. 4, we cannot guarantee exact preservation of the invariants at the interface between temporal intervals, since the reduced solution is projected into the space spanned by the updated basis. However, the preservation of the symplectic structure both in the reduction and in the discretization yields a good control on the Hamiltonian error, as it can be observed in Fig. 2.

Concluding remarks and future work
Nonlinear dynamical reduced basis methods for parameterized finite-dimensional Hamiltonian systems have been developed to mitigate the computational burden of large-scale, multi-query and long-time simulations. The proposed techniques provide an attractive computational approach to deal with the local low-rank nature of Hamiltonian dynamics while preserving the geometric structure of the phase space even at the discrete level.
Possible extensions of this work involve the numerical study of the proposed algorithm including high order splitting temporal integrators, numerical approximations ensuring the exact conservation of Hamiltonian, and restarting procedures of the Cayley RK algorithm. Moreover, the extension of dynamical reduced basis methods to Hamiltonian systems with a nonlinear Poisson structure would allow nonlinear structure-preserving model order reduction of a large class of problems, including Euler and Vlasov-Maxwell equations. Some of these topics will be investigated in forthcoming works.

A Exponential map
Let us consider Algorithm 1 with the exponential as coordinate map, namely ψ = exp : g 2N → g 2N . For any i m ∈ g 2N and U i m ∈ M, the matrix dexp −1 i m (L(U i m )) in Line 3 can be approximated by truncating the Baker-Campbell-Hausdorff (BCH) formula as where B k denotes the k-th Bernoulli number and ad 0 Observe that i ∈ R 2N ×2N belongs to the Lie algebra g 2N , and, hence, each i m is in g 2N and the solution of the RK-MK method remains on M, see e.g. [15,Theorem 8.4].
To assess the computational complexity of Algorithm 1, we need to consider two operations: the evaluation of i and the computation of U i m = exp( i m )U m , for any i = 2, . . . , s. To this aim, we first rewrite each commutator in (A.1) as a matrix polynomial.
Proof We proceed by induction on k.
The matrix polynomial form (A.2) allows to estimate the rank of the { i } s i=1 . Lemma A.2 Let i ∈ g 2N be defined as in (A.1). Then, Proof Using Lemma A.1, we have that, for any A, B ∈ g 2N , Since the rank of a matrix product is bounded by the minimum among the ranks of the factors, this implies that rank(C) ≤ (q + 1)rank(B) and, hence, rank( i ) ≤ (q + 1) rank(L(U i m )). We now prove that, for any 0 ≤ k ≤ q, there exists matrices E k , D k ∈ R 2N ×2N such that ad k A (B) = AE k + B D k . We proceed by induction on k. and, hence, rank(C) ≤ rank(A) + rank(B). This is equivalent to rank( i ) ≤ rank( i m ) + rank(L(U i m )). Using the definition of i m from Line 3 of Algorithm 1, the rank of i , for any i ≥ 2, can be bounded as rank( i ) ≤ rank(L(U i m )) + i−1 j=1 rank( j ).
Since rank( 1 ) = rank(L(U m )), it easily follows by induction that rank( i ) ≤ 2 i−1 rank(L(U i m )). Observe that the factorization (5.13) implies that each term {L(U i m )} s i=1 , with L defined in (5.8), has rank at most 4n. Therefore, from Lemma A.2, it follows that It can be inferred from (A.4) that the bound q + 1 is the one dominating in the computation of i m whenever the number s of RK stages is sufficiently large. An optimal number q opt of commutators to achieve the accuracy of the corresponding RK method can be derived as in [3,7]. We consider a few examples from [7, Table 3.1]: RKF45 has s = 6, q opt = 5, and hence q opt + 1 ≤ 2 i−1 for i ≥ 4; DVERK has s = 8, q opt = 10, and hence q opt + 1 ≤ 2 i−1 for i ≥ 5; Butcher7 has s = 9, q opt = 21, and hence q opt + 1 ≤ 2 i−1 for i ≥ 6. In light of these considerations, we consider the bound r i ≤ 4n(i − 1)(q + 1) although for small i this might not be sharp. With the estimate (A.4) on the rank of i m , we can assess the cost of computing U i m at each stage of the RK-MK Algorithm 1. The computation of the exponential of a matrix in g 2N requires O(N 3 ) operations, but this cost can be mitigated whenever the argument of the exponential is of low-rank. Similarly to Proposition 5.2, it can be shown that the cost to compute exp(uv )Y with u, v ∈ R 2N ×k and Y ∈ R 2N ×2n is O(k 3 + k 2 n + N nk) [8,Proposition 3]. In Algorithm 1 we need to evaluate the exponential of { i m } s i=1 , with rank( i m ) ≤ 4n(i − 1)(q + 1). Therefore, the computation of all U i m in Line 4, for 2 ≤ i ≤ s, requires O(N n 2 s 2 q + n 3 s 4 q 3 ).
The other contribution to the computational cost of Algorithm 1 comes from the evaluation of each i in (A.1). To estimate this cost, we resort to the polynomial expression (A.2) and the low-rank splitting of i m = α i β i , with α i , β i ∈ R 2N ×r i , and of L( The terms: Therefore, the overall computational cost to evaluate each i is O(Nr 2 i + q 2 r 3 i + nqr 2 i + N nr i ). Using r i = rank( i m ) ≤ 4n(i − 1)(q + 1) and summing over the stages of the RK scheme, all terms involved in Algorithm 1 can be evaluated with arithmetic complexity O(N n 2 q 2 s 3 +n 3 q 5 s 4 )+C F , where C F is the complexity of the algorithm to compute F(U ) in (5.2) at any given U ∈ M. The latter is, thus, the computational complexity of Algorithm 1 with ψ = exp. Since each i m can be written as the sum of elements in the Lie algebra g 2N , namely i m = i−1 j=1 α j β j with α j , β j ∈ R 2N ×r i , one might suggest to approximate exp( i m ) with E( i m ) := i−1 j=1 exp(α j β j ) in the spirit of [8]. However, such an approximation does not bring significant computational savings nor it is guaranteed to provide a good approximation of the exponential map.

A Efficient computation of the inverse tangent map
We propose an algorithm to solve (5.20) with a computational cost of order O (N n 2 ). We proceed exactly as in [10, Sect. 3.2.1] with the only difference that we consider any arbitrary S ∈ Sym(2n) ∩ sp(2n).
Using the definition of the derivative of the Cayley transform (5.4) we can recast (5.20) as Moreover, using the definition of R Q (V ) in (5.18) results in The skew-symmetric part of T 1 is then T 1 − T 1 = −(R Q (V ) Q + I 2n ) −1 (R Q (V ) + Q) T 2 . Therefore, It is straightforward to show that all operations involved in the computation of T 1 can be done with complexity of order O(N n 2 ).