Symplectic Model Order Reduction with Non-Orthonormal Bases
Patrick Buchfink Ashish Bhatt Bernard Haasdonk





Attributes-aided Part Detection and Refinement for Person Re-identification
    Shuzhao Li,
 Huimin Yu,
 Wei Huang,
 Jing Zhang 

 Zhejiang University 

 {leeshuz, yhm2005, huangwayne28, zj9301}@zju.edu.cn

    Received: 28 September 2018 / Accepted: 5 February 2019
==================================================================================================================================




Parametric high-fidelity simulations are of interest for a wide range of applications. But the restriction of computational resources renders such models to be inapplicable in a real-time context or in multi-query scenarios. Model order reduction (MOR) is used to tackle this issue. Recently, MOR is extended to preserve specific structures of the model throughout the reduction, e.g. structure-preserving MOR for Hamiltonian systems. This is referred to as symplectic MOR. It is based on the classical projection-based MOR and uses a symplectic reduced order basis (ROB). Such a ROB can be derived in a data-driven manner with the Proper Symplectic Decomposition (PSD) in the form of a minimization problem. Due to the strong nonlinearity of the minimization problem, it is unclear how to efficiently find a global optimum. In our paper, we show that current solution procedures almost exclusively yield suboptimal solutions by restricting to orthonormal ROBs. As new methodological contribution, we propose a new method which eliminates this restriction by generating non-orthonormal ROBs. In the numerical experiments, we examine the different techniques for a classical linear elasticity problem and observe that the non-orthonormal technique proposed in this paper shows superior results with respect to the error introduced by the reduction.






§ INTRODUCTION



Simulations enable researchers of all fields to run virtual experiments that are too expensive or impossible to be carried out in the real world. In many contexts, high-fidelity models are indispensable to represent the simulated process accurately. These high-fidelity simulations typically come with the burden of large computational cost such that an application in real-time or an evaluation for many different parameters is impossible respecting the given restrictions of computational resources at hand. Model order reduction (MOR) techniques can be used to reduce the computational cost of evaluations of the high-fidelity model by approximating these with a surrogate reduced-order model (ROM) <cit.>.

One class of high-fidelity models are systems of ordinary differential equations (ODEs) with a high order, i.e. a high dimension in the unknown variable. Such models typically arise from fine discretizations of time-dependent partial differential equations (PDEs). Since each point in the discretization requires one or multiple unknowns, fine discretizations with many discretization points yield a system of ODEs with a high order. In some cases, the ODE system takes the form of a finite-dimensional Hamiltonian system. Examples are linear elastic models <cit.> or gyro systems <cit.>.

Symplectic MOR <cit.> allows to derive a ROM for high-dimensional Hamiltonian systems by lowering the order of the system while maintaining the Hamiltonian structure. Thus, it is also referred to as structure-preserving MOR for Hamiltonian systems <cit.>. Technically speaking, a Petrov–Galerkin projection is used in combination with a symplectic reduced-order basis (ROB).

For a data-driven generation of the ROB, the conventional methods e.g. the Proper Orthogonal Decomposition (POD) <cit.> are not suited since they do not necessarily compute a symplectic ROB. To this end, the referenced works introduce the Proper Symplectic Decomposition (PSD) which is a data-driven basis generation technique for symplectic ROBs. Due to the high nonlineariy of the optimization problem, an efficient solution strategy is yet unknown for the PSD. The existing PSD methods (Cotangent Lift, Complex SVD, a nonlinear programming approach <cit.> and a greedy procedure introduced in <cit.>) each restrict to a specific subset of symplectic ROBs from which they select optimal solutions which might be globally suboptimal.

The present paper classifies the existing symplectic basis generation techniques in two classes of methods which either generate orthonormal or non-orthonormal bases. To this end, we show that the existing basis generation techniques for symplectic bases almost exclusively restrict to orthonormal bases. Furthermore, we prove that Complex SVD is the optimal solution of the PSD on the set of orthonormal, symplectic bases. During the proof, an alternative formulation of the Complex SVD for symplectic matrices is introduced. To leave the class of orthonormal, symplectic bases, we propose a new basis generation technique, namely the PSD SVD-like decomposition. It is based on an SVD-like decomposition of arbitrary matrices ∈^n × 2m introduced in <cit.>.

This paper is organized in the following way: <Ref> is devoted to the structure-preserving MOR for autonomous and non-autonomous, parametric Hamiltonian systems and thus, introduces symplectic geometry, Hamiltonian systems and symplectic MOR successively. The data-driven generation of a symplectic ROB with PSD is discussed in <Ref>. The numerical results are presented and elaborated in <Ref> exemplified by a Lamé–Navier type elasticity model which we introduce at the beginning of that section together with a short comment on the software that is used for the experiments. The paper is summarized and concluded in <ref>.




§ SYMPLECTIC MODEL REDUCTION


Symplectic MOR for autonomous Hamiltonian systems is introduced in <cit.>. We repeat the essentials for the sake of completeness and to provide a deeper understanding of the methods used. In the following ∈⊂^p describe p ∈ parameters of the system from the parameter set . We might skip the explicit dependence on the parameter vector  if it is not relevant in this specific context.



 §.§ Symplectic geometry in finite dimensions


[Symplectic form over ]
  Let  be a finite-dimensional vector space over . We consider a skew-symmetric and non-degenerate bilinear form : ×→ , i.e. for all _1,_2 ∈, it holds

    _1_2 = -_2_1        and        _2_3 = 0   ∀_3 ∈_3 = .

The bilinear form  is called symplectic form on  and the pair (, ) is called symplectic vector space.


It can be shown that  is necessarily of even dimension <cit.>. Thus,  is isomorphic to ^2n which is why we restrict to = ^2n and write [2n] instead of  in the following. In context of the theory of Hamiltonians, ^2n refers to the phase space which consists, in the context of classical mechanics, of position states = q_1, …, q_n∈^n of the configuration space and momentum states = p_1, …, p_n∈^n which form together the state = q_1, …, q_n, p_1, …, p_n∈^2n.

It is guaranteed <cit.> that there exists a basis _1, …, _n, _1, …, _n ⊂^2n such that the symplectic form takes the canonical structure

    [2n]_1_2 = _1 _2
      ∀_1, _2 ∈^2n,            := 
    [  n  n; -n  n ],

where n∈^n × n is the identity matrix, n∈^n × n is the matrix of all zeros and  is called Poisson matrix. Thus, we restrict to symplectic forms of the canonical structure in the following. For the Poisson matrix, it holds for any ∈^2n

    = 2n,         =  = -2n,         = 0.

These properties are intuitively understandable as the Poisson matrix is a 2n-dimensional, 90^∘ rotation matrix and the matrix -2n can be interpreted as a rotation by 180^∘ in this context.

[Symplectic map]
  Let A: ^2m→^2n, ↦, ∈^2n × 2m be a linear mapping for n,m ∈ and m≤ n. We call A a linear symplectic map and  a symplectic matrix with respect to [2n] and [2m] if
  
    = m.

  where [2m] is the canonical symplectic form on ^2m (and is equal to [2n] if n=m).
  
  Let ⊂^2m be an open set and :→^2n a differentiable map on . We call  a symplectic map if the Jacobian matrix () ∈^2n × 2m is a symplectic matrix for every ∈.


