Rank-adaptive structure-preserving model order reduction of Hamiltonian systems

This work proposes an adaptive structure-preserving model order reduction method for finite-dimensional parametrized Hamiltonian systems modeling non-dissipative phenomena. To overcome the slowly decaying Kolmogorov width typical of transport problems, the full model is approximated on local reduced spaces that are adapted in time using dynamical low-rank approximation techniques. The reduced dynamics is prescribed by approximating the symplectic projection of the Hamiltonian vector field in the tangent space to the local reduced space. This ensures that the canonical symplectic structure of the Hamiltonian dynamics is preserved during the reduction. In addition, accurate approximations with low-rank reduced solutions are obtained by allowing the dimension of the reduced space to change during the time evolution. Whenever the quality of the reduced solution, assessed via an error indicator, is not satisfactory, the reduced basis is augmented in the parameter direction that is worst approximated by the current basis. Extensive numerical tests involving wave interactions, nonlinear transport problems, and the Vlasov equation demonstrate the superior stability properties and considerable runtime speedups of the proposed method as compared to global and traditional reduced basis approaches.


Introduction
Hamiltonian systems describe conservative dynamics and non-dissipative phenomena in, for example, classical mechanics, transport problems, fluids and kinetic models.We consider finite-dimensional Hamiltonian systems, in canonical symplectic form, that depend on a set of parameters associated with the geometric configuration of the problem or which represent physical properties of the problem.The development of numerical methods for the solution of parametric Hamiltonian systems in many-query and long-time simulations is challenged by two major factors: the high computational cost required to achieve sufficiently accurate approximations, and the possible onset of numerical instabilities resulting from failing to satisfy the conservation laws underlying non-dissipative dynamics.Model order reduction (MOR) and reduced basis methods (RBM) provide an effective procedure to reduce the computational cost of such simulations by replacing the original high-dimensional problem with models of reduced dimensionality without compromising the accuracy of the approximation.The success of RBM relies on the assumption that the problem possesses a low-rank nature, i.e. that the set of solutions, obtained as time and parameters vary, is of low dimension.However, non-dissipative phenomena do not generally exhibit such global low-rank structure and are characterized by slowly decaying Kolmogorov n-widths.This implies that traditional reduced models derived via linear approximations are generally not effective.
In recent years, there has been a growing interest in the development of model order reduction techniques for transport-dominated problems to overcome the limitations of linear global approximations.A large class of methods consists in constructing nonlinear transformations of the solution manifold and to recast it in a coordinate framework where it admits a low-rank structure, e.g.[25,20,41,31,9,2,22,38].A second family of MOR techniques focuses on online adaptive methods that update local reduced spaces depending on parameter and time, e.g.[4,28,32].To the best of our knowledge, none of the aforementioned methods provides any guarantee on the preservation of the physical properties and the geometric structure of the problem considered, and they might therefore be unsuitable to treat non-dissipative phenomena.
In parametric dynamical systems, the state can be represented, at each time, as a matrix whose columns are the solution vectors associated with different parameter values.In this perspective, finding a low-dimensional space in which the solution state can be well approximated is strictly related to low-rank matrix approximations.In a time-dependent setting, dynamical low-rank approximation [21] provides a low-rank factorization updating technique to efficiently compute approximations of time-dependent large data matrices.This approach can be equivalently seen as a reduced basis method based on a modal decomposition of the approximate solution with dynamically evolving modes.A geometric perspective on the relation between dynamical low-rank approximation and model order reduction in the context of time-dependent matrices has been proposed in [10].To the best of our knowledge the only dynamical low-rank approximation methods able to preserve the geometric structure of Hamiltonian dynamics were proposed in [24] to deal with the spatial approximation of the stochastic wave equation and in [26] to deal with finite-dimensional Hamiltonian systems.The gist of these methods is to approximate the full model solution in a low-dimensional manifold that evolves in time and possesses the symplectic structure of the full phase-space.The reduced dynamics is then derived via a symplectic projection of the Hamiltonian vector field onto the tangent space of the reduced symplectic manifold at each reduced state.
Their success notwithstanding, traditional dynamical low-rank approximation techniques are based on a reduced (low-rank) space whose dimension is fixed at the beginning of the evolution.This is a major limitation since it frequently happens that the rank of the initial condition does not correctly reflect the effective rank of the solution at all times.Consider, as an example, a linear advection problem in 1D, where the parameter represents the transport velocity.It is clear that, if the initial condition does not depend on the parameter, its rank is equal to one.However, as the initial condition is advected in time with different velocities, its rank rapidly increases.Approximating such dynamics with a time-dependent sequence of reduced manifolds of rank-1 matrices yields poor approximations.Conversely, an overapproximation of the initial condition, and possibly of the solution at other times, could improve the accuracy but will inevitably yield situations of rank-deficiency, as observed in [21,Section 5.3].This example demonstrates that, in a dynamical reduced basis approach, it is crucial to accurately capture the rank of the full model solution at each time.This issue has, however, received little attention so far [34,7].In this work, we propose a novel dynamical low-rank approximation scheme for the solution of parametric Hamiltonian systems that combines adaptivity in the rank of the solution with preservation of the Hamiltonian structure of the dynamics.
The proposed rank-adaptive algorithm can be summarized as follows.
• Given a fixed partition of the temporal domain, we consider, in each temporal subinterval, the discretized reduced dynamical system obtained with the structure-preserving approach of [26].
While in [26] the rank of the approximate solution is fixed a priori, here we change the rank adaptively from one temporal interval to the next one.
• To this aim, a surrogate error based on a linearization of the problem residual is computed at chosen times and for all tested parameters.If the error indicator reveals, according to a specific criterion, that the current reduced space is too small to approximate the state, we augment it in the direction that is worst approximated by the current reduced basis.The reduced dynamical system is then evolved, in the subsequent temporal interval, in the augmented manifold.In case of overapproximation, the size of the reduced space is, instead, decreased.
• Two major difficulties are associated with this approach: (i) to maintain the global Hamiltonian structure of the dynamics while modifying the reduced phase space; and (ii) to evolve the system on the updated space starting from a rank-deficient initial condition.To address these problems, we devise a regularization of the velocity field of the reduced flow so that the resulting vector belongs to the tangent space of the updated reduced manifold, and the Hamiltonian structure is then preserved.
The remainder of the paper is organized as follows.In Section 2, we introduce parametrized Hamiltonian systems and describe their symplectic structure.In Section 3, we derive the evolution equations of the reduced Hamiltonian dynamics.The problem of overapproximation and rank-deficiency is discussed in Section 5, where the regularization algorithm is introduced.Section 4 deals with the numerical temporal integration of the reduced dynamics: first, we summarize the structure-preserving methods introduced in [26] for the evolution of the reduced basis, and then we design novel partitioned RK schemes that are accurate with order 2 and 3 and preserve the geometric structure of the evolution problem.Section 6 pertains to the rank-adaptive algorithm.We describe the major steps: computation of the error indicator, criterion for the rank update, and update of the reduced state.The computational complexity of the adaptive dynamical reduced basis algorithm is thoroughly analyzed in Section 7. Furthermore, an approach that combines tensorial and splitting techniques with coarsening strategies is proposed to efficiently deal with polynomial nonlinearities of the Hamiltonian gradient.Section 8 is devoted to extensive numerical simulations of the proposed algorithm and its numerical comparisons with non-adaptive and global reduced basis methods.Finally, Section 9 concludes with a few remarks.

Problem formulation
Let T := (t 0 , T ] be a temporal interval and let Γ ⊂ R d , with d ≥ 1, be a compact set of parameters.For each η ∈ Γ, we consider the Hamiltonian system described by the initial value problem: For where the dot denotes the derivative with respect to time t, V 2N is a 2N -dimensional vector space, and C 1 (T , V 2N ) denotes continuous differentiable functions in time taking values in V 2N .Moreover, the function H : V 2N × Γ → R is the Hamiltonian of the system, ∇ u is the gradient with respect to the state variable u, and J 2N is the so-called canonical symplectic tensor defined as with I N , 0 N ∈ R N ×N denoting the identity and zero matrices, respectively.The operator J 2N identifies a symplectic structure on the phase-space of the Hamiltonian system (2.1).Equivalently, the vector space V 2N admits a global basis that is symplectic and orthonormal according to the following definition.
Definition 2.1 (Orthosymplectic basis).The set of vectors {e i } 2N i=1 is said to be orthosymplectic in the 2N -dimensional vector space V 2N if e i J 2N e j = (J 2N ) i,j , and (e i , e j ) = δ i,j , ∀i, j = 1 . . ., 2N , where (•, •) is the Euclidean inner product and J 2N is the canonical symplectic tensor (2.2) on V 2N .

Dynamical reduced basis method for Hamiltonian systems
We are interested in solving the Hamiltonian system (2.1) for a given set of p vector-valued parameters {η j } p j=1 ⊂ Γ, that, with a small abuse of notation, we denote η h ∈ Γ h .Then, the state variable u in (2.1) can be thought of as a matrix-valued application u(•; η h ) : Throughout, for a given matrix R ∈ R 2N ×p , we denote with R j ∈ R 2N the vector corresponding to the j-th column of R, for any j = 1, . . ., p. Let [a 1 |a 2 | . . .|a r ] denote the matrix of size 2N × (m 1 + . . .+ m r ) resulting from the horizontal concatenation of the matrices a j ∈ R 2N ×mj for j = 1, . . ., r.The Hamiltonian system (2.1), evaluated at η h , can be recast as a set of ordinary differential equations in a 2N × p matrix unknown in V p 2N as follows.For R 0 (η h ) : where The function H j is the Hamiltonian of the dynamical system (2.1) corresponding to the parameter η j , for j = 1, . . ., p.We assume that, for a fixed sample of parameters η h ∈ Γ h , the vector field X H (•; η h ) ∈ V p 2N is Lipschitz continuous in the Frobenius norm • uniformly with respect to time, so that (3.1) is well-posed.
Let us consider the splitting of the time domain T into the union of intervals T τ := (t τ −1 , t τ ], τ = 1, . . ., N τ , with t 0 := t 0 and t Nτ := T , and let the local time step be defined as ∆t τ = t τ − t τ −1 for every τ .For the model order reduction of (3.1) we propose an adaptive dynamical scheme based on approximating the full model solution in a lower-dimensional space that is evolving, and whose dimension may also change over time.To this aim, we adopt a local perspective by considering, in each temporal interval, an approximation of the solution of (3.1) of the form where , . . ., p, and any t ∈ T τ .Here n τ ∈ N satisfies 2n τ ≤ p and n τ N , and is updated over time according to Algorithm 2 that we will thoroughly discuss in Section 6.With this notation, we introduce the collection of reduced spaces of 2N × p matrices having rank at most 2n τ , and characterized as where U represents the reduced basis and it is taken to be orthogonal and symplectic, while Z are the expansion coefficients in the reduced basis, i.e.
To approximate the Hamiltonian system (3.1) in T τ with an evolution problem on the reduced space M 2nτ we need to prescribe evolution equations for the reduced basis U (t) ∈ U τ and the expansion coefficients Z(t) ∈ Z τ .For this, we follow the approach proposed in [24] and [26], and derive the reduced flow describing the dynamics of the reduced state R in (3.2) by applying to the Hamiltonian vector field X H the symplectic projection Π T R(t) M2n τ onto the tangent space of the reduced manifold at the current state.The resulting local evolution problem reads: where we assume, for the time being, that the initial condition of (3.4) at time t τ −1 , τ ≥ 1, is given, and we refer to Section 6.3 for a complete description of how such an initial condition is prescribed.By exploiting the characterization of the projection operator Π T R(t) M2n τ in [26, Proposition 4.2], we obtain the local evolution equations for the factors U and Z in the modal decomposition of the reduced solution (3.2), as in [24, Proposition 6.9] and [26,Equation (4.10)].In more details, for any τ ≥ 1, given (U (t τ −1 ), where Y (t) := ∇H(R(t); η h ) ∈ V p 2nτ , and R(t) = U (t)Z(t) for all t ∈ T τ .Observe that the local expansion coefficients Z ∈ Z τ satisfy a Hamiltonian system (3.5a) of reduced dimension 2n τ , where the reduced Hamiltonian is defined as H U (Z; η h ) := H(U Z; η h ).
To compute the initial condition of the reduced problem at time t 0 we perform the complex SVD [29, Section 4.2] of R 0 (η h ) ∈ R 2N ×p in (3.1), truncated at the n 1 -th mode.Then, the initial reduced basis U 0 ∈ U 1 can be derived from the unitary matrix of left singular vectors of R 0 (η h ), via the isomorphism between U 1 and the Stiefel manifold of unitary N × n 1 complex matrices, cf.[24, Lemma 6.1].The expansion coefficients matrix is initialized as Z 0 = U 0 R 0 (η h ).

Partitioned Runge-Kutta methods
For the numerical time integration of the reduced dynamical system (3.4) we rely on partitioned Runge-Kutta (RK) methods.Partitioned RK methods were originally introduced to deal with stiff evolution problems by splitting the dynamics into a stiff and a nonstiff part so that the two subsystems could be treated with different temporal integrators.There are many other situations where a dynamical system possesses a natural partitioning, for example Hamiltonian or singularly perturbed problems, or nonlinear systems with a linear part.In our setting, the factorization of the reduced solution (3.2) into the basis U and the coefficients Z provides the natural splitting expressed in (3.5).
In this section we first consider structure-preserving numerical approximations of the evolution problems (3.5b) and (3.5a), treated separately.Then, for the numerical integration of the coupled system (3.5),we design partitioned RK schemes that are accurate with order 2 and 3 and preserve the geometric structure of each evolution problem.
Since the evolution equation (3.5a) is a Hamiltonian system (of reduced dimension) we can rely on symplectic methods for its temporal approximation, so that the symplectic properties of the flow are preserved at the discrete level, cf.e.g.[18].The evolution equation (3.5b) for the reduced basis is approximated using tangent methods analogous to the ones proposed in [26], and that we briefly summarize here.The idea of tangent methods for the solution of differential equations on manifolds, as introduced in [5], is to recast the local dynamics on the tangent space of the manifold, which is a linear space.The temporal approximation of (3.5b) by tangent methods allow to obtain, at a computational cost linear in N , a discrete reduced basis that is orthogonal and symplectic.Let F(•, •; η h ) : R 2N ×2nτ × Z τ → R 2N ×2nτ denote the velocity field of the evolution (3.5b) of the reduced basis, namely It can be easily shown that, for any Q ∈ U τ , F(Q, Z; η h ) belongs to the space This is a subspace of the tangent space of the manifold U τ of orthosymplectic 2N × 2n τ matrices at the point Q ∈ U τ .Let us assume to know, in each temporal interval T τ , the approximate solution . Then, any element of U τ , in a neighborhood of Q, can be expressed as the image of a vector V ∈ H Q via the retraction where cay is the Cayley transform, defined as cay(M ) = (I N −M/2) −1 (I N +M/2) for any skew-symmetric and Hamiltonian square matrix M ∈ R 2N ×2N .Since R Q is a retraction, rather than solving (3.5b) for U , one can derive the local behavior of U in a neighborhood of Q by evolving V (t), with U (t) = R Q (V (t)), in the space H Q .By computing the local inverse of the tangent map of the retraction R Q , the evolution problem for the vector V reads: for any t ∈ T τ , where We refer to [26,Section 5.3.1]for further details on the derivation of the function f τ .
The resulting set of evolution equations describes the reduced dynamics in each temporal interval T τ as: given where G := J 2n ∇H U (Z, η h ) from (3.5a) and f τ is defined in (4.4).
For the numerical approximation of (4.5), we derive partitioned Runge-Kutta methods.Let ) be the collection of coefficients of the Butcher tableau describing an s-stage symplectic RK method, and let ) be the set of coefficients of an s-stage explicit RK method.Then, the numerical approximation of (4.5) via partitioned RK integrators reads Runge-Kutta methods of order 2 and 3 with the aforementioned properties can be characterized in terms of the coefficients P Z and P U as in the following result.
Lemma 4.1.Consider the numerical approximation of (4.5) with the s-stage partitioned Runge-Kutta method (4.6) obtained by coupling the Runge-Kutta methods P Z = ({b i } s i=1 , {a ij } s i,j=1 ) and Then, the following statements hold.• Symplectic condition [18,Theorem VI.4.3].The Runge-Kutta method P Z is symplectic if • Order condition [17,Theorem II.2.13].The Runge-Kutta method P Z has order k, with Partitioned Runge-Kutta of order 2 and 3 can be derived as described in Appendix A.