For a linear map, it is easy to check that the condition <ref> is equivalent to the preservation of the symplectic form, i.e. for all _1, _2 ∈^2m

    [2n]_1 _2 
    = _1 _2
    = _1 m_2
    = [2m]_1_2.


Now we give the definition of the so-called symplectic inverse which will be used in <Ref>.

[Symplectic inverse]
	For each symplectic matrix ∈^2n × 2m, we define the symplectic inverse
	
    = m∈^2m × 2n.


The symplectic inverse  exists for every symplectic matrix and it holds the following inverse relation

    = m = mm = 2m.




 §.§ Finite-dimensional, autonomous Hamiltonian systems


To begin with, we introduce the Hamiltonian system in a finite-dimensional, autonomous setting.

[Finite-dimensional, autonomous Hamiltonian system]
  Let :^2n×→ be a scalar-valued function that we require to be continuously differentiable in the first argument and which we call Hamiltonian (function). Hamilton's equation is an initial value problem with the prescribed initial data ∈, () ∈^2n which describes the evolution of the solution (t, ) ∈^2n for all t ∈,  ∈ with
  
    (t, ) =    ((t, ), ) =: ((t, ), ),      (, ) =    ()

  where (∙, ) is called Hamiltonian vector field. The triple (, [2n], ) is referred to as Hamiltonian system. We denote the flow of a Hamiltonian system as the mapping : ^2n×→^2n that evolves the initial state () ∈^2n to the corresponding solution (t, ; , ()) of Hamilton's equation

    (, ) := (t, ; , ()),

where (t, ; , ()) indicates that it is the solution with the initial data ,().


The two characteristic properties of Hamiltonian systems are (a) the preservation of the Hamiltonian function and (b) the symplecticity of the flow.

[Preservation of the Hamiltonian]
	The flow of Hamilton's equation  preserves the Hamiltonian function .



	We prove the assertion by showing that the evolution over time is constant for any ∈^2n due to
	
    (())
    = (())()
    <ref>=(())(())
    <ref>= 0.



[Symplecticity of the flow]
	Let the Hamiltonian function be twice continuously differentiable in the first argument. Then, the flow (∙, ): ^2n→^2n of a Hamiltonian system is a symplectic map.



	See <cit.>.




 §.§ Symplectic model order reduction for autonomous Hamiltonian systems


The goal of MOR <cit.> is to reduce the order, i.e. the dimension, of high dimensional systems. To this end, we approximate the high-dimensional state (t) ∈^2n with

    (t, ) ≈(t, ) = (t, ),         =

with the reduced state (t) ∈^2k, the reduced-order basis (ROB) ∈^2n × 2k, the reconstructed state (t) ∈ and the reduced space ⊂^2n. The restriction to even-dimensional spaces ^2n and ^2k is not necessary for MOR in general but is required for the symplectic MOR in the following. To achieve a computational advantage with MOR, the approximation should introduce a clear reduction of the order, i.e. 2k ≪ 2n.

For Petrov–Galerkin projection-based MOR techniques, the ROB  is accompanied by a projection matrix ∈^2n × 2k which is chosen to be biorthogonal to , i.e. = 2k. The reduced-order model (ROM) is derived with the requirement that the residual (t, ) vanishes in the space spanned by the columns of the projection matrix, i.e. in our case

    (t, ) = (t, ) - ((t, ), ) ∈^2n,        (t, ) = 2k × 1,

where 2k × 1∈^2k is the vector of all zeros. Due to the biorthogonality, this is equivalent to

    (t, ) = ((t, ), ) = ((t, ), ),        (, ) = ().


In the context of symplectic MOR, the ROB is chosen to be a symplectic matrix <ref> which we call a symplectic ROB. Additionally, the transposed projection matrix is the symplectic inverse = and the projection in <ref> is called a symplectic projection or symplectic Galerkin projection <cit.>. The (possibly oblique) projection reads

    =  =  = .

In combination, this choice of  and  guarantees that the Hamiltonian structure is preserved by the reduction which is shown in the following proposition.

[Reduced autonomous Hamiltonian system]
	Let  be a symplectic ROB with the projection matrix =. Then, the ROM <ref> of a high-dimensional Hamiltonian system (^2n, [2n], ) is a Hamiltonian system (^2k, [2k], ) on ^2k with the canonical symplectic form [2k] and the reduced Hamiltonian function (, ) = (, ) for all ∈^2k.



	First, we remark that the symplectic inverse is a valid biorthogonal projection matrix since it fulfils =  = 2k. To derive the Hamiltonian form of the ROM in <ref>, we use the identity
	
    = (<ref>)=
    = - 
    = ,

	which makes use of the properties <ref> of the Poisson matrix. It follows with <ref>
	
    (t)
    = ((t))
    = ((t))
    = ((t))

	where the last step follows from the chain rule. Thus, the evolution of the reduced state takes the form of Hamilton's equation and the resultant ROM is equal to the Hamiltonian system (^2k, [2k], ).


[Linear Hamiltonian system]
  Hamilton's equation is a linear system in the case of a quadratic Hamiltonian (, ) = 12 ()  + () with () ∈^2n × 2n symmetric and () ∈^2n
	
    (t, ) = () (t,) + (),        () = (),        () = ().

	The evolution of the reduced Hamiltonian system reads
	
    (t, ) = () (t,) + (),        () =    () <ref>=() ,
    () =    () <ref>=() ,
          () =    () ,
    () =    ().

	with the reduced Hamiltonian function (, ) = 12 ()  + ().



We emphasise that the reduction of linear Hamiltonian systems follows the pattern of the classical projection-based MOR approaches <cit.> to derive the reduced model with = and = which allows a straightforward implementation in existing frameworks.


Since the ROM is a Hamiltonian system, it preserves its Hamiltonian. Thus, it can be shown that the error in the Hamiltonian (t, ) = ((t, ), ) - ((t, ), ) is constant <cit.>. Furthermore, there are a couple of results for the preservation of stability <cit.>, <cit.> under certain assumptions on the Hamiltonian function.























[Offline/online decomposition]
A central concept in the field of MOR for parametric systems is the so-called offline/online decomposition. The idea is to split the procedure in a possibly costly offline phase and a cheap online phase where the terms costly and cheap refer to the computational cost. In the offline phase, the ROM is constructed. The online phase is supposed to evaluate the ROM fast. The ultimate goal is to avoid any computations that depend on the high dimension 2n in the online phase.

For a linear system, the offline/online decomposition can be achieved if (), () and () allow a parameter-separability condition <cit.>.

For systems with non-linear parts, multiple approaches <cit.> exist to enable an offline/online decomposition by introducing an approximation of the non-linear terms. This allows an online-efficient MOR of non-linear systems. For symplectic MOR, the symplectic discrete empirical interpolation method (SDEIM) was introduced <cit.> to preserve the symplectic structure throughout the approximation of the non-linear terms.




 §.§ Finite-dimensional, non-autonomous Hamiltonian systems


Non-autonomous Hamiltonian systems can be redirected to the case of autonomous systems if differentiability with respect to the time is assumed for the Hamiltonian function. The concept of the extended phase space is used. We briefly introduce the approach and explain the link to the symplectic MOR.