Reduced dynamics under rank-deficiency
In Section 3 we have proposed to approximate the phase space of the full Hamiltonian system (3.1) by an evolving low-rank matrix manifold.Particular attention needs to be devoted to the case of overapproximation in which a full model solution with effective rank r < n is approximated by a rank-n matrix, as pointed out first in [21,Section 5.3].In this case, a rank-deficient reduced dynamical system needs to be solved and it is not clear how the effective rank of the reduced solution will evolve over time.Indeed, in each temporal interval T τ , the dynamics may not remain on the reduce manifold M 2nτ and the matrix S(Z) := ZZ + J 2nτ ZZ J 2nτ may become singular or severely ill conditioned.This happens, for example, when the full model state at time t 0 is approximated with a rank deficient matrix, or, as we will see in the rank-adaptive algorithm in Section 6, when the reduced solution at a fixed time is used as initial condition to evolve the reduced system on a manifold of states with increased rank.
In this section, we propose an algorithm to deal with the overapproximation while maintaining the geometric structure of the Hamiltonian dynamics and of the factors U and Z in (3.2).
Lemma 5.1 (Characterization of the matrix S).Let S := ZZ + J 2n ZZ J 2n ∈ R 2n×2n with Z ∈ R 2n×p and p ≥ 2n.S is symmetric positive semi-definite and it is skew-Hamiltonian, namely SJ 2n − J 2n S = 0.Moreover, if S has rank 2n then S is non-singular and S −1 is also skew-Hamiltonian.In particular, the null space of S is even dimensional and contains all pairs of vectors (v, J 2n v) ∈ R 2n × R 2n such that both v and J 2n v belong to the null space of Z .
Proof.It can be easily verified that S is symmetric positive semi-definite and skew-Hamiltonian.Any eigenvalue of a skew-Hamiltonian matrix has even multiplicity, hence the null space of S has even dimension.Since S is positive semi-definite, v ∈ ker(S) if and only if ZZ v = 0 and ZZ J 2n v = 0, that is ker(S) = ker(Z ) ∩ ker(Z J 2n ).Observe that all the elements v of the kernel of Z are such that J 2n v ∈ ker(Z J 2n ).
In addition to the algebraic limitations associated with the solution of a rank-deficient system, the fact that the matrix S might be singular or ill conditioned prevents the reduced basis from evolving on the manifold of the orthosymplectic matrices.As shown in [26,Proposition 4.3], if U (t τ −1 ) ∈ U τ then U (t) ∈ R 2N ×2nτ solution of (3.5b) in T τ satisfies U (t) ∈ U τ for all t ∈ T τ , owing to the fact that F(U, Z; η h ) belongs to the space H U in (4.2).

Lemma 5.2. The function F(•, •
and this is equal to X U J 2nτ if and only if J 2nτ S −1 = S −1 J 2nτ .This condition follows from Lemma 5.1.Lemma 5.2 can be equivalently stated by considering the velocity field F as a function of the triple (U, Z, S(Z)).Then F(U, Z, S(Z); η h ) belongs to H U if and only if U ∈ U τ , Z ∈ R 2nτ ×p and S(Z) is non-singular, symmetric and skew-Hamiltonian.If the matrix S is not invertible, i.e.Z / ∈ Z τ , its inverse needs to be replaced by some approximation S † .By Lemma 5.2, if S † is not symmetric skew-Hamiltonian, then F † (U, Z; η h ) := (I 2N − U U )AS † does no longer belong to the horizontal space H U .If, for example, S † is the pseudo inverse of S, then the above condition is theoretically satisfied, but in numerical computations only up to a small error, because, if S is rank-deficient, then its pseudoinverse corresponds to the pseudoinverse of the truncated SVD of S.
To overcome these issues in the numerical solution of the reduced dynamics (3.5), we introduce two approximations: first we replace the rank-deficient matrix S with an ε-regularization that preserves the skew-Hamiltonian structure of S and then, in finite precision arithmetic, we set as velocity field for the evolution of the reduced basis U an approximation of F in the space H U (t) , for all t ∈ T τ .The ε-regularization consists in diagonalizing S and then replacing, in the resulting diagonal factor, the elements below a certain threshold with a fixed factor ε ∈ R.This is possible since (real) symmetric matrices are always diagonalizable by orthogonal transformations.However, unitary transformations do not preserve the skew-Hamiltonian structure.We therefore consider the following Paige Van Loan (PVL) decomposition, based on symplectic equivalence transformations.

Lemma 5.3 ([40]
).Given a skew-Hamiltonian matrix S ∈ R 2n×2n there exists a symplectic orthogonal matrix W ∈ R 2n×2n such that W SW has the PVL form where S n ∈ R n×n is an upper Hessenberg matrix.
In our case, since the matrix S is symmetric, its PVL decomposition (5.1) yields tridiagonal matrices with identical blocks S nτ = S nτ .We further diagonalize S nτ using orthogonal transformations to obtain It can be easily verified that Q ∈ R 2nτ ×2nτ is orthogonal and symplectic.The PVL factorization Lemma 5.3 can be implemented as in, e.g., [1, Algorithms 1 and 2], with arithmetic complexity O(n 3 τ ).The factorization is based on orthogonal symplectic transformations obtained from Givens rotations [14] and symplectic Householder matrices, defined as the direct sum of Householder reflections [27].
Once the matrix S has been brought in the PVL form, we perform the ε-regularization.Introduce the diagonal matrix D nτ ,ε ∈ R nτ ×nτ defined as, and let us denote with D ε ∈ R 2nτ ×2nτ the diagonal matrix composed of two blocks, both equal to D nτ ,ε .The matrix S ε := QD ε Q ∈ R 2nτ ×2nτ is symmetric positive definite and skew-Hamiltonian.Its distance to S is bounded, in the Frobenius norm, as , where m ε is the number of elements of D nτ that are smaller than ε.Since the ε-regularized matrix S ε is invertible, S −1 ε exists and is skew-Hamiltonian.This property allows to construct the vector field with the property that F ε belongs to the tangent space of the orthosymplectic 2N × 2n τ matrix manifold.To gauge the error introduced by approximating the velocity field F in (4.1) with F ε , let us denote with L the operator Then, the error made in the evolution of the reduced basis (3.5b), by the ε-regularization, is Observe that the resulting vector field F ε belongs to the space H U by construction.However, in finite precision arithmetic, the distance of the computed F ε from H U might be affected by a small error that depends on the norm of the operators L and S ε .This rounding error can affect the symplecticity of the reduced basis over time, whenever the matrix S is severely ill conditioned.To guarantee that the evolution of the reduced basis computed in finite precision remains on the manifold of orthosymplectic matrices with an error of the order of machine precision, we introduce a correction of the velocity field It easily follows that, with either definitions, F ε, belongs to H U and the error in the Frobenius norm is We summarize the regularization scheme in Algorithm 1.

Rank-adaptivity
The dynamical reduced basis method that we have introduced in Section 3 is based on approximating the full model solution, in each temporal interval T τ , on a low-dimensional space of size n τ .The fact that the size of the reduced space can change over time allows to fully exploit the local low-rank nature of the solution.In this section, we propose an algorithm to detect when the reduced space needs to be enlarged or reduced and how this operation is performed.The method is summarized in Algorithm 2.
Here we focus on the case where the current rank of the reduced solution is too small to accurately reproduce the full model solution.In cases where the rank is too large, one can perform an ε-regularization following Algorithm 1 or decrease the rank by looking at the spectrum of the reduced state and remove the modes associated with the lowest singular values.

Error indicator
Error bounds for parabolic problems are long-established and have been widely used to certify global reduced basis methods, cf.e.g.[15,39].However, their extension to noncoercive problems often results in pessimistic bounds that cannot be used to properly assess the quality of the reduced approximation.Few works have focused on the development of error estimates (not bounds) for reduced solutions of advection-dominated problems.In this work, we propose an error indicator based on the linearized residual of the full model.A related approach, known as Dual-Weighted Residual method (DWR) [23], consists in deriving an estimate of the approximation error via the dual full model and the linearization of the error of a certain functional of interest (e.g.surface integral of the solution, stress, displacement, ...).Despite the promising results of this approach, the arbitrariness in the choice of the functional clashes with the goal of having a procedure as general as possible.
We begin with the continuous full model (3.1) and, for its time integration, we consider the implicit RK scheme used in the temporal discretization of the dynamical system for the expansion coefficients Z in (4.6), and having coefficients , The discrete residual operator, in the temporal interval T τ , is We consider the linearization of the residual operator (6.2) at (R τ , R τ −1 ), where R τ is the approximate reduced solution at time t τ , obtained from (4.6) as R τ = U τ Z τ ; thereby (6.3) Similar procedures have been adopted in the formulation of the piecewise linear methods for the approximation of nonlinear operators, providing accurate approximations in case of low-order nonlinearities.From the residual operator, an approximation of the local error R τ − R τ is given by the matrix-valued quantity E τ defined as with E 0 := R(t 0 ) − U 0 Z 0 .The quantity defined by (6.4) is the first order approximation of the error between the reduced and the full model solution.In particular, it quantifies the discrepancy due to the local approximation (3.2).Even if the linearization error is negligible, the computational cost related to the assembly of the entire full-order residual ρ and its Jacobian, together with the solution of a linear system for any instance of the p parameters η h , makes the indicator unappealing if used in the context of highly efficient reduced approximations.In [23], a hierarchical approach has been proposed to alleviate the aforementioned computational bottleneck but it relies on the offline phase to capture the dominant modes of the exact error.Instead, in this work, we solve (6.4) on a subset η h of the p vector-valued parameters η h of cardinality p p, and only each N E time steps during the simulation.To further reduce the computational cost, we compute (6.4) on a coarse mesh in the parameter domain, whenever possible, and then E τ is recovered on the original mesh via spline interpolation.Although the assembly and solution of the sparse linear system in (6.4) has, for example, arithmetic complexity O(N 1 2 ) [12] for problems originating from the discretization of two-dimensional PDEs, this sampling strategy allows to reduce the computational cost required by the error estimator as compared to the evolution of the reduced basis and the coefficients, as discussed in Section 8.

Criterion for rank update
Let E τ ∈ R 2N ×p be the error indicator matrix obtained in (6.4).To decide when to activate the rank update algorithm, we take into account that, for advection-dominated and hyperbolic problems discretized using spectral methods, the error accumulates, and the effect of unresolved modes on the resolved dynamic contributes to this accumulation [8].Moreover, it has been noticed [35] that, for many problems of practical interest, the modes associated with initially negligible singular values might become relevant over time, potentially causing a loss of accuracy if a reduced manifold of fixed dimension is employed.
Let us define t τ as the current time, t * as the last time at which the dimension of the reduced basis U was updated and let λ τ be the number of past updates at time t τ .At the beginning of the simulation t * = t 0 and λ 0 = 0.The rank update is performed if the ratio between the norms of error indicators at t τ and t * satisfies the criterion where r, c ∈ R are control parameters larger than 1.The ratio of the norms of the error indicator gives a qualitatively indication of how the error is increasing in time and (6.5) fixes a maximum acceptable growing slope.Deciding what represents an acceptable slope is a problem-dependent task but the numerical results in Section 8 show little sensitivity of the algorithm with respect to r and c.Moreover, the variable λ τ induces a frequent rank-update when n τ is small and vice versa when n τ is large, hence controlling both the efficiency and the accuracy of the updating algorithm.We postpone to future investigations greedy strategies for the selection of optimal control parameters.Note that other (combinations of) criteria are possible: one alternative is to check that the norm of the error indicator remains below a fixed threshold; another possibility is to control the norm of some approximate gradient of the error indicator, etc.By numerically testing these various criteria, we observe that, at least in the numerical simulations performed, the criterion (6.5) based on the ratio of error indicators is reliable and robust and gives the largest flexibility.

Update of the reduced state
If criterion (6.5) is satisfied, the rank adaptive algorithm updates the current reduced solution to a new state having a different rank.Specifically, assume that, in the time interval T τ −1 , we have solved the discrete reduced problem (4.6) to obtain the reduced solution As a first step, we derive an updated basis To this aim, we enlarge U τ −1 with two extra columns derived from an approximation of the error, analogously to a greedy strategy.In greater detail, with the algorithm described in Section 6.1, we derive the error matrix E τ associated with the reduced solution at the current time.Via a thin SVD, we extract the left singular vector associated with the principal component of the error matrix, and we normalize it in the 2-norm to obtain the vector e ∈ R 2N .We finally enlarge the basis U τ −1 with the two columns [e | J 2N e] ∈ R 2N ×2 .The rationale for this choice is that we seek to increase the accuracy of the low-rank approximation by adding to the reduced basis the direction that is worst approximated by the current reduced space.Numerical evidence of the improved quality of the updated basis in approximating the full model solution is provided in Section 8.1.
From the updated matrix , we construct an orthosymplectic basis in the sense of Definition 2.1, by performing a QR-like decomposition using symplectic unitary transformations.In particular, we employ a symplectic (modified) Gram-Schmidt algorithm [33], with the possibility of adding reorthogonalization [13] to enhance the stability and robustness of the algorithm.
Once the updated reduced basis U ∈ U τ is computed, we derive the matrix Z ∈ R 2nτ ×p by expanding the current reduced solution R τ −1 in the updated basis.Therefore, the updated Remark 6.1.Since the updated reduced state coincides with the reduced solution R τ −1 at time t τ −1 , all invariants of (3.1) preserved by the partitioned Runge-Kutta scheme (4.6) are conserved during the rank update.
Observe that, even if the current reduced state R τ −1 is in M 2nτ −2 , it does not belong to the manifold M 2nτ .Indeed, one easily shows that Z = U R τ −1 ∈ R 2nτ ×p does not satisfy the full-rank condition, As shown in Lemma 5.2, the fact that Z / ∈ Z τ implies that the velocity field F in (4.1), describing the evolution of the reduced basis, is not well-defined.Therefore, we need to introduce an approximate velocity field for the solution of the reduced problem (3.5) in the temporal interval T τ with initial conditions (U, Z) ∈ U τ × R 2nτ ×p .We refer to Section 5 for a discussion about this issue and the description of the algorithm designed to solve the rank-deficient reduced dynamics ensuing from the rank update.

Algorithm 2 Rank update
Compute the error indicator matrix E τ ∈ R 2N × p in (6.4)

Approximation properties of the rank-adaptive scheme
To gauge the local approximation properties of the rank-adaptive scheme for the solution of the reduced dynamical system (3.5),we consider the temporal interval T τ where the first rank update is performed.In other words, assume that is the numerical approximation of the solution R(t τ −1 ) ∈ M 2nτ−1 of the reduced dynamical system (3.4) at time t τ −1 with n τ −1 = n τ −2 = . . .= n 1 .After the rank update at time t τ −1 , the reduced state R satisfies the local evolution problem where (U nτ τ −1 , Z nτ τ −1 ) ∈ U τ × R 2nτ ×p are the rank-updated factors, and We make the assumption that the reduced problem (3.4) is well-posed.Let R(t) ∈ V p 2N be the full model solution of problem (3.1) in the temporal interval T τ with given initial condition R(t τ −1 ).The error between the approximate reduced solution of (6.6) and the full model solution at time t τ ∈ T is given by The quantity e τ A := R τ − R(t τ ) is the approximation error associated with the partitioned Runge-Kutta discretization scheme, and can be treated using standard convergence analysis techniques, in light of the fact that the retraction map is Lipschitz continuous in the Frobenius norm, as shown in [26,Proposition 5.7].The term e RA (t) := R(t) − R(t), for any t ∈ T τ , is associated with the rank update and can be bounded as , where L X H is the Lipschitz continuity constant of X H . Gronwall's inequality [16] gives, for all t ∈ T τ , Observe that the estimate (6.7) depends on the distance between the Hamiltonian vector field at the reduced state and its image under the map P ε R that approximates the orthogonal projection operator on the tangent space of M 2nτ .Although a rigorous bound for this term is not available, we expect that it can be controlled arbitrary well by increasing the size of the reduced basis, as will also be demonstrated in Section 8.Moreover, the estimate (6.7) on the whole temporal interval T depends exponentially on the final time T .A linear dependence on T can be obtained only in special cases, for example when ∇ R H is uniformly negative monotone.

Computational complexity of the rank-adaptive algorithm
In this section we discuss the computational cost required for the numerical solution of the reduced problem (3.5) with the rank-adaptive algorithm introduced in Section 6.
In each temporal interval T τ , the algorithm consists of two main steps: the evolution step, which entails the repeated evaluation of the velocity fields F and G in (4.6) at each stage of the Runge-Kutta temporal integrator, and the rank update step, which requires the evaluation of the error indicator and the update of the approximate reduced solution at the current time step.
The rank update strategy introduced in Section 6, and summarized in Algorithm 2, has an arithmetic complexity of O(N p 2 ) + O(N n 2 τ ) + O(N pn τ ), and the computational bottleneck is the computation of the error indicator.As suggested in Section 6.1, sub-sampling techniques and mesh coarsening can be employed to overcome this limitation.The evolution step consists in solving the discrete reduced system (4.6) in each temporal interval.To understand the computational complexity of this step, we neglect the number of nonlinear iterations required by the implicit temporal integrators for the evolution of the coefficients Z.The solution of (4.6) requires the evaluation of four operators: the velocity fields G and F, the retraction R and its inverse tangent map f τ .The algorithms proposed in [26, Section 5.3.1]for the computation of R and f τ have arithmetic complexity O(N n 2 τ ).We denote with C H = C H (N, n τ , p) the computational cost to evaluate the gradient of the reduced Hamiltonian at the reduced solution.Finally, the velocity field F is computed via Algorithm 1 with a computational complexity of O(N n τ p) + O(N n 2 τ ) + O(pn 2 τ ) + O(n 3 τ ), while C H is the cost to evaluate Y .It follows that the rank-adaptive algorithm for the solution of the reduced system (4.5) with a partitioned Runge-Kutta scheme has a computational complexity being at most linear in the dimension of the full model N , provided the computational cost C H to evaluate the Hamiltonian vector field at the reduced solution has a comparable cost.Concerning the latter, observe that the assembly of the reduced state R from the factors U and Z and the matrix-vector multiplication U ∇ R H(R; η h ) require O(N pn τ ) operations.Therefore, the computational bottleneck of the algorithm is associated with the evaluation of the Hamiltonian gradient at the reduced state R.
This problem is well-known in model order reduction and emerges whenever reduced models involve non-affine and nonlinear operators, cf.e.g.[30,Chapters 10 and 11].Several hyper-reduction techniques have been proposed to mitigate or overcome this limitation, resulting in approximations of nonlinear operators that can be evaluated at a cost independent of the size of the full model.However, we are not aware of any hyper-reduction method able to exactly preserve the Hamiltonian phase space structure during model reduction.Furthermore, hyper-reduction methods entail an offline phase to learn the low-rank structure of the nonlinear operators by means of snapshots of the full model solution.Compared to traditional global model order reduction, in a dynamical reduced basis approach the constraints on the computational complexity of the reduced operators is less severe since we allow the dimension of the full model to enter, albeit at most linearly, the computational cost of the operations involved.This means that the dynamical model order reduction can accommodate Hamiltonian gradients where each vector entry depends only on a few, say k N , components of the reduced solution, with a resulting computational cost of C H = O(N pn τ ) + O(kN p).This is the case when, for example, the dynamical system (2.1) ensues from a local discretization of a partial differential equation in Hamiltonian form.Note that this assumption is also required for the effective application of discrete empirical interpolation methods (DEIM) [6].
When dealing with low-order polynomial nonlinearities of the Hamiltonian vector field, we can use tensorial techniques to perform the most expensive operations only once and not at each instance of the parameter, as discussed in the following.

Efficient treatment of polynomial nonlinearities
Let us consider the explicit expression of the cost C H for different Hamiltonian functions H.If the Hamiltonian vector field X H in (3.1) is linear, then where A ∈ R 2N ×2N is a given linear application, associated with the spatial discretization of the Hamiltonian function H. Standard matrix-matrix multiplication to compute G has arithmetic complexity where k is the number of nonzero entries of the matrix A. The computational complexity of the algorithm is therefore still linear in N provided the matrix A is sparse.This is the case in applications we are interested in where the Hamiltonian system (3.1)ensues from a local spatial approximation of a partial differential equation.
In case of low-order polynomial nonlinearities, we use the tensorial representation [36] of the nonlinear function and rearrange the order of computing.The gist of this approach is to exploit the structure of the polynomial nonlinearities to separate the quantities that depend on the dimension of the full model from the reduced variables, by manipulating the order of computation of the various factors.Consider the evolution equations for the coefficients Z in (3.5a) for a single value η j of the parameter η h ∈ Γ h .The corresponding reduced Hamiltonian vector can be expressed in the form where Z j ∈ Z τ with p = 1, q ∈ N is the polynomial degree of the nonlinearity, A i ∈ R 2N ×2N are sparse discrete differential operators, G {q} represents the matricized q-order tensor and ⊗ denotes the Kronecker product.The last expression in (7.1) allows to separate the computations involving factors of size N from the reduced coefficients Z, so that the matrix G U ∈ R 2nτ ×(2nτ ) q can be precomputed during the offline phase.
In the case of the proposed dynamical reduced basis method, we employ the tensorial POD approach to reduce the computational complexity of the evaluation of G, the RHS of (3.5a), and its Jacobian needed in the implicit symplectic integrator at each time step of the numerical integrator.We start by noticing that a straightforward calculation of the second expression in (7.1) suggests O(cN pn τ ) + O(cpqk) + O(cN pq) operations, where the first term is due to the reduced basis ansatz and the Galerkin projection, the second term to the multiplication by the sparse matrices A i and the third term to the evaluation of a polynomial of degree q for each entry of a 2N × p matrix.The constant c represents the number of iterations of the Newton solver and k := max i k i , where k i is the number of nonzero entries of A i .Moreover, in each iteration we evaluate not only the nonlinear term but also its Jacobian, with an additional cost of O(cN p(q − 1)) + O(cpk G n τ ) + O(cN pn 2 τ ) operations, with k G being the number of nonzero entries of the full-order Jacobian.These terms represent, respectively, the operations required to evalute the polynomial functions in the Jacobian, the assembly of the Jacobian matrix and its Galerkin projection onto the reduced basis.This high computational cost can again be mitigated by resorting to the second formula in (7.1),where the term G U is precomputed at each iteration, for each stage of the partitioned RK integrator (4.6).To estimate the computational cost of the procedure we resort to the multi-index notation by introducing n := (n τ , . . ., n τ ) ∈ R n and hence G U Z in (7.1) can be recast as The arithmetic complexity of this step is O(qkn τ ) + O((q − 1)N n q τ ) + O(N n q+1 τ ), where the first term is due to the matrix multiplication of the q matrices A i U in (I), the second term to the pointwise and diagonal matrices multiplications involved in the computations of (II) and the third term to the multiplications by U T J 2N in (III).We stress that the cost required to assemble G U is independent of the number of parameters p and the number of iterations of the nonlinear solver.Once G U has been precomputed, the evaluation of the reduced RHS has a computational cost of O(cpn q+1 τ ) [36].The same splitting technique is exploited for each evaluation of the reduced Jacobian and most of the precomputed terms in (7.2) can be reused.The proposed treatment of polynomial nonlinearities results in an effective reduction of the computational cost in case of low-order polynomial nonlinearity (q = 2, 3), a large set of vector-valued parameters (p 10) and a moderate number n τ of basis vectors.

Numerical tests
To assess the performance of the proposed adaptive dynamical structure preserving reduced basis method, we consider finite-dimensional parametrized Hamiltonian dynamical systems arising from the spatial approximation of partial differential equations.Let Ω ⊂ R d be a continuous domain and let u : T × Ω × Γ → R m belong to a Sobolev space V endowed with the inner product •, • .A parametric evolutionary PDE in Hamiltonian form can be written as with suitable boundary conditions prescribed at the boundary ∂Ω.Here, the dot denotes the derivative with respect to time, and δ denotes the variational derivative of the Hamiltonian H defined as In the numerical tests, we consider, for any fixed value of the parameter η j ∈ Γ h , numerical spatial approximations of (8.1) that yield a 2N -dimensional Hamiltonian system in canonical form where u h belongs to a finite 2N -dimensional subspace of V, ∇ u is the gradient with respect to the state variable u h and H h : R 2N → R is such that ∆x 1 . . .∆x d H h is a suitable approximation of H. Testing (8.2) for p values η h = {η j } p j=1 of the parameter, yields a matrix-valued ODE of the form (3.1), where the j-th column of the unknown matrix R(t) ∈ R 2N ×p is equal to u h (t, η j ) for all j = 1, . . ., p.
We validate our adaptive dynamical reduced basis method on several representative Hamiltonian systems of the form (8.2), of increasing complexity, and compare the quality of the adaptive dynamical approach with a reduced model with a global basis.The proposed approach, including all the steps introduced in the previous sections, is summarized in Algorithm 3.For the global model, we consider the method proposed in [29,Section 4.2], where a reduced basis is built via a complex SVD of a suitable matrix of snapshots and the reduced model is derived via symplectic Galerkin projection onto the space spanned by the global basis.We analyze and compare the accuracy, conservation properties and efficiency of the reduced models by monitoring the various quantities.To assess the approximation properties of the reduced model, we track the error, in the Frobenius norm, between the full model solution R and the reduced solution R at any time t ∈ T , namely Moreover, we study the conservation of the Hamiltonian via the relative error in the 1 -norm in the parameter space Γ h , that is Finally, we monitor the computational cost of the different reduction strategies.Throughout, the runtime is defined as the sum of the lengths of the offline and online phases in the case of the complex SVD (global method); while, for the dynamical approaches it is the time required to evolve basis and coefficients (4.5) plus the time required to compute the error indicator and update the dimension of the approximating manifold, in the adaptive case.The adaptive dynamical reduced basis method is numerically tested on two nonlinear problems, the shallow water and Schrödinger equations in one and two dimensions.Finally, we consider a preliminary application to particle simulations of plasma physics problem with the reduction of the Vlasov equation with a forced external electric field, modeling the evolution of charged particle beams.All numerical simulations are performed using Matlab computing environment on computer nodes with Intel Xeon E5-2643 (3.40GHz).The code and the data supporting the findings of this study are available from the authors upon request.
Algorithm 3 Rank-adaptive reduced basis method Compute U 0 ∈ U 1 via complex SVD of R 0 (η h ) truncated at the n 1 -th mode, and Initialize the error indicator matrix • Use the tensorial POD approach (7.1) to assemble the operator G • Use the retraction map given in (4.3) to compute R Uτ−1 • Compute f τ according to (4.4), using Regularization (Algorithm 1), with parameter ε as input, to assemble F 6: if mod(τ, N E ) = 0 then 7: Compute the error indicator matrix and check the rank update criterion using Rank_update (Algorithm 2) as end if end for 10: end procedure

Shallow water equations
The shallow water equations (SWE) describe the kinematic behaviour of a thin inviscid single fluid layer flowing over a variable topography.In the setting of irrotational flows and flat bottom topography, the fluid is described by a scalar potential φ and the canonical Hamiltonian formulation (8.1) is recovered [37].The resulting time-dependent nonlinear system of PDEs is defined as with spatial coordinates x ∈ Ω, time t ∈ T , state variables h, φ : Ω × T → R, ∇• and ∇ divergence and gradient differential operators in x, respectively.The variable φ is the scalar potential of the fluid and h represents the height of the free-surface, normalized by its mean value.The system is coupled with periodic boundary conditions for both the state variables.The evolution problem (8.5) admits a canonical symplectic Hamiltonian form (8.1) with the Hamiltonian We consider numerical simulations in d = 1 and d = 2 dimensions on rectangular spatial domains.The domain Ω is partitioned using a Cartesian mesh in M − 1 equispaced intervals in each dimension, having mesh width ∆x and ∆y, when d = 2.As degrees of freedom of the problem we consider the nodal values of the height and potential, i.e. u h (t; η h ) := (h h , φ h ) = (h 1 , . . ., h N , φ 1 , . . ., φ N ), for all t ∈ T and η h ∈ Γ h , where N := M d , h m = h i,j with m := (j − 1)M + i, and i, j = 1, . . ., M .In 1D, N = M , and the index j is dropped.We consider second order accurate central finite difference schemes to discretize the differential operators in (8.5), and denote with D x and D y the discrete differential operators acting in the x-and y-direction, respectively.The semi-discrete formulation of (8.5) represents a canonical Hamiltonian system with the gradient of the Hamiltonian function with respect to u h given by where is the Hadamard product between two vectors.The discrete Hamiltonian is In the one-dimensional case, the operator D y vanishes.

One-dimensional shallow water equations (SWE-1D)
For this example, we set Ω = [−10, 10] and we consider the parameter domain Γ = 1 10 , 1 7 × 2 10 , 15  10 .The discrete set of parameters Γ h is obtained by uniformly sampling Γ with 10 samples per dimension, for a total of p = 100 different configurations.Problem (8.5) is completed with the initial condition with η h = (α, β), where α controls the amplitude of the initial hump in the depth h and β describes its width.We consider a partition of the spatial domain Ω into N − 1 equispaced intervals with N = 1001.
The full model solution u h (t; η h ) is computed using a uniform step size ∆t = 10 −3 in the time interval T = (0, T := 7].We use the implicit midpoint rule as time integrator because, being symplectic, it preserves the geometrical properties of the flow of the semi-discrete equation associated to (8.7).To study the reducibility properties of the problem, we explore the solution manifold and collect the solutions to the high-fidelity model in different matrices.The global snapshot matrix S ∈ R 2N ×(Nτ p) contains the snapshots associated with all sampled parameters η h and time steps, while, for any τ = 1, . . ., N τ , the matrix S τ ∈ R 2N ×p collects the full model solutions at fixed time t τ .In Figure 1a, we compare the normalized singular values of S and S τ , averaged over time for the latter.Although, in both cases, the exponential decay of the spectrum suggests the existence of reduced approximation spaces, the decay of the singular values of the averaged S τ is roughly 5 times faster than that of S. This difference suggests that a low-rank dynamical approach may be beneficial to reduce the computational cost and to increase the accuracy of the solution of the reduced model compared to a method with a global basis.Furthermore, the evolution of the numerical rank of S τ over time, reported in Figure 1b, shows a rapid growth during the first steps, followed by a mild increase in the remaining part of the simulation.This is compatible with the observations, made in Section 6.2, about the behavior of the singular value spectrum for advection dominated problems.
In order to compare the performances of local and global model order reduction, we consider, as global reduced method, the complex SVD approach [29] with reduced dimension 2n ∈ {10, 20, 30, 40, 60, 80}.This is used to generate a symplectic reduced basis from the solution of the high-fidelity model (8.5) obtained every 10 time steps and by uniformly sampling Γ with 4 samples per dimension.The reduced system is solved using the implicit midpoint rule with the same time step ∆t used for the full order model.quadratic operator, describing the evolution of (8.5), is reduced by using the approach described in Section 7.1 and the reduced operators are computed once during the offline stage.
Concerning the adaptive dynamical reduced model, we evaluate the initial condition (8.9) at all values η h ∈ Γ h and compute the matrix S 1 ∈ R 2N ×p having as columns each of the evaluations.As initial condition for the reduced system (3.5),we use where U 0 ∈ R 2N ×2n1 is obtained using the complex SVD applied to the snapshot matrix S 1 .System (3.5) is then evolved using the 2-stage partitioned Runge-Kutta method described in (A.1).For the following numerical experiments, we consider 2n 1 ∈ {6, 8, 10, 12} as initial dimensions of the approximating reduced manifolds.As control parameters for the rank update criterion of Algorithm 2, we fix the value c = 1.2 and study examples with r ∈ {1.02, 1.05, 1.1, 1.2}.Moreover, we examine the case in which the rank-updating algorithm is never triggered, i.e., the basis U (t) evolves in time but its dimension is fixed (n τ = n 1 for all τ ).In the adaptive case, the error indicator E τ in (6.4) is computed every 100 iterations using a coarse mesh with 500 equispaced intervals on the subset η h obtained by sampling 5 parameters per dimension from Γ h .In Figure 2, we compare the global reduced model, the dynamical models for different values of r, and the high-fidelity model in terms of total runtime and accuracy at the final time T by monitoring the error (8.3).The results show that, as we increase the dimension of the global reduced basis, the global reduced model provides accurate approximations but the runtime becomes larger than the one required to solve the high-fidelity problem.Hence, the global method loses the efficiency.The adaptive dynamical reduced approach outperforms the global reduced method by reaching comparable levels of accuracy at a computational time which is one order of magnitude smaller than the one required by the global reduction.Compared to the high-fidelity solver, the adaptive dynamical reduced method achieves an accuracy of E(T ) = 2.55 • 10 −5 with a speedup up of 42, in the best-case scenario.For this numerical experiment, the effectiveness of the rank update algorithm is limited by the error introduced in the approximation of the initial condition via a reduced basis.While the error is reduced from a factor of 4 in the case of 2n 1 = 8 to a factor of 20 in the case of 2n 1 = 12, compared to the non adaptive method, the accuracy is not significantly improved when 2n 1 = 6.We note that, when the adaptive algorithm is effective, the additional computational cost associated with the evaluation of the error indicator and the evolution of a larger basis is balanced by a considerable error reduction.3), at time T = 7, as a function of the for the complex SVD method ( ), the dynamical RB method ( ) and the adaptive dynamical RB method for different values of the control parameters r and c ( , ).For the sake of comparison, we report the runtime required by the high-fidelity solver ( ) to compute the numerical solutions for all values of the parameter η h ∈ Γ h .
To better gauge the accuracy properties of the adaptive dynamical reduced basis method, we compare its error with the error given by the high-fidelity solver for the same initial condition.The solution to the full model, with the projection of (8.9) onto the column space of U 0 as the initial condition, is the target of the adaptive reduced procedure, which aims at correctly representing the high-fidelity solution space at every time step.The importance of having a reduced space that accurately reproduces the initial condition can be inferred from Figure 3a: the error associated with a poorly resolved initial condition dominates over the remaining sources of error, and adapting the dimension of the reduced basis is not beneficial in terms of accuracy.As noted above, increasing 2n 1 not only improves the performance of the non adaptive reduced dynamical procedure but also boosts the potential gain, in terms of relative error reduction, of the adaptive method, as can be seen in Figure 3e.
In Figures 3 we report the growth of the dimension of the reduced basis for different initial dimension 2n 1 .For the evolution of the error, we do not notice any significant difference as the parameter r for the adaptive criterion (6.5) varies.Ideally, within each temporal interval, the reduced solution is close, in the Frobenius norm, to the best rank-2n τ approximation of the full model solution.To verify this property for the adaptive dynamical reduced basis method, we monitor the evolution of the error E ⊥ between the full model solution R, at the current time and for all η h ∈ Γ h , and its projection onto the space spanned by the current reduced basis evolved following (3.5b),namely E ⊥ (t) = R(t) − U (t)U (t) R(t) .
In Figure 4, the projection error is shown for different values of 2n 1 (Figures 4a and 4b) and the corresponding evolution of the reduced basis dimension is reported (Figures 4c and 4d).We notice that, when the dimension of the basis U is not adapted, the projection error tends to increase in time.This can be ascribed to the fact that the effective rank of the high-fidelity solution is growing and the reduced basis is no longer large enough to capture the rank-increasing solution.Adapting 2n τ during the simulation results in a zero-growth scenario, with local negative peaks when the basis is enlarged.This indicates that the strategy of enlarging the reduced manifold in the direction of the larger error (see Section 6.3) yields a considerable improvement of the approximation.In Figure 5 we show the relative error in the conservation of the Hamiltonian for different dimensions of the reduced manifold, and values of the control parameters r and c.As the Hamiltonian (8.8) is a cubic quantity, we do not expect exact conservation associated with the proposed partitioned Runge-Kutta temporal integrators.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 Figure 5.