[Finite-dimensional, non-autonomous Hamiltonian system]
  Let :×^2n×→ be a scalar-valued function function that is continuously differentiable in the second argument. A non-autonomous (or time-dependent) Hamiltonian system (^2n, [2n], ) is of the form
  
    (t, ) = (t, (t, ), ).

  We therefore call (t, ) a time-dependent Hamiltonian function.


A problem for non-autonomous Hamiltonian systems occurs as the explicit time dependence of the Hamiltonian function introduces an additional variable, the time, and the carrier manifold becomes odd-dimensional. As mentioned in <ref>, symplectic vector spaces are always even-dimensional which is why a symplectic description is no longer possible. Different approaches are available to circumvent this issue.



As suggested in <cit.>, we use the methodology of the so-called symplectic extended phase space <cit.> to redirect the non-autonomous system to an autonomous system. The formulation is based on the extended Hamiltonian function : ^2n+2→ with

    () = (, ) + ,       =    ∈^2n + 2,       =  ∈^2n,      ,∈.

Technically, the time is added to the extended state  with = t and the corresponding momentum = -(t, (t)) is chosen such that the extended system is an autonomous Hamiltonian system.

This procedure requires the time-dependent Hamiltonian function to be differentiable in the time variable. Thus, it does for example not allow for the description of loads that are not differentiable in time in the context of mechanical systems. This might, e.g., exclude systems that model mechanical contact since loads that are not differentiable in time are required.



 §.§ Symplectic model order reduction of non-autonomous Hamiltonian systems



For the MOR of the, now autonomous, extended system, only the original phase space variable ∈^2n is reduced. The time and the corresponding conjugate momentum , are not reduced. To preserve the Hamiltonian structure, a symplectic ROB ∈^2n × 2k is used for the reduction of ∈^2n analogous to the autonomous case. The result is a reduced extended system which again can be written as a non-autonomous Hamiltonian system (^2k, [2k], ) with the time-dependent Hamiltonian (t,,) = (t,, ) for all (t, ) ∈ [, ] ×^2k.

An unpleasant side effect of the extended formulation is that the linear dependency on the additional state variable  (see <ref>) implies that the Hamiltonian cannot have strict extrema. Thus, the stability results listed in <cit.> and <cit.> do not apply if there is a true time-dependence in the Hamiltonian (t, ). Nevertheless, symplectic MOR in combination with a non-autonomous Hamiltonian system shows stable results in the numerical experiments.

Furthermore, it is important to note that only the extended Hamiltonian  is preserved throughout the reduction. The time-dependent Hamiltonian (·, t) is not necessarily preserved throughout the reduction, i.e. ((t)) = ((t)) but potentially ((t),t) ≠((t),t).



§ SYMPLECTIC BASIS GENERATION WITH THE PROPER SYMPLECTIC DECOMPOSITION (PSD)


We yet require a symplectic ROB for symplectic MOR. In the following, we pursue the approach of a ROB generated from a set of snapshots of the system. A snapshot is an element of the so-called solution manifold  that is approximated with a low-dimensional surrogate 

    := (t, )  |  t ∈, ∈⊂^2n,       := (t, )  |  t ∈, ∈≈.

In <cit.>, the Proper Symplectic Decomposition (PSD) is proposed as a snapshot-based basis generation technique for symplectic ROBs. The idea is to derive the ROB from a minimization problem which is suggested in analogy to the very well established Proper Orthogonal Decomposition (POD, also Principal Component Analysis) <cit.>. 

Classically, the POD chooses the ROB  to minimize the sum over squared norms of all ∈ residuals ( 2n -  ) _i of the orthogonal projection _i of the 1≤ i≤ single snapshots _i ∈ measured in the 2-norm ∙ with the constraint that the ROB  is orthogonal, i.e.

    ∈^2n × 2k∑_i=1^2n - _i^2      subject to   = 2k.

In contrast, the PSD requires the ROB to be symplectic instead of orthogonal which is expressed in the reformulated constraint. Furthermore, the orthogonal projection is replaced by the symplectic projection _i which results in

    ∈^2n × 2k∑_i=1^(2n - ) _i^2      subject to   = .

We summarize this in a more compact (matrix-based) formulation in the following definition.

[Proper Symplectic Decomposition (PSD)]
  Given  snapshots _1, …, _∈, we denote the snapshot matrix as = [_1, …, _] ∈^2n ×. Find a symplectic ROB ∈^2n × 2k which minimizes
	
    ∈^2n × 2k(2n - ) ^2      subject to   = ,

	We denote the minimization problem <ref> in the following as (), where  is the given snapshot matrix.


The constraint in <ref> ensures that the ROB  is symplectic and thus, guarantees the existence of the symplectic inverse . Furthermore, the matrix-based formulation <ref> is equivalent to the vector-based formulation presented in <ref> due to the properties of the Frobenius norm ∙.



 §.§ Symplectic, orthonormal basis generation



The foremost problem of the PSD is that there is no explicit solution procedure known so far due to the high nonlinearity and possibly multiple local optima. This is an essential difference to the POD as the POD allows to find a global minimum by solving an eigenvalue problem <cit.>.

Current solution procedures for the PSD restrict to a certain subset of symplectic matrices and derive an optimal solution for this subset which might be suboptimal in the class of symplectic matrices. In the following, we show that this subclass almost exclusively restricts to symplectic, orthonormal ROBs.

[Symplectic, orthonormal ROB]
  We call a ROB ∈^2n × 2k symplectic, orthonormal (also orthosymplectic, e.g. in <cit.>) if it is symplectic w.r.t. [2n] and [2k] and is orthonormal, i.e. the matrix  has orthonormal columns
  
    =       and       = 2k.



In the following, we show an alternative characterization of a symplectic and orthonormal ROB. Therefore, we extend the results given e.g. in <cit.> for square matrices ∈^2n × 2n in the following <ref> to the case of rectangular matrices ∈^2n × 2k.  This was also partially addressed in <cit.>.

[Characterization of a symplectic matrix with orthonormal columns]
The following statements are equivalent for any matrix ∈^2n × 2k

  
  *  is symplectic with orthonormal columns,
  
  *  is of the form
  
    =    =: ∈^2n × 2k,      ∈^2n × k,       = k,       = k,

  
  *  is symplectic and it holds =.