Two-dimensional shallow water equations (SWE-2D)
We set Ω = [−4, 4] 2 as the spatial domain and Γ = 1 5 , 1 2 × 11 10 , 17 10 as the domain of parameters.We consider 10 uniformly spaced values of the parameter for each dimension of Γ to define the discrete subset Γ h .As initial condition, we consider The results presented in Figures 9 on the evolution of the error E(t) for 2n 1 = {4, 6, 8}, corroborate the conclusions, already drawn from the 1D test case, regarding the effect of a poorly approximated initial condition on the performances of the adapting procedure.The evolution of the basis dimension is reported in Figures 9b, 9d and 9f for different values of r, c and 2n 1 .

Two-dimensional nonlinear Schrödinger equation
The nonlinear Schrödinger equation (NLS) is used to model, among others, the propagation of light in nonlinear optical fibers and planar waveguides and to describe the Bose-Einstein condensates in a macroscopic gaseous superfluid wave-matter state at ultra-cold temperature.In the 2D setting, we test the adaptive strategy in the case of a Fourier mode cascade, where, starting from an initial condition represented by few low Fourier modes, the energy exchange to higher modes quickly complicates the dynamic of the problem [3].More specifically, in the spatial domain Ω, we consider the cubic Schrödinger equation with periodic boundary conditions, and vector-valued parameter η.By writing the complex-valued solution u in terms of its real and imaginary parts as u = q + iv, (8.12) can be written as a Hamiltonian system in canonical symplectic form with Hamiltonian Let us consider the spatial domain Ω = [−2π, 2π] 2 and the set of parameters Γ = [0.97,1.03] 2 .We seek the numerical solution to (8.12), for p = 64 uniformly sampled parameters η h := (α, β) ∈ Γ h entering the initial condition u 0 (x, y; η h ) = (1 + α sin x) (2 + β sin y) .(8.13)This problem is characterized by an energy exchange between Fourier modes.Although this process is local, it is not well understood how the energy exchange mechanism is influenced by the problem dimension and parameters.In particular, although the values of α and β have a limited impact on the low-rank structure of the initial condition (8.13), the explicit effect of their variation on the energy exchange process is not known.We use a centered finite difference scheme to discretize the Laplacian operator.The domain Ω is partitioned using M = 101 nodes per dimension, for a total of N = 10000 points and with ∆x = ∆y = 4π • 10 −2 .Let u h (t; η h ), for all t ∈ T and η h ∈ Γ h , be the vector collecting the degrees of freedom associated with the nodal approximation of u.The semi-discrete problem is canonically Hamiltonian with the discrete Hamiltonian function with periodic boundary conditions for q i,j and v i,j .We consider N τ = 12000 time steps in the interval T = (0, T := 3] so that ∆t = 2.5 • 10 −4 .As in the previous examples, the implicit midpoint rule is used as numerical in the high-fidelity solver.The reduced dynamics (3.5) is integrated using the 2-stage partitioned RK method.
To assess the reducibility of the problem, we collect in S ∈ R 2N ×(Nτ p) the snapshots associated with all parameters η h and times t τ , and in S τ ∈ R 2N ×p the snapshots associated with all parameters η h at fixed time t τ , with τ = 1, . . ., N τ .The slow decay of the singular values of S, reported in Figure 14a, suggests that a global reduced basis approach is not viable for model order reduction.The growing complexity of the high-fidelity solution, associated with different values of α and β, is reflected by the growth of the numerical rank shown in Figure 14b.Hence, despite the exponential decay of the singular values of S τ , Figure 14b indicates that this test represents a challenging problem even for the adaptive algorithm and a balance between accuracy and computational cost is necessary while adapting the dimension of the reduced manifold.