We remark that these matrices are characterized in <cit.> to be elements in Sp(2k, ^2n) ∩ V_k(^2n) where Sp(2k, ^2n) is the symplectic Stiefel manifold and V_k(^2n) is the Stiefel manifold.


  “(i)  (ii)”: Let ∈^2n × 2k be a symplectic matrix with orthonormal columns. We rename the columns to = [  ] with = [_1, …, _k] and = [_1, …, _k]. The symplecticity of the matrix written in terms of  and  reads
  
    = 
    [  ;   ]
    =
    [  k  k; -k  k ] =    = k,
    
    	- =    = k.

  Expressed in terms of the columns _i,_i of the matrices ,, this condition reads for any 1≤ i,j≤ k
  
    _i _j = 0,       _i _j = δ_ij,       _i _j = -δ_ij,       _i _j = 0,

  and the orthonormality of the columns of  implies
  
    _i _j = δ_ij,      _i _j = δ_ij.

  For a fixed i ∈{1, …, k}, it is easy to show with = 2n that _i is of unit length
  
    1 = δ_ii = _i _i = _i _i = _i^2.

  Thus, _i and _i are both unit vectors which fulfill _i _i = _i_i_^2n = 1. By the Cauchy-Schwarz inequality, it holds _i_i = _i_i if and only if the vectors are parallel. Thus, we infer _i = _i which is equivalent to _i = _i. Since this holds for all i ∈{1, …, k}, we conclude that  is of the form proposed in <ref>.

  “(ii)  (iii)”:
  Let  be of the form <ref>. Direct calculation yields
  
    = [ ;  ][   ]
    = [    ; -   ](<ref>)=[  k  k; -k  k ]
    =

  which shows that  is symplectic. Thus, the symplectic inverse  exists. The following calculation shows that it equals the transposed 
	
    = 
    	= [ ;  ]
    	= [ -;   ]
    	= [ -;   ]
    	= [ ;  ]
    	= .

	
	“(iii)  (i)”: Let  be symplectic with =. Then, we know that  has orthonormal columns since
	
    k =  = .



<Ref> essentially limits the symplectic, orthonormal ROB  to be of the form <ref>. Later in the current section, we see how to solve the PSD for ROBs of this type. In <Ref>, we are interested in ridding the ROB  of this requirement to explore further solution methods of the PSD.

As mentioned before, the current solution procedures for the PSD almost exclusively restrict to the class of symplectic, orthonormal ROBs introduced in <ref>. This includes the Cotangent Lift <cit.>, the Complex SVD <cit.>, partly the non-linear programming algorithm from <cit.> and the greedy procedure presented in <cit.>. We briefly review these approaches in the following proposition.

[Symplectic, orthonormal basis generation]
  The Cotangent Lift (CT), Complex SVD (cSVD) and the greedy procedure for symplectic basis generation all derive a symplectic and orthonormal ROB. The non-linear programming (NLP) admits a symplectic, orthonormal ROB if the coefficient matrix  in <cit.> is symplectic and has orthonormal columns, i.e. it is of the form _ = [  ]. The methods can be rewritten with = [  ], where the different formulations of  read

    _CT = [      ; n × k ]      _cSVD = [ ;  ],      _greedy = [_1, …, _k],      _NLP =

  where 
  
    
  * , , ∈^n × k are matrices that fulfil
	  
    = k,       +  = k,       = ,

	  which is technically equivalent to = k and = k (see <ref>) for _CT and _cSVD,
	  
  * _1, …, _k ∈^2n are the basis vectors selected by the greedy algorithm,
	  
  * ∈^2n × 2k is a ROB computed from CT or cSVD and ∈^2k × r, r ≤ k, stems from the coefficient matrix _ = [  ] computed by the NLP algorithm.
  



  All of the listed methods derive a symplectic ROB of the form = [  ] which satisfies <ref>. By <ref>, these ROBs are each a symplectic, orthonormal ROB. 


In the following, we show that PSD Complex SVD is the solution of the PSD in the subset of symplectic, orthonormal ROBs. This was partly shown in <cit.> which yet lacked the final step that, restricting to orthonormal, symplectic ROBs, a solution of ([   -]) solves () and vice versa. This proves that the PSD Complex SVD is not only near optimal in this set but indeed optimal. Furthermore, the proof we show is alternative to the original and naturally motivates an alternative formulation of the PSD Complex SVD which we call the POD of  in the following. To begin with, we reproduce the definition of PSD Complex SVD from <cit.>.

[PSD Complex SVD]
  We define the complex snapshot matrix
  
    = [_1 + _1, …, _ + _] ∈^n ×,      _j = [ _j; _j ]for all  1≤ j≤

  which is derived with the imaginary unit .
  The PSD Complex SVD is a basis generation technique that requires the auxiliary complex matrix ∈^n × k to fulfil
  
    ∈^n × k - ^2      subject to   = k

  and builds the actual ROB ∈^2n × 2k with
  
    = [  ],       = [ ;  ].

  The solution of <ref> is known to be based on the left-singular vectors of  which can be explicitly computed with a complex version of the SVD.


We emphasize that we denote this basis generation procedure as PSD Complex SVD in the following to avoid confusions with the usual complex SVD algorithm.

[Minimizing PSD in the set of symplectic, orthonormal ROBs]
  Given the snapshot matrix ∈^2n × we augment this with “rotated” snapshots to = [  ]. We assume that 2k is such that we obtain a gap in the singular values of , i.e. σ_2k() > σ_2k+1(). Then, minimizing the PSD in the set of symplectic, orthonormal ROBs is equivalent to the following minimization problem
  
    ∈^2n × 2k(2n - ) [   ]^2      subject to   = 2k.

  Clearly, this is equivalent to the POD <ref> applied to the snapshot matrix . We, thus, call this procedure the POD of  in the following. A minimizer can be derived with the SVD as it is common for POD <cit.>.



The proof proceeds in three steps: we show

  
  * that (,) is a pair of left- and right-singular vectors of  to the singular value σ if and only if (, ) also is a pair of left- and right-singular vectors of  to the same singular value σ,
  
  * that a solution of the POD of  is a symplectic, orthonormal ROB, i.e. =  = [  ],
  
  * that the POD of  is equivalent to the PSD for symplectic, orthonormal ROBs.

  We start with the first step (i). Let (,) be a pair of left- and right-singular vectors of  to the singular value σ. We use that the left-singular (or right-singular) vectors of  are a set of orthonormal eigenvectors of  (or ). To begin with, we compute
  
    =     ( + ) 
    =  + 
    = ,
    
    =    [  ;   ]
    = [   -; -   ]
    =

where we use = -. Thus, we can reformulate the eigenvalue problems of  and, respectively,  as
  
    σ
    =    
    =       ·|   σ
    =    <ref>=
    σ
    =    
    =       ·|   σ
    =    <ref>=.

  Thus, (, ) is necessarily another pair of left- and right-singular vectors of  with the same singular value σ. We infer that the left-singular vectors _i, 1≤ i≤ 2n, ordered by the magnitude of the singular values in a descending order can be written as

    = [_1   _1   _2   _2   …_n   _n] ∈^2n × 2n.

  
  For the second step (ii), we remark that the solution of the POD is explicitly known to be any matrix which stacks in its columns 2k left-singular vectors of the snapshot matrix  with the highest singular value <cit.>. Due to the special structure <ref> of the singular vectors for the snapshot matrix , a minimizer of the POD of  necessarily adopts this structure. We are allowed to rearrange the order of the columns in this matrix and thus, the result of the POD of  can always be rearranged to the form

    = [  ],       = [_1   _2   …  _k],       = [_1   _2   …  _k].

Note that it automatically holds that = k and () = k since, in both products, we use the left-singular vectors from the columns of the matrix  from <ref> which is known to be orthogonal from properties of the SVD. Thus, <ref> holds and we infer from <ref> that the POD of  indeed is solved by a symplectic, orthonormal ROB.
  
  For the final step (iii), we define the orthogonal projection operators
  
    =  =  + ,       = 2n - .

  Both are idempotent and symmetric, thus =  =. Due to = 2n, it further holds 
  
    =  =  -  -  =  = .

  Thus, it follows
  
    ^2
    = 
    = 
    = ^2

  and with = [  ]
  
    2 ^2 = ^2 + ^2
    =  [  ]^2
    = ^2,

  where we use in the last step that for two matrices ∈^2n × u, ∈^2n × v for u,v ∈, it holds ^2 + ^2 = [  ]^2 for the Frobenius norm ∙.
  
  Since it is equivalent to minimize a function f: ^2n × 2k→ or a multiple c f of it for any positive constant c ∈, minimizing ^2 is equivalent to minimizing 2^2 = ^2. Additionally, for a ROB of the form = [  ] the constraint of orthonormal columns is equivalent to the requirements in <ref>. Thus, to minimize the PSD in the class of symplectic, orthonormal ROBs is equivalent to the POD of  <ref>.


 
We remark that in the same fashion as the proof of step (iii) in <ref>, it can be shown that, restricting to symplectic, orthonormal ROBs, a solution of ([  ]) is a solution of () and vice versa, which is one detail that was missing in <cit.> to show the optimality of PSD Complex SVD in the set of symplectic, orthonormal ROBs.


We next prove that PSD Complex SVD is equivalent to POD of  from <ref> and thus, also minimizes the PSD in the set of symplectic, orthonormal bases. To this end, we repeat the optimality result from <cit.> and extend it with the results of the present paper.

[Optimality of PSD Complex SVD]

Let 𝕄_2 ⊂^2n × 2k denote the set of symplectic bases with the structure = [  ]. The PSD Complex SVD solves ([   -]) in 𝕄_2.



  See <cit.>.


[Equivalence of POD of  and PSD Complex SVD]
  PSD Complex SVD is equivalent to the POD of . Thus, PSD Complex SVD yields a minimizer of the PSD for symplectic, orthonormal ROBs.


  By <ref>, PSD Complex SVD minimizes <ref> in the set 𝕄_2 of symplectic bases with the structure = [  ]. Thus, <ref> holds with = which is equivalent to the conditions on  required in <ref>. By <ref>, we infer that 𝕄_2 equals the set of symplectic, orthonormal bases.
  
  Furthermore, we can show that, in the set 𝕄_2, a solution of ([   -]) is a solution of () and vice versa (see <ref>). Thus, PSD Complex SVD minimizes the PSD for the snapshot matrix  in the set of orthonormal, symplectic matrices and PSD Complex SVD and the POD of  solve the same minimization problem.


We emphasize that the computation of a minimizer of <ref> via PSD Complex SVD requires less memory storage than the computation via POD of . The reason is that the complex formulation uses the complex snapshot matrix ∈^n × which equals 2 · n · floating point numbers while the solution with the POD of  method artificially enlarges the snapshot matrix to ∈^2n × 2 which are 4 · n · floating point numbers. Still, the POD of  might be computationally more efficient since it is a purely real formulation and thereby does not require complex arithmetic operations.



 §.§ Symplectic, non-orthonormal basis generation


In the next step, we want to give an idea how to leave the class of symplectic, orthonormal ROBs. We call a basis generation technique symplectic, non-orthonormal if it is able to compute a symplectic, non-orthonormal basis.

In <ref>, we briefly showed that most existing symplectic basis generation techniques generate a symplectic, orthonormal ROB. The only exception is the NLP algorithm suggested in <cit.>. It is able to compute a non-orthonormal, symplectic ROB. The algorithm is based on a given initial guess _0 ∈^2n × 2k which is a symplectic ROB e.g. computed with PSD Cotangent Lift or PSD Complex SVD. Nonlinear programming is used to leave the class of symplectic, orthonormal ROBs and derive an optimized symplectic ROB = _0 with the symplectic coefficient matrix ∈^2k × 2r for some r ≤ k. Since this procedure searches a solution spanned by the columns of _0, it is not suited to compute a global optimum of the PSD which we are interested in the scope of this paper.

In the following, we present a new basis generation technique that is based on an SVD-like decomposition for matrices ∈^2n × m presented in <cit.>. To this end, we introduce this decomposition in the following.

[SVD-like decomposition <cit.>]
Any real matrix ∈^2n × m can be decomposed as the product of a symplectic matrix ∈^2n × 2n, a sparse and potentially non-diagonal matrix ∈^2n × m and an orthogonal matrix ∈^m × m with

    = ,       =  
      
        >c*3>c<
        >c
      
        p     q     p     m-2p-q 
    
          (>cccc<)
          >c
                     p 
                 q 
                 n-p-q 
                 p 
                 q 
                 n-p-q 
    ,
          = (_1, …, _p) ∈^p × p,
    _i > 0   for  1≤ i≤ p.

with p,q ∈ and () = 2p+q and where we indicate the block row and column dimensions in  by small letters. The diagonal entries _i, 1 ≤ i≤ p, of the matrix  are related to the pairs of purely imaginary eigenvalues λ_j(), λ_p+j() ∈ of = ∈^m × m with 

    λ_j() = -(_j)^2,      λ_p+j() = (_j)^2,      1 ≤ j≤ p.



[Singular values]
We call the diagonal entries _i, 1 ≤ i≤ p, of the matrix  from <ref> in the following the symplectic singular values. The reason is the following analogy to the classical SVD.

The classical SVD decomposes ∈^2n × m as = where ∈^2n × 2n, ∈^m × m are each orthogonal matrices and ∈^2n × m is a diagonal matrix with the singular values σ_i on its diagonal () = [σ_1, …, σ_r, 0, …, 0] ∈^min(2n,m), r = (). The singular values are linked to the real eigenvalues of = with λ_i() = σ_i^2. Furthermore, for the SVD, it holds due to the orthogonality of  and , respectively, = ^2 and = ^2.

A similar relation can be derived for an SVD-like decomposition from <ref>. Due to the structure of the decomposition <ref> and the symplecticity of , it holds

    =    ^=
    
    	  =    ,
    	    =
      
        >c*3>c<
        >c
      
        p     q     p     m-2p-q    
    
          (>cccc<)
          >c
              ^2        p
                 q
    
    			  -^2              p
                 m-2p-q
    .

This analogy is why we call the diagonal entries _i, 1 ≤ i≤ p, of the matrix  symplectic singular values.


The idea for the basis generation now is to select k ∈ pairs of columns of  in order to compute a symplectic ROB. The selection should be based on the importance of these pairs which we characterize by the following proposition by linking the Frobenius norm of a matrix with the symplectic singular values.


Let ∈^2n × m with an SVD-like decomposition = with p,q ∈ from <ref>. The Frobenius norm of  can be rewritten as

    ^2 = ∑_i=1^p+q (_i)^2,      _i = _i √(_i^2 + _n+i^2),     1≤ i≤ p,
    _i,     p+1 ≤ i≤ p+q

where _i ∈^2n is the i-th column of  for 1≤ i≤ 2n. In the following, we refer to each _i as the weighted symplectic singular value.



  We insert the SVD-like decomposition = and use the orthogonality of  to reformulate
  
    ^2 
    =    ^2
    = ^2
    = ()
    = ∑_i=1^p (_i)^2 _i _i + ∑_i=1^p (_i)^2 _n+i_n+i + ∑_i=1^q _p+i_p+i
    
      =    ∑_i=1^p (_i)^2 _i^2 + _n+i^2 
    + ∑_i=1^q _p+i^2

  which is equivalent to <ref>.