Vlasov-Poisson plasma model with forced electric field
The Vlasov-Poisson system describes the dynamics of a collisionless magnetized plasma under the action of a self-consistent electric field.The evolution of the plasma at any time t ∈ T ⊂ R is described in terms of the distribution function f s (t, x, v) (s denotes the particle species) in the Cartesian phase space domain (x, v) ∈ Ω := Ω x × Ω v ⊂ R 2 .In this work, we consider the one-species (s = 1) paraxial approximation of the Vlasov-Poisson equation, used in the study of long and thin beams of particles [19].More specifically, we assume that the beam has reached a stationary state, the longitudinal length of the beam is the predominant spatial scale and the velocity along the longitudinal direction is constant.Moreover, we look at the case in which the effects of the self-consistent electric field E are negligible compared to the ones caused by an external electric field that we denote by Ξ.The external electric field is assumed to be independent of time and periodic with respect to the longitudinal dimension.Using the scaling argument proposed in [11] and the aforementioned assumptions, the problem is: For f 0 ∈ V |t=0 , find f ∈ C 1 (T ; L 2 (Ω)) ∩ C 0 (T ; V ) such that where the electric field Ξ is prescribed at all t ∈ T , x ∈ Ω x , the parameter ν ∈ R represents a spatial scaling and the Vlasov equation has been normalized so that mass and charge are set to m = q = 1.In (8.14), since we are considering stationary states, the variable t can be interpreted as the longitudinal coordinate and T as the longitudinal spatial domain.
For the semi-discrete approximation of the Vlasov equation in (8.14) we consider a particle method: The distribution function f is approximated by the superposition of P ∈ N computational macro-particles each having a weight ω , so that where X(t) and V (t) are the vector of the position and velocity of the macro-particles, respectively, and S is a compactly supported shape function, here chosen to be the Dirac delta.The idea of particle methods is to derive the time evolution of the approximate distribution function f h by advancing the macro-particles along the characteristics of the Vlasov equation.Particle methods, like particle-in-cell (PIC), are widely use in the numerical simulation of plasma problems.However, the slow convergence requires the use of many particles to achieve sufficient accuracy and therefore PIC methods are expensive.Model order reduction, in the number of macro-particles, of these semi-discrete schemes can be crucial and potentially extremely beneficial.
The particle approximation of problem (8.14) yields a Hamiltonian system where the unknowns are the vectors of position X and velocity V of the particles with the discrete Hamiltonian reads Here φ denotes the potential, defined as Ξ(x) = −∂ x φ(x), for all x ∈ Ω x , W p := diag(ω 1 , . . ., ω P ), and diag(d) denotes the diagonal matrix with diagonal elements given by the vector d.
For this test we consider N = 1000 particles with uniform unitary weight, ω i = 1, for all i = 1, . . ., N .The external electric field is given as Ξ(t, x) = −x 3 for all t ∈ T and x ∈ Ω x .The entries of the initial position X(0) and velocity V (0) vectors are independently sampled from the perturbed Maxwellian