It proves true in the following <ref> that we can delete single addends _i in <ref> with the symplectic projection used in the PSD if we include the corresponding pair of columns in the ROB. This will be our selection criterion in the new basis generation technique that we denote PSD SVD-like decomposition.

[PSD SVD-like decomposition]
We compute an SVD-like decomposition <ref> as = of the snapshot matrix ∈^2n × and define p,q ∈ as in <ref>. In order to compute a ROB  with 2k columns, find the k indices i ∈ = {i_1, …, i_k}⊂{1, …, p+q} which have large contributions _i in <ref> with

    = _ℐ⊂{1, …, p+q}
    ℐ = k∑_i ∈ℐ_i^2 .

To construct the ROB, we choose the k pairs of columns _i ∈^2n from  corresponding to the selected indices  such that

    = [_i_1, …, _i_k, _n+i_1, …, _n+i_k] ∈^2n × 2k.



The special choice of the ROB is motivated by the following theoretical result which is very analogous to the results known for the classical POD in the framework of orthogonal projections.

[Projection error by neglegted weighted symplectic singular values]
Let ∈^2n × 2k be a ROB constructed with the procedure described in <ref> to the index set ⊂{1, …, p+q} with p,q ∈ from <ref>. The PSD functional can be calculated by

    (2n - ) ^2
    = ∑_i ∈{1, …, p+q}∖_i ^2,

which is the cumulative sum of the squares of the neglected weighted symplectic singular values.



  Let ∈^2n × 2k be a ROB constructed from an SVD-like decomposition = of the snapshot matrix ∈^2n × 2k with the procedure described in <ref>. Let p,q ∈ be defined as in <ref> and = {i_1, …, i_k}⊂{1, …, p+q} be the set of indices selected with <ref>.
  
  For the proof, we introduce a slightly different notation of the ROB . The selection of the columns _i of  is denoted with the selection matrix ∈^2n × 2k based on
  
    _α,β =
        
          1,   α=i_β∈
    
          0,   else   for      1≤α≤ 2n,
       1≤β≤ k,
            = [,  ].

  which allows us to write the ROB as the matrix–matrix product =. Furthermore, we can select the neglected entries with 2n -.

We insert the SVD-like decomposition and the representation of the ROB introduced in the previous paragraph in the PSD which reads
  
    (2n - ) ^2
    = (2n - ) ^2
    = ‖ (2n - ^= ) ‖^2_F

where we use the orthogonality of  and the symplecticity of  in the last step. We can reformulate the product of Poisson matrices and the selection matrix as

    = 
    [ ;  ]
    =
    [  k -k;  k  k ][  ; - ]
    = .

Thus, we can further reformulate the PSD as

    (2n - ) ^2
    = 2n - ^2
    = ∑_i ∈{1, …, p+q}∖_i ^2

where _i are the weighted symplectic singular values from <ref>. In the last step, we use that the resultant diagonal matrix in the braces sets all rows of  with indices i, n+i to zero for i ∈. Thus, the last step can be concluded analogously to the proof of <ref>.


A direct consequence of <ref> is that the decay of the PSD functional is proportional to the decay of the sum over the neglected weighted symplectic singular values _i from <ref>. In the numerical example <ref>, we observe an exponential decrease of this quantities which induces an exponential decay of the PSD functional.

[Computation of the SVD-like decomposition]
To compute an SVD-like decompostion <ref> of , several approaches exist. The original paper <cit.> derives a decomposition based on the product  which is not good for a numerical computation since errors can arise from cancellation. In <cit.>, an implicit version is presented that does not require the computation of the full product  but derives the decomposition implicitly by transforming . Furthermore, <cit.> introduces an iterative approach to compute an SVD-like decomposition which computes parts of an SVD-like decomposition with a block-power iterative method. In the present case, we use the implicit approach <cit.>.




 §.§ Interplay of non-orthonormal and orthonormal ROBs

We give further results on the interplay of non-orthonormal and orthonormal ROBs. The fundamental statement in the current section is the Orthogonal SR decomposition <cit.>.

[Orthogonal SR decomposition] 
For each matrix ∈^2n × m with m ≤ n, there exists a symplectic, orthogonal matrix ∈^2n × 2n, an upper triangular matrix _11∈^m × m and a strictly upper triangular matrix _21∈^m × m such that

    = [       _11; (n-m) × m;       _21; (n-m) × m ]
    = [_m   _m] [ _11; _21 ],      _i =     [_1, …, _i],    1≤ i≤ n,
     =     [_1, …, _n, _1, …, _n].


We remark that a similar result can be derived for the case m > n <cit.> but it is not introduced since we do not need it in the following.


Let ∈^2n × m with m ≤ n. We consider the QR decomposition

    = [           ; (2n-m) × m ]

where ∈^2n × 2n is an orthogonal matrix and ∈^2n × m is upper triangular. The original Orthogonal SR decomposition <cit.> for the square matrix states that we can decompose ∈^2n × 2n as a symplectic, orthogonal matrix ∈^2n × 2n, an upper triangular matrix _11∈^n × n, a strictly upper triangular matrix _21∈^n × n and two (possibly) full matrices _12, _22∈^n × n

    = [ _11 _12; _21 _22 ]      and thus      
    = [ _11 _12; _21 _22 ][           ; (2n-m) × m ]
    = [ _11; _21 ][          ; (n-m) × m ].

Since  is upper triangular, it does preserve the (strictly) upper triangular pattern in _11 and _21 and we obtain the (strictly) upper triangular matrices _11, _21∈^m × m from

    [       _11; (n-m) × m;       _21; (n-m) × m ]
    = [ _11; _21 ][          ; (n-m) × m ].



Based on the Orthogonal SR decomposition, the following two propositions prove bounds for the projection errors of PSD which allows an estimate for the quality of the respective method. In both cases we require the basis size to satisfy k ≤ n or 2k ≤ n, respectively. This restriction is not limiting in the context of symplectic MOR as in all application cases k ≪ n.


Let ∈^2n × k be a minimizer of POD with k ≤ n basis vectors and ∈^2n × 2k be a minimizer of the PSD in the class of orthonormal, symplectic matrices with 2k basis vectors. Then, the orthogonal projection errors of  and  satisfy

    (2n - ) ^2 ≤2n - ^2.




The Orthogonal SR decomposition (see <ref>) guarantees that a symplectic, orthogonal matrix ∈^2n × 2k and ∈^2k × k exist with =. Since both matrices  and  are orthogonal and () ⊂(), we can show that  yields a lower projection error than  with

    2n - ^2
      =    2n - 2n - ^2
      = ∑_i=1^2n - 2n - _i^2
    ≤    2n - ^2_≤ 1 ∑_i=1^2n - _i^2
      ≤2n - ^2


Let ∈^2n × 2k be a minimizer of the PSD in the class of symplectic, orthonormal ROBs. By definition of , it yields a lower projection error than . Since both ROBs are symplectic, orthonormal, we can exchange the symplectic inverse with the transposition (see <ref>, (iii)). This proves the assertion with

    2n - ^2
    ≥2n - ^2
    ≥2n - ^2.