Concluding remarks
We have considered parametrized non-dissipative problems in their canonical symplectic Hamiltonian formulation.For their model order reduction, we propose a nonlinear structure-preserving reduced basis method consisting in approximating the problem solution with a modal decomposition where both the expansion coefficients and the reduced basis are evolving in time.Moreover, the dimension of the reduced basis is updated in time according to an adaptive strategy based on an error indicator.The resulting reduced models allow to achieve stable and accurate results with small reduced basis even for problems characterized by a slowly decaying Kolmogorov n-width.The strength is the combination of the dynamical adaptivity of the reduced basis and the preservation of the geometric structure underlying key physical properties of the dynamics, illustrated by examples.
The study of efficient and structure-preserving algorithms for general nonlinear Hamiltonian vector fields and the development of partitioned Runge-Kutta methods that ensure the exact preservation of (at least linear and quadratic) invariants are still open problems and provide interesting directions of investigation.Moreover, the application of our rank-adaptive reduced basis method to fully kinetic plasma models, like the Vlasov-Poisson problem, might also be subject of future studies.P U the explicit midpoint method, namely the non-zero coefficients have values b 2 = b 2 = 1, a 22 = a 21 = 1/2.The resulting partitioned RK method has order of accuracy 2 and the numerical integrator P Z is symplectic.
To derive a partitioned Runge-Kutta method of order 3, we take P Z to be the 2-stage Gauss-Legendre (GL) method of order 4 enlarged with a fictitious stage.Starting from the enlarged 4-stage GL scheme in Table 1 (left), we derive an explicit RK method of order 3 by imposing the conditions (4.9) and (4.11).The resulting scheme is described by the Butcher tableau in Table 1 on the right.By construction and in view of Lemma 4.1, the following result holds.
Lemma A.2.The 3-stage partitioned Runge-Kutta method characterized by the set of coefficients P Z and P U in Table 1 has order of accuracy 3 and the numerical integrator P Z is symplectic.Table 1: Butcher tableau for the Gauss-Legendre scheme of order 4, on the left, and for the explicit 3-stage Runge-Kutta method of order 3, on the right.
We construct a partitioned RK scheme of order 3 with a larger region of absolute stability by including a further stage.This can be obtained by coupling the Gauss-Legendre scheme of order 6, suitably enlarged with a fictitious stage, to an explicit RK method, as in Table 2. Table 2: Butcher tableau for the Gauss-Legendre scheme of order 6, on the left and explicit 4-stage Runge-Kutta method of order 3, on the right.
The nine unknown coefficients of the explicit third order scheme in Table 2 are obtained by solving the underdetermined system derived by imposing the eight order conditions (4.9) and (4.11).A further equation can be imposed by adding a constraint on the region of absolute stability of the scheme: this is given by {z ∈ C : |R(z)| < 1} where the stability function is R(z) = 1 + z + z 2 /2 + z 3 /6 + Kz 4 , with K := b 4 a 43 a 32 a 21 .