<ref> proves that we require at most twice the number of basis vectors to generate a symplectic, orthonormal basis with an orthogonal projection error at least as small as the one of the classical POD. An analogous result can be derived in the framework of a symplectic projection which is proven in the following proposition.

 
Assume there exists a minimizer ∈^2n × 2k of the general PSD for a basis size 2k ≤ n with potentially non-orthonormal columns. Let ∈^2n × 4k be a minimizer of the PSD in the class of symplectic, orthogonal bases of size 4k. Then, we know that the symplectic projection error of  is less than or equal to the one of , i.e.
    (2n - ) ^2 ≤(2n - ) ^2.




  Let ∈^2n × 2k be a minimizer of PSD with 2k ≤ n. By <ref>, we can determine a symplectic, orthogonal matrix ∈^2n × 4k and ∈^4k × 2k with =. Similar to the proof of <ref>, we can bound the projection errors. We require the identity
  
    (2n - ) (2n - )
    = 2n -  -  + ^=4k_=_=
    = 2n - .

  With this identity, we proceed analogously to the proof of <ref> and derive for a minimizer ∈^2n × 4k of PSD in the class of symplectic, orthonormal ROBs
  
    (2n - ) ^2
        ≤    (2n - ) ^2
    = (2n - ) (2n - ) ^2
    ≤    (2n - )^2_≤ 1(2n - ) ^2
    ≤(2n - ) ^2.



<ref> proves that we require at most twice the number of basis vectors to generate a symplectic, orthonormal basis with a symplectic projection error at least as small as the one of a (potentially non-orthonormal) minimizer of PSD.



§ NUMERICAL RESULTS



The numerical experiments in the present paper are based on a two-dimensional plane strain linear elasticity model which is described by a Lamé–Navier equation

    t^2(, t, ) -  Δ_(, t, ) + ( + ) [] [](, t, ) =  (, t).

for ∈⊂^2 and t ∈ with the density ∈, the Lamé constants = (, ) ∈^2, the external force : ×→^2 and Dirichlet boundary conditions on ⊂ := ∂ and Neumann boundary conditions on ⊂. We apply non-dimensionalization (e.g. <cit.>), apply the Finite Element Method (FEM) with first-order Lagrangian elements on a triangular mesh and rewrite the system as first-order system to arrive at the parametric linear system <ref> with

    (t, ) = [ (t, ); (t, ) ],      () = [ ()  n;  n    ],      (t) = [  -(t); n × 1 ]

where (t, ) ∈^n is the vector of displacement DOFs, (t, ) ∈^n is the vector of linear momentum DOFs, () ∈^n × n is the stiffness matrix, ∈^n × n is the inverse of the mass matrix and (t, ) is the vector of external forces.

We remark that a Hamiltonian formulation with the velocity DOFs (t) = t(t) ∈^n instead of the linear momentum DOFs (t) is possible if a non-canonical symplectic structure is used. Nevertheless, in <cit.> it is suggested to switch to a formulation with a canonical symplectic structure for the MOR of Hamiltonian systems.

In order to solve the system <ref> numerically with a time-discrete approximation _i() ≈(t_i, ) for each of ∈ time steps t_i ∈, 1≤ i≤, a numerical integrator is required. The preservation of the symplectic structure in the time-discrete system requires a so-called symplectic integrator <cit.>. In the context of our work, the implicit midpoint scheme is used in all cases.

[Modified Hamiltonian]
We remark that, even though the symplectic structure is preserved by symplectic integrators, the Hamiltonian may be modified in the time-discrete system compared to the original Hamiltonian. In the case of a quadratic Hamiltonian (see <ref>) and a symplectic Runge-Kutta integrator, the modified Hamiltonian equals the original Hamiltonian since these integrators preserve quadratic first integrals. For further details, we refer to <cit.> or <cit.>.


The model parameters are the first and second Lamé constants with = (, ) ∈ = [35 · 10^9, 125 · 10^9]  ^2× [35 · 10^9, 83 · 10^9]  ^2 which varies between cast iron and steel with approx. 12% chromium <cit.>. The density is set to = 7856  ^3. The non-dimensionalization constants are set to =  = 81 · 10^9 ^2, ξ = 1, g = 9.81 ^2. The geometry is a simple cantilever beam clamped on the left side with a force applied to the right boundary. The time interval is chosen to be t ∈ [, ] with = 0 and = 7.2 · 10^-2 which is one oscillation of the beam. For the numerical integration = 151 time steps are used.



The symplectic MOR techniques examined are PSD Complex SVD (<ref>), the greedy procedure <cit.> and the newly introduced PSD SVD-like decomposition (<ref>). The MOR techniques that do not necessarily derive a symplectic ROB are called non-symplectic MOR techniques in the following. The non-symplectic MOR techniques investigated in the scope of our numerical results are the POD applied to the full state (t, ) (POD full state) and a POD applied to the displacement (t, ) and linear momentum states (t, ) separately (POD separate states). To summarize the basis generation methods, let us enlist them in <Ref> where (∙) and (∙) denote the SVD and the complex SVD, respectively.



All presented experiments are generalization experiments, i.e. we choose 9 different training parameter vectors ∈ on a regular grid to generate the snapshots and evaluate the reduced models for 16 random parameter vectors that are distinct from the 9 training parameter vectors. Thus, the number of snapshots is = 9 · 151 = 1359. The size 2k of the ROB  is varied in steps of 20 with 2k ∈ 20, 40, …, 280, 300.