Figure 1 :
Figure 1: SWE-1D: (a) Singular values of the global matrix S and time average of the singular values of the local trajectories matrix Sτ .The singular values are normalized using the largest singular value for each case.(b) -rank of the local trajectories matrix Sτ for different values of .

Figure 2 :
Figure2: SWE-1D: Error(8.3), at time T = 7, as a function of the for the complex SVD method ( ), the dynamical RB method ( ) and the adaptive dynamical RB method for different values of the control parameters r and c ( , ).For the sake of comparison, we report the runtime required by the high-fidelity solver ( ) to compute the numerical solutions for all values of the parameter η h ∈ Γ h .

Figure 3 :
Figure 3: SWE-1D: On the left column, we report the evolution of the E(t) (8.3) for the adaptive and non adaptive dynamical RB methods for different values of the control parameter r and different dimensions 2n1 of the approximating manifold of the initial condition.The target error is obtained by solving the full model with initial condition obtained by projecting (8.9) onto a symplectic manifold of dimension 2n1.On the right column, we report the evolution of the dimension of the dynamical reduced basis over time.The adaptive algorithm is driven by the error indicator (6.5), while in the non adaptive setting, the dimension does not change with time.We consider the cases 2n1 = 6 (Figs.(a)-(b)), 2n1 = 8 (Figs.(c)-(d)), 2n1 = 10 (Figs. (e)-(f)).

Figure 4 :
Figure 4: SWE-1D: In Figs.(a) and (b), we report the evolution of the projection error E ⊥ (t) for different values of the initial dimension 2n1 of the reduced manifold.Figs.(c) and (d), we report the corresponding evolution of the dimension of the reduced manifolds.

Figure 9 :
Figure 9: SWE-2D: On the left column, we report the evolution of the error E(t) (8.3) for the adaptive and non adaptive dynamical RB methods for values the control parameters r and c, and for different dimensions 2n1 of the initial reduced manifold.The target error is obtained by solving the full model with initial condition obtained by projecting (8.11) onto a symplectic manifold of dimension 2n1.On the right column, we report the evolution of the dimension of the dynamical reduced basis over time.The adaptive algorithm is driven by the error indicator (6.5), while in the non adaptive setting, the dimension does not change with time.We consider the cases 2n1 = 4 (Figs.(a)-(b)), 2n1 = 6 (Figs.(c)-(d)) and 2n1 = 8 (Figs.(e)-(f)).

Figure 10 :
Figure 10: NLS-2D: (a) Singular values of the global snapshots matrix S and of the time average of the local trajectories matrix Sτ .The singular values are normalized using the largest singular value for each case.(b) -rank of the local trajectories matrix Sτ for different values of .

Figure 17 :
Figure 17: Vlasov 1D1V: On the left, we show the evolution of the error E(t) (8.3) for the adaptive and non adaptive dynamical RB methods for different values of the control parameters r and c, and for different dimensions 2n1 of the initial reduced manifold.On the right, we show the evolution of the dimension of the dynamical reduced basis.We consider the cases 2n1 = 4 (Figs.(a)-(b)), 2n1 = 8 (Figs.(c)-(d)) and 2n1 = 12 (Figs.(e)-(f)).