The software used for the numerical experiments is RBmatlab[https://www.morepas.org/software/rbmatlab/] which is an open-source library based on the proprietary software package MATLAB® and contains several reduced simulation approaches. An add-on to RBmatlab is provided[https://doi.org/10.5281/zenodo.2578078] which includes all the additional code to reproduce the results of the present paper. The versions used in the present paper are RBmatlab 1.16.09 and MATLAB® 2017a.



 §.§ Autonomous beam model

In the first model, we load the beam on the free end (far right) with a constant force which induces an oscillation.  Due to the constant force, the discretized system can be formulated as an autonomous Hamiltonian system. Thus, the Hamiltonian is constant and its preservation in the reduced models can be analysed. All other reduction results are very similar to the non-autonomous case and thus, are exclusively presented for the non-autonomous case in the following <ref>.



  §.§.§ Preservation over time of the modified Hamiltonian in the reduced model


In the following, we investigate the preservation of the Hamiltonian of our reduced models. With respect to <ref>, we mean the preservation over time of the modified Hamiltonian. Since the Hamiltonian is quadratic in our example and the implicit midpoint is a symplectic Runge-Kutta integrator, the modified Hamiltonian equals the original which is why we speak of “the Hamiltonian” in the following.

We present in <ref> the count of the total 240 simulations which show a preservation (over time) of the reduced Hamiltonian in the reduced model. The solution  of a reduced simulation preserves the reduced Hamiltonian over time if ((t_i), ) - ((t_0), ) () < 10^-10 for all discrete times t_i ∈ [, ], 1≤ i≤ where () > 0 is a parameter-dependent normalization factor. The heat map shows that no simulation in the non-symplectic case preserves the Hamiltonian whereas the symplectic methods all preserve the Hamiltonian which is what was expected from theory.

In <ref>, we exemplify the non-constant evolution of the reduced Hamiltonian for three non-symplectic bases generated by POD separate states with different basis sizes and one selected test parameter (λ, μ) ∈. It shows that in all three cases, the Hamiltonian starts to grow exponentially.





 §.§ Non-autonomous beam model


The second model is similar to the first one. The only difference is that the force on the free (right) end of the beam is loaded with a time-varying force. The force is chosen to act in phase with the beam. The time dependence of the force necessarily requires a non-autonomous formulation which requires in the framework of the Hamiltonian formulation a time-dependent Hamiltonian function which we introduced in <ref>.

We use the model to investigate the quality of the reduction for the considered MOR techniques. To this end, we investigate the projection error, i.e. the error on the training data, the orthogonality and symplecticity of the ROB and the error in the reduced model for the test parameters.



  §.§.§ Projection error of the snapshots and singular values


The projection error is the error on the training data collected in the snapshot matrix , i.e.
    (2k) = (2n - ) ^2,   POD:  =    ,
    PSD:  =     (= for orthosymplectic ROBs, <ref>).

It is a measure for the approximation qualities of the ROB based on the training data. <ref> (left) shows this quantity for the considered MOR techniques and different ROB sizes 2k. All basis generation techniques show an exponential decay. As expected from theory, POD full state minimizes the projection error for the orthonormal basis generation techniques (see <ref>). PSD SVD-like decomposition shows a lower projection error than the other PSD methods for 2k ≥ 80 and yields a similar projection error for k ≤ 60. Concluding this experiment, one might expect the full-state POD to yield decent results or even the best results. The following experiments prove this expectation to be wrong.

The decay of (a) the classical singular values σ_i, (b) the symplectic singular values _i (see <ref>) and (c) the weighted symplectic singular values _i (see <ref>) sorted by the magnitude of the symplectic singular values  is displayed in <ref> (right). All show an exponential decrease. The weighting introduced in <ref> for _i does not influence the exponential decay rate of _i. The decrease in the classical singular values is directly linked to the exponential decrease of the projection error of POD full state due to properties of the Frobenius norm (see <cit.>). A similar result was deduced in the scope of the present paper for PSD SVD-like decomposition and the PSD functional (see <ref>).





  §.§.§ Orthonormality and symplecticity of the bases

To verify the orthonormality and the symplecticity numerically, we consider the two functions

    (2k) =  - 2k,      (2k) =  - 2k

which are zero / numerically zero if and only if the basis is orthonormal or symplectic, respectively. In <ref>, we show both values for the considered basis generation techniques and RB sizes.

The orthonormality of the bases is in accordance with the theory. All procedures compute symplectic bases except for PSD SVD like-decomposition. PSD greedy shows minor loss in the orthonormality which is a known issue for the -orthogonalization method used (modified symplectic Gram-Schmidt procedure with re-orthogonalization <cit.>). But no major impact on the reduction results could be attributed to this deficiency in the scope of this paper.

Also the symplecticity (or -orthogonality) of the bases behaves as expected. All PSD methods generate symplectic bases whereas the POD methods do not. A minor loss of symplecticity is recorded for PSD SVD-like decomposition which is objected to the computational method that is used to compute an SVD-like decomposition. Further research on algorithms for the computation of an SVD-like decomposition should improve this result. Nevertheless, no major impact on the reduction results could be attributed to this deficiency in the scope of this paper.





  §.§.§ Relative error in the reduced model

We investigate the error introduced by MOR in the reduced model. The error is measured in the relative ∞-norm ∙ in time and space

    (2k, ) :=
    max_i∈ 1,…,(t_i, ) - (t_i, )/max_i∈ 1,…,(t_i, ),

where 2k indicates the size of the ROB ∈^2n × 2k and ∈ is one of the test parameters. To display the results for all 16 test parameters at once, we use box plots in <ref>. The box represents the 25%-quartile, the median and the 75%-quartile. The whiskers indicate the range of data points which lay within 1.5 times the interquartile range (IQR). The crosses show outliers. For the sake of a better overview, we truncated relative errors above 10^0 = 100%.

The experiments show that the non-symplectic MOR techniques show a strongly non-monotonic behaviour for an increasing basis size. For many of the basis sizes, there exists a parameter which shows crude approximation results which lay above 100% relative error. The POD full state is unable to produce results with a relative error below 2%.

On the other hand, the symplectic MOR techniques show an exponentially decreasing relative error. Furthermore, the IQRs are much lower than for the non-symplectic methods. We stress that the logarithmic scale of the y axis distorts the comparison of the IQRs – but only in favour of the non-symplectic methods. The low IQRs for the symplectic methods show that the symplectic MOR techniques derive a reliable reduced model that yields good results for any of the 16 randomly chosen test parameters. Furthermore, none of the systems shows an error above 0.19% – for PSD SVD-like decomposition this bound is 0.018%, i.e. one magnitude lower.

In the set of the considered symplectic, orthogonal MOR techniques, PSD greedy shows the best result for most of the considered ROB sizes. This superior behaviour of PSD greedy in comparison to PSD complex SVD is unexpected since PSD greedy showed inferior results for the projection error in <ref>. This was also observed in <cit.>.

Within the set of investigated symplectic MOR techniques, PSD SVD-like decomposition shows the best results followed by PSD greedy and PSD complex SVD. While the two orthonormal procedures show comparable results, PSD SVD like-decomposition shows an improvement in the relative error. Comparing the best result of either PSD greedy or PSD complex SVD with the worst result of PSD SVD-like decomposition considering the 16 different test parameters for a fixed basis size – which is pretty much in favour of the orthonormal basis generation techniques –, the improvement of PSD SVD-like decomposition ranges from factor 3.3 to 11.3 with a mean of 6.7.













§ SUMMARY AND CONCLUSIONS




We gave an overview of autonomous and non-autonomous Hamiltonian systems and the structure-preserving model order reduction (MOR) techniques for these kinds of systems <cit.>. Furthermore, we classified the techniques in orthonormal and non-orthonormal procedures based on the capability to compute a symplectic, (non-)orthonormal reduced order basis (ROB). To this end, we introduced a characterization of rectangular, symplectic matrices with orthonormal columns. Based thereon, an alternative formulation of the PSD Complex SVD <cit.> was derived which we used to prove the optimality with respect to the PSD functional in the set of orthonormal, symplectic ROBs. As a new method, we presented a symplectic, non-orthonormal basis generation procedure that is based on an SVD-like decomposition <cit.>. First theoretical results show that the quality of approximation can be linked to a quantity we referred to as weighted symplectic singular values.

The numerical examples show advantages for the considered linear elasticity model for symplectic MOR if a symplectic integrator is used. We were able to reduce the error introduced by the reduction with the newly introduced non-orthonormal method.

We conclude that non-orthonormal methods are able to derive bases with a lower error for both, the training and the test data. Yet, it is still unclear if the newly introduced method computes the global optimum of the PSD functional. Further work should investigate if a global optimum of the PSD functional can be computed with an SVD-like decomposition.









 




























This research was partly funded by the German Research Foundation (DFG) grant number HA5821/5-1 and within the GRK 2198/1.


We thank the German Research Foundation (DFG) for funding this work. The authors thank Dominik Wittwar for inspiring discussions.

























































plain














