PRESERVING LAGRANGIAN STRUCTURE IN
NONLINEAR MODEL REDUCTION WITH
APPLICATION TO STRUCTURAL DYNAMICS

arXiv:1401.8044v1 [] 31 Jan 2014

KEVIN CARLBERG∗ , RAY TUMINARO† , AND PAUL BOGGS‡
Abstract. This work proposes a model-reduction methodology that preserves Lagrangian structure (equivalently Hamiltonian structure) and achieves computational efficiency in the presence of high-order nonlinearities and arbitrary parameter
dependence. As such, the resulting reduced-order model retains key properties such as energy conservation and symplectic timeevolution maps. We focus on parameterized simple mechanical systems subjected to Rayleigh damping and external forces,
and consider an application to nonlinear structural dynamics. To preserve structure, the method first approximates the system’s ‘Lagrangian ingredients’—the Riemannian metric, the potential-energy function, the dissipation function, and the external
force—and subsequently derives reduced-order equations of motion by applying the (forced) Euler–Lagrange equation with these
quantities. From the algebraic perspective, key contributions include two efficient techniques for approximating parameterized
reduced matrices while preserving symmetry and positive definiteness: matrix gappy POD and reduced-basis sparsification
(RBS). Results for a parameterized truss-structure problem demonstrate the importance of preserving Lagrangian structure
and illustrate the proposed method’s merits: it reduces computation time while maintaining high accuracy and stability, in
contrast to existing nonlinear model-reduction techniques that do not preserve structure.
Key words. nonlinear model reduction, structure preservation, Lagrangian dynamics, Hamiltonian dynamics, structural
dynamics, positive definiteness, matrix symmetry

1. Introduction. Computational modeling and simulation for simple mechanical systems characterized
by a Lagrangian formalism has become indispensable across a variety of industries. Such simulations enable
the understanding of complex systems, reduced design costs, and improved reliability for a wide range of
applications. For example, computational structural dynamics tools have become widely used in applications
ranging from aerospace to biomedical-device design; molecular-dynamics simulations have gained popularity
in materials science and biology, in particular. However, the high computational cost incurred by simulating
large-scale simple mechanical systems can result in simulation times on the order of weeks, even when using
high-performance computers. As a result, these simulation tools are impractical for time-critical applications
that demand the accuracy provided by large-scale, high-fidelity models. In particular, applications such as
nondestructive evaluation for structural health monitoring, multiscale modeling, embedded control, design
optimization, and uncertainty quantification require highly accurate results to be obtained quickly.
In this work, we consider models that depend on a set of parameters, e.g., design variables, operating
conditions. In this context, model-reduction methods present a promising approach for addressing timecritical problems. During the offline stage, these methods perform computationally expensive ‘training’
tasks, which may include evaluating the high-fidelity model for several instances of the system parameters
and computing a low-dimensional subspace for the solution. Then, during the inexpensive online stage, these
methods quickly compute approximate solutions for arbitrary values of the system parameters. To accomplish
this, they reduce the dimension of the high-fidelity model by restricting solutions to lie in the low-dimensional
subspace that was computed offline; they also introduce other approximations when nonlinearities are present.
Thus, the reduced-order model used online is characterized by a low-dimensional dynamical system that
arises from a projection process on the high-fidelity-model equations. This offline/online strategy is effective
primarily in two scenarios: ‘many query’ problems (e.g., Bayesian inference), where the high offline cost is
amortized over many online evaluations, and real-time problems (e.g., control) characterized by stringent
constraints on online evaluation time.
Generating a reduced-order model that preserves the Lagrangian structure intrinsic to mechanical systems is not a trivial task. Such structure is critical to preserve, as it leads to fundamental properties such as
energy conservation (in the absence of non-conservative forces), conservation of quantities associated with
symmetries in the system, and symplectic time-evolution maps. In fact, the class of structure-preserving
time integrators (e.g., geometric integrators [15], variational integrators [19]) has been developed to ensure
∗ Harry

S. Truman Fellow, Quantitative Modeling & Analysis Department ktcarlb@sandia.gov
Analysis and Applications Department, rstumin@sandia.gov
‡ Quantitative Modeling & Analysis Department (retired), ptboggs@sandia.gov

† Numerical

1

that the discrete solution to the high-fidelity computational model associates with the time-evolution map
of a (modified) Lagrangian system.
Lall et al. [18] show that performing a Galerkin projection on the Euler–Lagrange equation—as opposed
to the first-order state-space form—leads to a reduced-order model that preserves Lagrangian structure.
However, the computational cost of assembling the associated low-dimensional equations of motion scales
with the dimension of the high-fidelity model. For this reason, this approach is efficient only when the lowdimensional operators can be assembled a priori ; this occurs only in very limited cases e.g., when operators
have a low-order polynomial dependence on the state and are affine in functions of the parameters [22].
Several methods have been developed in the context of nonlinear-ODE model reduction that can reduce
the computational cost of assembling the low-dimensional equations of motion. However, these methods destroy Lagrangian structure when applied to simple mechanical systems. For example, collocation approaches
[4, 23] perform a Galerkin projection on only a small subset of the full-order equations characterizing the
high-fidelity model. Although this method works well for some nonlinear ODEs, it destroys Lagrangian
structure. The discrete empirical interpolation method (DEIM) [9, 14, 11] and gappy proper orthogonal
decomposition (POD) reconstruction methods [13, 6, 7] compute a few entries of the vector-valued nonlinear
functions, and then approximate the uncomputed entries by interpolation or least-squares regression using
an empirically derived basis. Galerkin projection can then be performed with the approximated nonlinear
function. Again, this technique destroys Lagrangian structure.
The goal of this work is to devise a reduced-order model for nonlinear simple mechanical systems with
general parameter dependence that leads to computationally inexpensive online solutions and preserves Lagrangian structure. We focus particularly on parameterized structural-dynamics models under Rayleigh
damping and external forces. The methodology we propose constructs a reduced-order model by first approximating the ‘Lagrangian ingredients’ (i.e., quantities defining the problem’s Lagrangian structure) and
subsequently deriving the equations of motion by applying the Euler–Lagrange equation to these ingredients.
The method approximates the Lagrangian ingredients as follows:
I. Configuration space. The low-dimensional configuration space is derived using standard dimensionreduction techniques, e.g., proper orthogonal decomposition, modal decomposition.
II. Riemannian metric. The Riemannian metric is defined by a low-dimensional symmetric positive-definite
matrix. We propose two efficient methods for approximating this low-dimensional matrix that preserve
symmetry and positive definiteness.
III. Potential-energy function. The potential energy function is approximated by employing the original potential-energy function, but with the low-dimensional reduced-basis matrix replaced by a lowdimensional sparse matrix with only a few nonzero rows. This sparse matrix is computed online by
matching the gradient of the reduced potential to first order about the equilibrium configuration.
IV. Dissipation function. The damping matrix associated with the Rayleigh dissipation function is a linear
combination of the mass matrix (which defines the Riemannian metric) and the Hessian of the potential.
Thus, we form the approximated Rayleigh dissipation function in the same fashion, but employ the
approximated mass matrix from ingredient II and approximated potential from ingredient III.
V. External force. The external force is derived by applying the Lagrange–D’Alembert principle with
variations in the configuration space. We approximate this by applying gappy POD reconstruction to
the external force as expressed in the original coordinates. As a result, the external force appearing in
the reduced-order equations of motion can be derived by applying the Lagrange–D’Alembert principle
to this modified external force with variations restricted to the low-order configuration space.
We note that a structure-preserving method [5] has been recently proposed for nonlinear port-Hamiltonian
systems, which are generalizations of Hamiltonian systems. While this technique guarantees that properties
such as stability and passivity are preserved, it does not in fact preserve Lagrangian or classical Hamiltonian structure. In particular, the resulting reduced-order equations of motion cannot be derived from
approximated ingredients such as those enumerated above; as a consequence, the approach does not ensure
symplecticity or energy conservation for conservative systems, for example.
As hinted above, preserving structure for Lagrangian ingredient II is equivalent to efficiently approximating a low-dimensional reduced matrix while preserving symmetry and positive definiteness. This algebraic
task is relevant to a broad scope of applications, e.g., approximating extreme eigenvalues/eigenvectors of a
parameterized matrix, preserving Hessian positive definiteness in optimization algorithms. For this reason,
2

Section 2 presents approximation techniques for Lagrangian ingredient II in a stand-alone algebraic setting
that does not rely on the Lagrangian formalism. Similarly, Section 3 considers Lagrangian ingredient III in
a purely alegraic context that does not depend on Lagrangian dynamics.
The remainder of the paper is organized as follows. Section 4 introduces the Lagrangian-mechanics
formulation. Section 5 outlines existing model-reduction techniques and highlights the need for an efficient,
structure-preserving method. Section 6 presents the proposed method. Section 7 presents numerical experiments applied to a simple mechanical system from structural dynamics. Finally, Section 8 summarizes the
contributions and suggests further research.
The structure-preserving model-reduction methods proposed in this work also preserve Hamiltonian
structure when the Hamiltonian formulation of classical mechanics is used. Appendix A provides this connection.
2. Preserving matrix symmetry and positive definiteness. This section presents approximation
techniques for Lagrangian ingredient II in an algebraic setting. First, to establish notation, denote the system
parameters by µ ∈ D, where D represents the parameter domain. Let A (µ) denote an N × N parameterized
symmetric positive-definite (possibly dense) matrix. Finally, let V denote a dense, parameter-independent,
full-column-rank N × n matrix with n  N whose columns can be interpreted as a reduced basis spanning
an n-dimensional subspace of RN . We consider the following online problem:
(P1) At a cost independent of N , compute a symmetric positive-definite matrix Ã (µ? ) that is
appropriately close to the matrix VT A (µ? )V for any specified online point µ? ∈ D.
Note that due to the density of V, directly computing VT A (µ? )V requires computing all O(N ) entries of
the matrix A (µ? ). Recall that the offline/online strategy we adopt permits expensive offline operations that
facilitate the solution to online problem (P1). These operations may include collecting p ‘snapshots’ of the
matrix A µi , i = 1, . . . , p, where µi ∈ D denotes the ith instance of the training set.
We assume that computing a single entry of A (µ? ) for any specified online point µ? ∈ D is inexpensive,
i.e., the number of floating-point operations (flops) is independent of N . However, we make no other
assumptions regarding the parameters or the matrix. In particular, we do not assume affine parametric
dependence of the matrix, and we view µ 7→ A (µ) simply as a mechanism for generating instances of the
matrix A.
We now present two methods for solving online problem (P1). Method 1 approximates the reduced
matrix by projecting the full matrix onto a sparse basis, while Method 2 approximates the reduced matrix as
a linear combination of pre-computed reduced matrices. Later, Section 5.2 constrasts the proposed methods
with existing model reduction approaches such as DEIM, gappy POD, and collocation. These existing
methods apply a one-sided sampling operator to the matrix, i.e., they replace VT with some type of sparse
matrix. While this leads to an inexpensive approximation, it gives rise to a non-symmetric reduced matrix
approximation. This destroys the underlying problem structure and therefore fails to meet the requirements
of online problem (P1).
2.1. Reduced-basis sparsification (RBS). We first consider a strategy that ‘injects sparseness’ into
the matrix V. That is, we replace V by UA ∈ RN ×n , which has full column rank and only m rows (with
n ≤ m  N ) containing nonzero entries:
Ã (µ) = UA T A (µ)UA .

(2.1)

This sparse matrix may be expressed as UA ≡ PUA , where P ∈ {0, 1}N ×m is a ‘sampling matrix’ consisting
is a dense matrix with full column rank
of m selected columns of the N × N identity matrix, UA ∈ Rm×n
∗
n, and Rm×n
denotes
the
noncompact
Stiefel
manifold:
the
set
of
full-rank
m × n matrices. Clearly, Ã (µ)
∗
is symmetric positive definite if A (µ) is symmetric positive definite; thus, the approximation defined by
(2.1) preserves the requisite structure. Note that this approximation will also preserve structure in cases
where A (µ) is symmetric positive semidefinite or simply symmetric. Further, the (online) operation count
for computing Ã (µ? ) for online point µ? ∈ D is independent of N . Computing PT A (µ? )P is equivalent
2
to computing only
entries of A (µ? ) and entails O(m2 ) flops; subsequently computing
 T m ?(symmetric)

T
?
Ã (µ ) = UA P A (µ )P UA requires O(m2 n + mn2 ) flops.
3

Given a sampling matrix P, the matrix UA can be computed offline to minimize the average approximation error over the snapshots, i.e., according to the following optimization problem:
UA = arg

min

m×n
X∈R∗

p
X



XT PT A µi PX − VT A µi V

2
,
F

(2.2)

i=1

where the subscript F denotes the Frobenius norm. To handle the fact that R∗m×n is an open set, optimization
problem (2.2) can first be solved over Rm×n and the solution can be subsequently projected onto Rm×n
, which
∗
is analogous to the approach taken by Vandereycken [25, Algorithm 6]. We note that other objective functions
may be considered for specialized online analyses, e.g., the approximation error of extreme eigenvalues or
eigenvectors over the matrix snapshots. Note that problem (2.2) is a small-scale optimization problem, as
m, n  N . It can be solved at a cost independent of N (for each optimization iteration) during
the offline

stage after the matrix snapshots A µi , i = 1, . . . , p and their reduced counterparts VT A µi V, i = 1, . . . , p
have been computed.
Procedure 1 provides the offline and online steps required to implement the RBS approximation.
Procedure 1: Reduced-basis sparsification for symmetric matrices
Offline stage
1
2
3
4


Collect matrix snapshots A µi , i =1, . . . , p.
Form reduced the matrices VT A µi V, i = 1, . . . , p.
Choose the sample matrix P.
Determine UA as the solution to problem (2.2).
Online stage (given µ? )

1
2

Compute PT A (µ? )P.


Form Ã (µ? ) = UA T PT A (µ? )P UA .

Remark. This paper does not focus on methods for selecting the sampling matrix P, a task that is typically
carried out during the offline stage and uses the collected snapshots. All numerical experiments presented in
Section 7 use the GNAT greedy approach [6] for this purpose. This method has proven to be quite robust,
even when applied to the approximation techniques proposed in this paper. A more careful study of sampling
algorithms will be addressed in future work, where ideas of tailoring the sampling approach to the specific
reduced-order modeling approximation will be explored.
2.1.1. Exactness conditions. In the full-sampling case where m = N , the approximation is exact if
problem (2.2) is solved via a gradient-based method with an initial guess of X(0) = PT V; we do this in
practice. Under these conditions, UA = V and so Ã (µ) = VT A (µ)V.
In the general case where m < N , it is possible to show that the approximation is exact if the matrix
is parameter-independent (i.e., A (µ) = A) and m ≥ n. This situation is considered in the discussion that
follows Theorem 2.1 below. It is also possible to prove a more general exactness result in cases where the
parametric dependence of A (µ) is relatively simple. In particular, consider the parametric form
A (µ) = h1 (µ) A1 + h2 (µ) A2 ,

(2.3)

where A1 and A2 are N × N symmetric positive-definite matrices and h1 , h2 : D → R. It can then be
shown that a sparse reduced basis exists that exactly captures VT A (µ)V under conditions related to how
well the eigenvalues of the sampled matrix PT A(µ)P encompass (or surround) those of the reduced matrix
VT A(µ)V. Loosely stated, the encompassing conditions amount to how well the sampled matrix captures the
behavior of the reduced matrix. Formally, the following theorem makes precise the notion of encompassing
using eigenvalue interlacing ideas from classical linear algebra.
Theorem 2.1. Let A (µ) have the form given by Eq. (2.3).
4

Then,
UA T PT A (µ) P UA = VT A (µ) V,

∃ UA ∈ Rm×n such that

∀µ ∈ D

(2.4)

if and only if the generalized eigenvalues of (VT A2 V, VT A1 V) interlace the generalized eigenvalues of
(PT A2 P, PT A1 P), i.e.,
(s)

λi

(r)

≤ λi

(s)

≤ λi+m−n

f or i = 1, ..., n,

(2.5)

with
 T
 (r)
 (r)
(r) 
V A2 V xi = λi VT A1 V xi ,

i = 1, . . . n

and
 T
 (s)
 (s)
(s) 
P A2 P xi = λi PT A1 P xi ,

i = 1, . . . , m.

Note that the eigenvalues are sorted in order of increasing magnitude. Appendix D.1 contains the proof.
Here, we discuss the theorem’s implications.
When A is independent of µ, we can choose h1 = h2 = 1 and A1 = A2 . The interlacing property is then
(s)
(r)
trivially satisfied for m = n with λi = λi = 1, and so the equality in (2.4) always holds. When instead
(s)
(r)
(s)
A1 6= A2 and m=n+1, the interlacing definition is quite restrictive, as it implies that λk ≤ λk ≤ λk+1 .
We would not generally expect the eigenvalues of the sampled and reduced matrices to have this property.
However, when m  n + 1, each interval width is (much) larger and so the condition is not nearly as
(s)
(r)
(s)
restrictive. For example, if n = 100 and m = 300, then interlacing implies that λi ≤ λi ≤ λi+200 . Thus,
we generally expect the conditions of the theorem to be satisfied for sufficiently large m, though this is not
guaranteed and depends on matrix spectrum. As a final note, although the theorem assumes a specific form
of A (µ), it should characterize the rough behavior of a more general A (µ) that does not vary ‘too much’
from the affine functional form (2.3).
2.2. Matrix gappy POD. An alternative structure-preserving approximation applicable to problem
(P1) assumes the following form:
Ã (µ) =

nA
X

i
ξA
(µ)VT Ai V.

(2.6)

i=1

Here, the matrices Ai , i = 1, . . . , nA are N × N symmetric matrices that are computed offline and define a
basis for the matrix A (µ). Due to the symmetry of Ai , i = 1, . . . , nA , the
 approximation Ã (µ) will always
nA
1
be symmetric. The parameter-dependent coefficients ξA ≡ ξA
, . . . , ξA
are computed online in an efficient
manner that ensures Ã (µ) is positive definite and thereby preserves requisite structure.
The next sections describe procedures for computing the matrix basis and coefficients. We refer to this
method as ‘matrix gappy POD’, as it amounts to the gappy POD procedure [13] applied to matrix data
with modifications to preserve positive definiteness. The approach, which we originally proposed [8], is a
more general formulation of the ‘matrix DEIM’ approach [26] (or ‘multi-component EIM’ [24] in the context
of the reduced-basis method applied to parametrized non-affine elliptic PDEs), as it permits least-squares
reconstruction (not simply interpolation). Further, it is equipped with a mechanism to maintain positive
definiteness.
2.2.1. Offline computation: matrix basis. To obtain the matrix basis, we propose applying a
vectorized POD method, wherein the basis can be considered a set of ‘principal matrices’ that optimally
i
represent1 the matrix A over the training
 set {µ }. The (offline) steps for this method are as follows:
i
1. Collect matrix snapshots A µ , i = 1, . . . , p.
1 These matrices are optimal in the sense that they minimize the average projection error (as measured in the Frobenius
norm) of the matrix snapshots.

5


2
2
2. Vectorize the snapshots ai ≡ v A µi ∈ RN , i = 1, . . . , p, where the function v : RN ×N → RN
vectorizes a matrix.
3. Compute an nA -dimensional (with nA ≤ p) POD basis of the vectorized snapshots


2
Wa ≡ a1 · · · anA ∈ RN ×nA

(2.7)

using vectorized snapshots {ai }pi=1 and an ‘energy criterion’ ηA ∈ [0, 1] as inputs to Algorithm 1 of
Appendix B.
4. Transform these POD basis vectors into their matrix counterparts:

Ai = v −1 ai , i = 1, . . . , nA .
(2.8)
Each matrix Ai , i = 1, . . . , nA is guaranteed to be symmetric, as Algorithm 1 forms this basis by taking
linear combinations of symmetric matrices.
2.2.2. Online computation: coefficients. The approximation error can be bounded as follows:
kVT A(µ)V − Ã (µ) kF = kVT A(µ)V −

nA
X

i
ξA
(µ)VT Ai VkF

i=1
nA
X

≤ kVk2F A (µ) −

i
ξA
(µ)Ai

F

(2.9)
(2.10)

i=1

where kVk2F = n if V is orthogonal. This leads to a natural choice for the scalar coefficients based on
minimizing the upper bound (2.10). In particular, we compute coefficients ξA (µ? ) online for a specific
µ? ∈ D as the solution to
nA
X
xi PT Ai Pk2F
minimize kPT A (µ? )P −
(x1 ,...,xnA )
i=1
nA
X
xi VT Ai V > 0.
subject to

(2.11)

i=1

Note that the coefficients are computed to match (as closely as possible) the full matrix and the linear
combination of pre-computed full matrices at a few entries. The constraints amount to a strict linear-matrixinequality, where A > 0 denotes a generalized inequality that indicates A is a positive-definite matrix. This
constraint ensures that structure is preserved. Note that the constraint can be modified (resp. dropped) in
cases where positive semidefiniteness (resp. simply symmetry) aims to be preserved.
Problem (2.11) is equivalent to a linear least-squares problem with nonlinear constraints; this can be
seen from its vectorized form:
minimize

T

x=[x1 ··· xnA ]

subject to

kP̄T v (A (µ? )) − P̄T Wa xk22
nA
X

(2.12)
T

xi V Ai V > 0.

i=1

Here, P̄ is an alternate form of the sampling matrix that can be applied to vectorized matrices, i.e.
P̄T v (A (µ? )) = v PT A (µ? )P .2
The objective function is equivalent to that of the gappy POD method [13]—which will be further
discussed in Section 5.2.2—applied to matrix data. Note that this optimization problem is solved online
h
i
2
2
2
symmetry, this sampling matrix can be expressed as P̄ ≡ p̄1 · · · p̄(m +m)/2 ∈ {0, 1}N ×(m +m)/2 , where
   


2
T
= v pi pj
for i = 1, . . . , j and j = 1, . . . , m and P ≡ p1 · · · pm . From the the definition of p̄i+(j −j)/2 , it

2 Exploiting

p̄i+(j

2

−j)/2

follows that p̄i+(j

2

−j)/2

extracts the (i, j)th entry from the vectorized form of a matrix.
6

using the online-sampled data PT A (µ? )P; Appendix C describes a method for solving this optimization
problem. In practice, we usually observe the constraints to be inactive at the unconstrained solution.
Therefore, typically the constraints need not be handled directly, and solving problem (2.11) amounts to
solving a small-scale linear least-squares problem characterized by an (m2 + m)/2 × nA matrix. To ensure
a unique solution to problem (2.11), the matrix P̄T Wa must have full column rank. This can be achieved
by enforcing (m2 + m)/2 ≥ nA as well as mild conditions on the sampling matrix P.
Procedure 2 describes the offline and online stages for implementing the matrix gappy POD approximation.
Procedure 2: Matrix gappy POD
Offline stage
1

2

3
4

Compute the basis matrices Ai , i = 1, . . . , nA using the vectorized POD approach described in
Section 2.2.1.
Determine the sampling matrix P which gives rise to a full column rank matrix P̄T Wa and with m
chosen so that (m2 + m)/2 ≥ nA .
Compute low-dimensional matrices VT Ai V, i = 1, . . . , nA .
Retain the sampled entries of the matrix basis PT Ai P, i = 1, . . . , nA ; discard other entries.
Online stage (given µ? )

1
2
3

Compute PT A(µ? )P.
Solve the small-scale optimization problem (2.11) for coefficients ξA (µ? ).
Assemble the low-dimensional matrix Ã (µ? ) by Eq. (2.6).

2.2.3. Exactness conditions. Theorem 2.2. The matrix gappy POD approximation is exact for any
specified online parameters µ? ∈ D if
1. v (A (µ? )) ∈ range (Wa ) and
2. P̄T Wa has full column rank.
See Appendix D.2 for the proof.
Condition 1 holds, e.g., when µ? ∈ {µi } and nA = p. Condition 2 can be straightforwardly enforced by
the choice of P and automatically holds in the case of full sampling, i.e., m = N .
3. Preserving potential-energy structure. This section presents a technique for approximating
Lagrangian ingredient III within an algebraic setting. To begin, define a parameterized scalar-valued function
V : RN ×D → R with (q; µ) 7→ V that is nonlinear in both arguments and can be interpreted as a Lagrangian
dynamical system’s (parameterized) potential energy. Here, q ∈ RN denotes the system’s configuration
variables and µ ∈ D denotes the system parameters that belong to parameter domain D. Unlike the matrix
approximations of the previous section, the nonlinear dependence on the configuration variables q introduces
additional challenges that must be considered carefully.
We aim to devise an offline method—which may entail expensive operations—for constructing a scalarvalued function Ṽr : Rn × D → R. This function will be used online and should satisfy the demands of online
problem (P2):
(P2) Compute the gradient vector ∇qr Ṽr (q?r ; µ? ) at a cost independent of N . Given any online
parameters µ? ∈ D, this vector should be appropriately close to VT ∇q V (q̄ (µ? ) + Vq?r ; µ? )
for all coordinates q?r ∈ Rn .
As before, V represents a dense, parameter-independent, full-column-rank N × n matrix. We denote by
q̄ : D → RN a parameterized reference configuration about which the low-dimensional reduced configuration
space is centered. Notice that this problem is concerned with approximating the gradient of the scalar-valued
function as opposed to the function itself. As will be discussed in Section 5, this problem arises in model
reduction of parameterized Lagrangian-dynamics systems, where the gradient of the potential appears in the
equations of motion.
7

In certain specialized cases, the above approximation can be simplified considerably. For example, when
q̄(µ) = 0, ∀µ ∈ D and the function V (q; µ) is purely quadratic in its first argument, then VT ∇q V (q̄(µ) +
Vqr ; µ) = VT A (µ)Vqr , where A (µ) is a symmetric Hessian matrix; in this case, one of the approximation
techniques described in Section 2 can be straightforwardly
R applied. Alternatively, if the potential energy
is defined by the integral over a domain (i.e., V (q; µ) = Ω V(X, q; µ)dΩX ), a sparse cubature method [2]
can be used to achieve computational efficiency and structure preservation. In more general cases, however,
another approach is needed. In the following, we develop a method that makes no simplifying assumptions
about the dependence of the potential V on the configuration variables or parameters.
Due to the density of the matrix V, the most straightforward approach of setting Ṽr (qr ; µ) = V (q̄(µ) +
Vqr ; µ) leads to expensive online operations: computing the gradient ∇qr Ṽr (q?r ; µ? ) = VT ∇q V (q̄ (µ? ) +
Vq?r ; µ? ) requires first computing all N entries of the gradient vector ∇q V (q̄ (µ? ) + Vq?r ; µ? ). To rectify
this, we revisit the RBS technique proposed in Section 2.1 and introduce some minor modifications. In
×n
with only
particular, we replace V by a sparse parameter-dependent matrix UV (µ) ≡ PUV (µ) ∈ RN
∗
m×n
m  N nonzero rows, where UV (µ) ∈ R∗
is a dense matrix. That is, we approximate the potential
energy as
Ṽr (qr ; µ) ≡V (q̄(µ) + UV (µ)qr ; µ).

(3.1)

This approximation preserves structure, as Ṽr remains a parameterized scalar-valued function. Now, we wish
T
to compute UV such that ∇qr Ṽr (q?r ; µ? ) = UV (µ? ) ∇q V (q̄ (µ? ) + UV (µ? )q?r ; µ? ) is as close as possible
to VT ∇q V (q̄ (µ? ) + Vq?r ; µ? ) for any online point µ? ∈ D and any q?r ∈ Rn . One can imagine a variety of
methods for computing UV toward this stated goal. For example, one can formulate an optimization problem
to match the potential gradient at training points [8]; this effectively leads to a parameter-independent matrix
UV . However, we found this approach to lead to significant errors for many problems. Instead, we pursue
an idea motivated by the analysis in Section 3.1, which centers on the first two terms in a Taylor expansion
of VT ∇q V (q̄ (µ? ) + Vq?r ; µ? ) about the reference configuration.
In practice, we often find that the trajectories of dynamical systems are localized in the configuration
space. This is particularly true for mechanical oscillators often encountered in structural dynamics, where
the trajectory does not deviate drastically from the equilibrium configuration. Using this observation, we
focus our approximation efforts on accurately capturing the behavior of the potential in a neighborhood
of the online reference configuration q̄ (µ? ). Implicitly, this assumes that the online configurations do not
greatly diverge from this point. To this end, consider computing UV (µ? ) online such that the approximation
T
UV (µ? ) ∇q V (q̄ (µ? ) + UV (µ? )q?r ; µ? ) matches VT ∇q V (q̄ (µ? ) + Vq?r ; µ? ) to first order about the reference
configuration:
T

T

UV (µ? ) ∇q V (q̄ (µ? ) ; µ? ) + UV (µ? ) ∇qq V (q̄ (µ? ) ; µ? )UV (µ? )q?r
= VT ∇q V (q̄ (µ? ) ; µ? ) + VT ∇qq V (q̄ (µ? ) ; µ? )Vq?r ,

∀q?r ∈ Rn .

(3.2)

Notice that the high-order terms amount to approximating a reduced Hessian (defined via the dense matrix
V) by a second reduced Hessian (defined via the sparse matrix UV (µ? )). This is equivalent to online problem (P1) presented in Section 2 that was addressed by the RBS algorithm (as well as a matrix gappy POD
approach). This RBS algorithm is supported by Theorem 2.1, which shows that an exact approximation of
the reduced Hessian is possible under certain assumptions. While these assumptions do not always hold, the
theorem gives an expectation that a good approximation can be found under more general circumstances.
Unfortunately, the presence of the low-order terms in Eq. (3.2) alters the character of the reduced approximation and so Theorem 2.1 no longer applies. In this case, the matrix UV (µ? ) must serve to capture both
gradient and Hessian information simultaneously, which introduces restrictive assumptions in order to obtain
an equivalent result to Theorem 2.1; this will be shown in Lemma 1.
To avoid the limitations associated with these restrictions, we choose the reference configuration to be
equilibrium, i.e., q̄(µ) = q0 (µ) with equilibrium defined as ∇q V (q0 (µ) ; µ) = 0. This forces the low-order
Taylor terms to zero and simplifies Eq. (3.2) to
T

UV (µ? ) ∇qq V (q0 (µ? ) ; µ? )UV (µ? ) = VT ∇qq V (q0 (µ? ) ; µ? )V.
8

(3.3)

Now, Theorem 2.1 holds, implying that Equation (3.3) can be exactly solved when m = n. For this reason,
we compute UV (µ? ) online for each µ? to satisfy (3.3) using n sample indices. Specifically, we define it
according to


X
?
,
(3.4)
UV (µ ) =
0(m−n)×n
where X is given by solving
LT1 X = LT2 ,
L2 ∈ Rn×n denotes the lower-triangular Cholesky factor of VT ∇qq V (q0 (µ? ) ; µ? )V, L1 ∈ Rn×n denotes
the lower-triangular Cholesky factor of P1 T ∇qq V (q0 (µ? ) ; µ? )P1 , and P1 represents the first n columns
of P. We defer discussing the computational cost for this approach to Section 3.2, and now return to
the previously alluded difficulties associated with solving (3.2) when the reference configuration does not
correspond to equilibrium.
3.1. Solvability of the two-term Taylor equation. The method presented in the previous section
was motivated by difficulties in inexpensively approximating the reduced gradient of a nonlinear function.
In this section, we give some insight into these difficulties by investigating a much easier situation: the
solvability of the two-term Taylor equation (3.2), which we write in matrix/vector form as
UV T PT c + UV T PT APUV q?r = VT c + VT AVq?r ,

∀q?r ∈ Rn .

(3.5)

Here, we have set c = ∇q V (q̄ (µ? ) ; µ? ) and A = ∇qq V (q̄ (µ? ) ; µ? ). We have also dropped dependence on
µ? such that c and A are parameter independent in the following analysis; this is equivalent to restricting
equation (3.2) to a single instance of µ? . This is somewhat less than ideal in that we would normally wish to
minimize online costs by computing a single UV during the offline phase that is then valid for all subsequent
online calculations. However, what we now show is that it is not always possible to satisfy equation (3.2)
even when one is restricted to finding a UV for a single instance of µ? .
As (3.5) must hold for all q?r , we have the following two necessary and sufficient conditions
UV T PT APUV = VT AV

and UV T PT c = VT c.

(3.6)

m×n
g
It is possible to show that satisfying these conditions is equivalent to finding a U
such that
V ∈R
T

g
g
U
V UV = I

T

T
g
eT
and U
V P c̃ = V c̃

(3.7)

where
eTV
e = I.
V
g
e
The definitions of U
V , V, and c̃ are given below. The key point is that the necessary and sufficient conditions
T
g
g T
eT
for equation (3.7) amount to finding an orthogonal matrix, U
V , such that UV P c̃ = V c̃ for a given
e and a given vector, c̃. In the simple case when A is the identity and V is orthogonal,
orthogonal matrix V,
e = V. More generally, we have the following definitions:
g
we have UV = UV and V
−T
T T
g
U
V = P L PUV Lφ ,

e = LT VL−T ,
V
φ

c̃ = L−1 c,

where L is the lower-triangular Cholesky factor of A, and Lφ is the lower-triangular Cholesky factor of
VT AV. The above equivalence hinges on the identities PT LPPT LT P = PT AP and PT LPPT L−1 = PT .
These hold due to the lower-triangular form of the matrix L.
The following lemma addresses the conditions under which Eq. (3.7) or equivalently Eq. (3.5) hold.
Lemma 1. Consider the equations
T

g
g
U
V UV = I

and
9

T

T
g
e T c̃
U
c=V
V P e

(3.8)

with the following matrices given: P ∈ {0, 1}N ×m consists of selected columns of the identity matrix (see
m×n
e ∈ RN ×n with V
eTV
e = I, and c̃ ∈ RN ×1 . Then, assuming that m ≥ n, some U
g
prior definition), V
V ∈R
exists such that Eq. (3.8) is satisfied if and only if
e T c̃||2 = ||PT c̃||2 and m = n
||V

or

e T c̃||2 ≤ ||PT c̃||2 and m > n.
||V

(3.9)

See Appendix D.3 for the proof.
Obviously, equation (3.9) is satisfied if either c̃ = 0 (i.e., the equilibrium configuration is taken as the
reference configuration) or if m = N . Unfortunately, however, equation (3.9) is not guaranteed to be satise T c̃||2 ≤ ||PT c̃||2 corresponds
fiable in more general situations. Specifically, when m > n, the condition ||V
to comparing the magnitude of a vector of length n obtained by rotating and dropping components with a
e and c̃ are not correlated, then
second vector of length m obtained by simply dropping components. If V
one could perhaps hope that on average the vector with more components would generally have a larger
e and all
magnitude. However, when c̃ lies completely within the subspace spanned by the columns of V
T
e
components of c̃ are non-zero, then ||V c̃||2 = ||c̃||2 and so satisfying the necessary and sufficient conditions
requires ‘full sampling’ m = N . While this scenario may be considered pessimistic, one can expect that a
e In general, there is
very large value of m will be required when c̃ lies primarily in the range space of V.
no guarantee that even the simplified (i.e., parameter-independent) form of the two-term Taylor equation
is solvable. When one also considers that Eq. (3.5) corresponds to the restriction of Eq. (3.2) to a single
instance of µ? , the above result should be seen as quite discouraging.
For this reason, we abandon any attempt at computing a parameter-independent sparse matrix UV
during the offline phase that can serve to approximate the reduced gradient for all online points µ? ∈ D.
Instead, we limit ourselves to the online computation of a parameter-dependent matrix UV (µ? ) that is only
valid for a single point µ? but can be used for all reduced configuration variables q?r ∈ Rn that may arise
during the online evaluation, e.g., at each nonlinear iteration and time instance considered while numerically
solving the equations of motion. Additionally, we set the reference configuration to equilibrium, which results
in c̃ = 0 and guarantees solvability of the the two-term Taylor expression with m = n.
3.2. Implementation and cost. Procedure 3 summarizes the offline/online strategy for implementing
the RBS strategy for approximating the potential energy. This method satisfies the online computational
Procedure 3: Reduced-basis sparsification for potential energy
Offline stage
1

Determine the sampling matrix P.
Online stage (given µ? )

1
2
3
4

Compute VT ∇qq V (q0 (µ? ) ; µ? ) V.
Compute P1 T ∇qq V (q0 (µ? ) ; µ? ) P1 .
Solve Equation (3.4) for UV (µ? ).
For any q?r ∈ Rn , set Ṽr (q?r ; µ? ) = V (q0 (µ? ) + PUV (µ? ) q?r ; µ? ), and compute the gradient as
T
∇qr Ṽr (q?r ; µ? ) = UV (µ? ) PT ∇q V (q0 (µ? ) + PUV (µ? ) q?r ; µ? )

cost requirements of problem (P2) with one exception: online step 1 incurs an N -dependent operation count.
However, online steps 1–3 depend only on the online point µ? and not on the reduced configuration variables
q?r . Thus, these steps are performed only once per parameter instance, and their cost can be amortized over
all online-queried values of q?r . As a result, this does not preclude significant computational savings, as will
be shown in the numerical results reported in Section 7. Note that online step 2 is equivalent to computing
just O(n2 ) entries of ∇qq V , which can be completed at a cost independent of N . Step 3 requires O(n3 )
operations.
Remark. Most nonlinear reduced-order modeling methods [4, 23, 9, 14, 11, 6, 7] assume ‘H-independence’
[11], which states that the Jacobian of the vector-valued nonlinear function is sparse; in the present context,
this corresponds to sparsity of the matrix ∇qq V . When this assumption holds, the proposed methodology
10

incurs low online computational cost. This efficiency results from the fact that computing PT ∇q V (q0 (µ? ) +
PUV (µ? ) q?r ; µ? ) in Step 4 of Procedure 3 requires that only m components of the gradient ∇q V be evaluated;
if H-independence holds, then these m components depend on only O(m) components of the argument
q0 (µ? ) + PXqr , leading to an N -independent operation count.
Unfortunately, H-independence does not hold for some problems in Lagrangian dynamics. For example
molecular-dynamics models can be characterized by a potential that includes interaction terms between
all particles, resulting in a dense matrix ∇qq V . Here, the proposed method can still achieve efficiency by
‘centering’ the configuration space at equilibrium such that q0 (µ) = 0, ∀µ ∈ D. In this case, the method
requires computing only m components of the argument q0 (µ? ) + PXqr in Step 4 of Procedure 3 regardless
of the sparsity of the matrix ∇qq V . This efficiency is achievable due to the fact that the method injects
‘sparsification’ in the argument of the nonlinear function. This ability to achieve an N -independent operation
count when H-independence is violated distinguishes this method from others in the literature.
4. Lagrangian dynamics formulation. We have now developed techniques to approximate parameterized reduced symmetric-positive-definite matrices and potential functions. In this section, we show how
these methods enable us to achieve the objective of this work: preserving Lagrangian structure in model
reduction for nonlinear mechanical systems. We begin by presenting the Lagrangian-dynamics formulation
for such systems and highlighting critical problem structure. Later, Section 5 describes existing nonlinear
model-reduction techniques and explains how they destroy structure in this context. Section 6 presents
the proposed structure-preserving methodology, which employs the approximation techniques proposed in
Sections 2 and 3.
We consider parameterized, nonlinear simple mechanical systems, with a particular focus on structuraldynamics models constructed by a finite-element formulation. Such models are defined by a triple (Q, g, V )
parameterized by system parameters µ ∈ D. The parameters may describe variations in shape and material
properties, for example. The triple is composed of:
• A differentiable configuration manifold Q. We take Q = RN where N denotes the number of degrees
of freedom in the model, considered to be ‘large’ in this work.
• A parameterized Riemannian metric g(v, w; µ), where v and w belong to the tangent bundle of
Q. We take g(v, w; µ) = vT M (µ)w, where M (µ) denotes the N × N parameterized symmetric
positive-definite mass matrix.
• A parameterized potential-energy function V : Q × D → R.
The kinetic energy of a simple mechanical system can be expressed as T (q̇; µ) = 21 g(q̇, q̇; µ) =
1 T
2 q̇ M (µ)q̇, where q : [0, T] → Q denotes the time-dependent configuration variables and T denotes the
final time. This leads to the following expression for the Lagrangian, which represents the difference between
the kinetic and potential energies:
1
g (q̇, q̇; µ) − V (q; µ)
2
1
= q̇T M (µ)q̇ − V (q; µ).
2

L(q, q̇; µ) =

(4.1)
(4.2)

In many cases, the non-conservative forces3 consist of an applied external force and a dissipative force
arising from Rayleigh viscous damping. This dissipative force derives from a positive-semidefinite dissipation
function4
F (q̇; µ) ≡

1 T
q̇ C (µ)q̇,
2

(4.3)

where C (µ) = αM (µ) + β∇qq V (q0 (µ) ; µ) denotes a parameterized N × N symmetric positive-semidefinite
matrix with α ∈ R and β ∈ R. Here, q0 : D → RN denotes the (parameterized) equilibrium configuration
such that ∇q V (q0 (µ) ; µ) = 0. So, we consider non-conservative forces of the form f (q, q̇, t; µ) − ∇q̇ F (q̇; µ),
where f denotes the external force that is derived from the Lagrange–D’Alembert variational principle.
3 Conservative

forces can be handled by directly including them in the Lagrangian.
damped systems can also often be derived by a positive-semidefinite dissipation function [1].

4 Non-viscously

11

Given the Lagrangian (4.2), one can derive the equations of motion for a simple mechanical system
subject to an external force and Rayleigh viscous damping from the forced Euler–Lagrange equation
d
∇q̇ L(q, q̇; µ) − ∇q L(q, q̇; µ) = f (q, q̇, t; µ) − ∇q̇ F (q̇; µ).
dt

(4.4)

Substituting Eqs. (4.2) and (4.3) into Eq. (4.4) leads to the familiar equations of motion
M (µ)q̈ + C (µ)q̇ + ∇q V (q; µ) = f (q, q̇, t; µ) .

(4.5)

Conservative mechanical systems, where f (q, q̇, t; µ) = 0 and C (µ) = 0, exhibit important properties and
can be characterized using the Hamiltonian formulation of classical mechanics discussed in Appendix A. For
example, these systems conserve energy and quantities associated with symmetry, and their time-evolution
maps are symplectic. Because these properties are intrinsic characteristics of the mechanical systems, it
is desirable for numerical methods to preserve these properties. As mentioned in the introduction, the
class of structure-preserving time integrators has been developed for this purpose. This class of integrators
ensures that the numerical solution preserves essential properties such as energy conservation, momentum
conservation, and symplecticity [15, 19].
For this reason, we aim to develop a reduced-order model that preserves the structure of the mechanical
system, yet is computationally inexpensive to simulate. This will ensure that the reduced-order model
preserves these characteristic properties. Further, the reduced-order equations of motion for these can be
solved with a structure-preserving time integrator; this will ensure that the numerical solution computed
using the reduced-order model will also preserve these properties. The properties of the system we seek to
preserve are those enumerated in Section 1: a configuration space, a parameterized Riemannian metric, a
parameterized potential-energy function, a parameterized positive-semidefinite dissipation function, and an
external force derived from the Lagrange–D’Alembert principle. The first three properties constitute the
parameterized triple that ensures the model describes a simple mechanical system; the last two characterize
the non-conservative forces.
5. Existing model-reduction techniques. Model-reduction techniques aim to generate a lowdimensional model that is inexpensive to evaluate, yet captures the essential features of the high-fidelity
model. These methods first conduct a computationally expensive offline stage during which they perform
analyses (e.g., solving the equations of motion, modal analyses) for a training set {µi }pi=1 ⊂ D. Then, these
methods employ the data generated during these analyses to define a configuration manifold of reduced
dimension, as well as other approximations to achieve efficiency in the presence of nonlinearities or arbitrary
parameter dependence. This low-dimensional configuration manifold is subsequently employed to generate
a low-dimensional model that can be used to perform inexpensive analyses for any specified point µ? ∈ D
during the online stage.
When the configuration space is Euclidean (as is the case for the models considered herein), the configuration space of reduced dimension n  N can be expressed as
Qr ≡ {q̄(µ) + Vqr | qr ∈ Qr },

(5.1)

where q̄(µ) : D → RN denotes the (parameterized) reference configuration about which the affine reduced
subspace is centered, Qr = Rn , and V ∈ R∗N ×n defines the reduced basis represented as a (typically dense)
matrix. This leads to the following expression for the generalized coordinates and their derivatives:
q = q̄(µ) + Vqr
q̇ = Vq˙r

(5.2)

q̈ = Vq¨r .

(5.4)

(5.3)

Thus, the low-dimensional configuration space can be described in terms of low-dimensional generalized
coordinates qr ∈ Qr or in terms of original coordinates by Eq. (5.2). The basis V can be determined by a
variety of techniques, including proper orthogonal decomposition and modal decomposition.
12

5.1. Galerkin projection. Model reduction based on Galerkin projection preserves Lagrangian structure. As pointed out by Lall et al. [18], the Galerkin projection must be carried out on the Euler–Lagrange
equation (4.4)—not the first-order state-space form—in order to preserve this structure.
Following their approach, Galerkin-projection-based methods replace the original configuration space Q
by the reduced-order configuration space Qr and subsequently derive the equations of motion in the usual
way using a set of lower-dimensional generalized coordinates. In this way, the resulting model has an identical
structure to the original problem.
For simple mechanical systems subject to non-conservative forces, this amounts to defining the Lagrangian as
Lr (qr , q˙r ; µ) ≡ L(q̄(µ) + Vqr , Vq˙r ; µ)
1
= q˙r T VT M (µ)Vq˙r − V (q̄(µ) + Vqr ; µ)
2

(5.5)
(5.6)

and the dissipation function as
Fr (q˙r ; µ) ≡ F (Vq˙r ; µ)
1
= q˙r T VT C (µ)Vq˙r .
2

(5.7)
(5.8)

The external force, which is derived based on the Lagrange–D’Alembert variational principle, is transformed
by relation (5.2) into
fr (qr , q˙r , t; µ) ≡ VT f (q̄(µ) + Vqr , Vq˙r , t; µ) .

(5.9)

Following Section 4, the forced Euler–Lagrange equation applied to the Lagrangian Lr , the dissipation
function Fr , and the external force fr leads to the reduced-order equations of motion
d
∇q˙ Lr (qr , q˙r ; µ) − ∇qr Lr (qr , q˙r ; µ) + ∇q˙r Fr (q˙r ; µ) = fr (qr , q˙r , t; µ) .
dt r

(5.10)

This can be rewritten as
VT M (µ)Vq¨r + VT C (µ)Vq˙r + VT ∇q V (q̄(µ) + Vqr ; µ) = VT f (q̄(µ) + Vqr , Vq˙r , t; µ) .

(5.11)

Note that Eq. (5.11) could have also been derived by applying Galerkin projection to the original Euler–
Lagrange equation (4.5), i.e., making substitutions (5.2)–(5.4) and left multiplying the system of equations
by VT .
Thus, the Galerkin reduced-order model preserves the problem structure because it preserves all five
Lagrangian properties:
I. a configuration space Qr = Rn , which relates to the original configuration space by Eq. (5.1),
II. a parameterized Riemannian metric gr (vr , wr ; µ) = vrT VT M (µ)Vwr ,
III. a parameterized potential-energy function Vr (qr ; µ) = V (q̄(µ) + Vqr ; µ),
IV. a parameterized positive-semidefinite dissipation function Fr (q˙r ; µ) = 12 q˙r T VT C (µ)Vq˙r , and
V. an external force fr derived from applying the Lagrange–D’Alembert principle to the original external
force f , but restricted to variations in the configuration space Qr .
5.1.1. Computational bottleneck. Although the equations of motion (5.11) are low dimensional,
they remain computationally expensive to solve when the operators exhibit arbitrary parameter dependence
and the potential is nonlinear. The reason is simple: computing the low-dimensional components of (5.11)
incurs large-scale operations due to the density of V. For example, the following steps are required to
compute VT M (µ? )V for a specific µ? ∈ D during the online stage:
i. Compute M (µ? ), which incurs O(N ω) flops, where ω denotes the average number of nonzeros per row
of the matrix M (µ? ).
ii. Compute the product M (µ? )V, which incurs O(N ωn) flops.
iii. Compute the product VT (M (µ? )V), which incurs O(N n2 ) flops.
13

Thus, the cost scales with the large dimension N of the original configuration manifold. The same analysis
holds for the product VT C (µ? )V.
If the potential energy V exhibits a (general) nonlinear dependence on coordinates q, the situation
worsens. In this case, the vector ∇q V (q̄(µ) + Vqr ; µ? ) and product VT ∇q V must be computed for every
instance of qr . Similarly, VT f (q, q̇, t; µ? ) must be computed for every time instance. Thus, a dimension
reduction is generally insufficient to generate models with computational complexity
independent of N .
P
Remark. If the mass matrix is affine in functions of the parameters M (µ) = i αi (µ)M
P i with αi : D → R
and Mi ∈ RN ×N , then products VT Mi V can be assembled offline, and VT M (µ? )V = i αi (µ? ) VT Mi V
can be computed in O(n2 ) floating-point operations (flops) during the online stage [17, 21]. Similar lowcomplexity results can be obtained for the other terms if they can be similarly expressed in separable form.
However, affine parameter dependence is a quite limiting scenario and does not generally hold.
5.2. Complexity reduction. Several techniques have been developed to mitigate the computational
bottleneck described in Section 5.1.1. Before applying projection, these methods compute (or sample) only
a few entries of the vector-valued functions; other entries are not computed. In effect, this complexityreduction strategy is equivalent to employing a sparse left-projection test basis. Such methods have been
successfully applied to ODEs that do not exhibit particular structure. However, when applied to mechanical
systems described by Lagrangian mechanics, these techniques destroy Lagrangian structure.
5.2.1. Collocation. Collocation approaches [4, 23] compute only a subset of the full-order equations
of motion (4.5) before applying Galerkin projection. That is, the reduced-order equations of motion (5.11)
are approximated by
VT PPT M (µ)Vq¨r + VT PPT C (µ)Vq˙r + VT PPT ∇q V (q̄(µ) + Vqr ; µ)
= VT PPT f (q̄(µ) + Vqr , Vq˙r , t; µ).

(5.12)

Recall that the sampling matrix P consists of m selected columns of the identity matrix. If one considers the
matrix PPT V as defining a basis for a test space, Eq. (5.12) can be viewed as a Petrov–Galerkin projection.
Computing the components of Eq. (5.12) is inexpensive in the case of H-independence, i.e., when the
matrices
 M, C, ∇qqV , ∇q f , and ∇q̇ f are sparse. To see this, consider the first term in Eq. (5.12): computing
VT P PT M (µ? ) V for specific online point µ? ∈ D incurs O(mωn + mn2 ) flops when operations are
carried out in the order implied by the parentheses. This cost is small if the sparsity measure of M is small,
i.e., ω  N .
However, this cost-reduction approach destroys the problem’s structure, as it does not preserve the
following Lagrangian properties described in Section 4:
II. The approximated reduced mass matrix VT PPT M (µ)V is not symmetric, so it does not define a
metric.
III. The term VT PPT ∇qq V (q̄(µ) + Vqr ; µ)V is not symmetric, so it cannot be the Hessian of a
potential-energy function.
IV. The approximated reduced damping matrix VT PPT C (µ)V is not symmetric, so it does not derive
from a dissipation function.
Note that Property I is trivially satisfied, as the configuration space can be described as Qr = Rn and
relates to the original configuration space by Eq. (5.1). Further, Property V is satisfied, because the nonconservative forces can be derived by applying the Lagrange–D’Alembert variational principle to a modified
external force PPT f (q, q̇, t; µ), but restricted to variations in the (true) configuration space Qr .
5.2.2. DEIM/gappy POD. Methods based on the discrete empirical interpolation method [9, 14, 11]
or gappy POD [13, 6, 7] approximate via least-squares regression or interpolation the nonlinear vector-valued
functions appearing in Eq. (4.5); these include M (µ)q̈, C (µ)q̇, ∇q V (q; µ), and f (q, q̇, t; µ). Because these
approaches construct a separate approximation for each term in the governing equations, they often achieve
higher accuracy than collocation.
During the offline stage, these methods construct an orthogonal basis Wθ ∈ RN ×nθ with nθ ≤ m for each
nonlinear function θ(t; µ) appearing in the equations of motion. The basis Wθ can be computed empirically
via proper orthogonal decomposition (POD), in which case the approximation technique is referred to as
‘gappy POD’ [13]. This consists of two steps: 1) collect snapshots Xθ = {θ(t; µ) | t ∈ Tsample (µ), µ ∈ {µi }},
14

where Tsample (µ) ⊂ [0, T] designates the time instances taken by the time-integration method for the training
simulation; and 2) compute Wθ by Algorithm 1 of Appendix B using Xθ and an energy criterion ηθ ∈ [0, 1]
as inputs.
During the online stage, these methods approximate the nonlinear function as
θ(t; µ) ≈ Wθ [PT Wθ ]+ PT θ(t; µ)

(5.13)

where a superscript + denotes the Moore–Penrose pseudoinverse and [PT Wθ ]+ PT θ is simply the solution
to the linear least-squares problem
minimize
kPT θ − PT Wθ θ r k22 .
n
θ r ∈R

θ

(5.14)

Notice that when nθ = m, the least-squares residual is zero (assuming the PT Wθ has full column rank) and
so the above procedure corresponds to interpolation.
As with collocation, this approximation technique leads to computational-cost savings during the online
stage if computing PT θ(t; µ) incurs a flop count independent of N , i.e., θ(t; µ) exhibits H-independence.
Substituting least-squares approximations for the nonlinear functions into Eq. (5.11) yields the approximated
reduced-order equations of motion
YMq¨r M (µ)Vq¨r + YCq˙r C (µ)Vq˙r + Y∇q V ∇q V (q̄(µ) + Vqr ; µ) = Yf f (q0 (µ) + Vqr , Vq˙r , t; µ) . (5.15)
Here, we have used the notation

+
Yθ ≡ VT Wθ PT Wθ PT ,

(5.16)

and the subscript of Y and W denotes the function for which the approximation has been constructed.
Unfortunately, this approximation method also destroys the Lagrangian structure. As before, Lagrangian
properties II–IV are lost because the reduced mass, stiffness, and damping matrices are not symmetric.
However, Property I is preserved. Property V is also preserved, because the non-conservative external force
can be derived by the Lagrange–D’Alembert principle applied to the modified external force f̃ (q, q̇, t; µ) =

+
Wθ PT Wf PT f (q, q̇, t; µ) with variations restricted to the configuration space Qr .
6. Efficient, structure-preserving model reduction. The main idea of the proposed approach is
to directly approximate the quantities defining the Lagrangian structure of the Galerkin-projection reducedorder model, and subsequently derive the equations of motion. Section 5.1 enumerates these quantities for the
simple mechanical systems considered herein: the Riemannian metric gr , the potential-energy function Vr , the
semidefinite dissipation function Fr , and the external force fr . Approximations to these ingredients should
1) preserve salient properties, 2) lead to computationally inexpensive reduced-order-model simulations, and
3) incur minimal approximation error.
To this end, we propose a model defined by
I. a configuration space Qr = Rn , which relates to the original coordinates by Eq. (5.1),
II. an approximated Riemannian metric g̃r ,
III. an approximated potential-energy function Ṽr ,
IV. an approximated positive-semidefinite dissipation function F̃r , and
V. an approximated external force f˜r derived from applying the Lagrange–D’Alembert principle to an
approximated force f̃ represented in the original coordinates, but limited to variations in the reduced
configuration space Qr .
We can derive the equations of motion by applying the forced Euler–Lagrange equation with these approximations:
d
∇q˙ L̃r (qr , q˙r ; µ) − ∇qr L̃r (qr , q˙r ; µ) + ∇q˙r F̃r (q˙r ; µ) = f˜r (qr , q˙r , t; µ) ,
dt r

(6.1)

where the approximated Lagrangian is defined as
L̃r (qr , q˙r ; µ) ≡

1
g̃r (q˙r , q˙r ; µ) − Ṽr (qr ; µ).
2
15

(6.2)

reduced
Lagrangian
ingredients

reduced
Lagrangian
ingredients

apply
Euler–Lagrange
equation

reduced-order
equations of
motion

approximated
reduced
Lagrangian
ingredients
apply
Euler–Lagrange
equation

approximated
reduced-order
equations of
motion

approximated
reduced-order
equations of
motion

(a) Existing complexity-reduction methods (see Section 5.2). By approximating the equations of motion, such
methods destroy Lagrangian structure.

(b) Proposed approach. By approximating Lagrangian
ingredients before deriving the equations of motion, the
approach preserves Lagrangian structure.

Fig. 6.1. Comparing existing complexity-reduction approaches with the proposed approach. A dashed arrow implies a
complexity-reduction approximation.

Note that Eq. (6.1) approximates Eq. (5.10), while Eq. (6.2) approximates Eq. (5.6).
Figure 6.1 depicts the strategy graphically. The next sections describe two proposed methods that align
with this strategy for structure preservation. For reference, Table 6.1 reports components of the equations
of motion for these methods, as well as for the model-reduction methods discussed in the previous sections.

Galerkin

mass
matrix
VT M (µ)V

damping
matrix
VT C (µ)V

potential-energy
gradient
VT ∇q V (q0 (µ) + Vqr ; µ)

external
force
VT f

struct.
pres.?
yes

low
cost?
no

collocation

VT PPT M (µ)V

VT PPT C (µ)V

VT PPT ∇q V (q0 (µ) + Vqr ; µ)

VT PPT f

no

yes

gappy POD

YMq¨r M (µ)V

YCq˙r C (µ)V

Y∇q V ∇q V (q0 (µ) + Vqr ; µ)

Yf f

no

yes

proposal 1

UM T M (µ)UM

αUM T M (µ)UM +
βUV T ∇qq V (q0 (µ) ; µ)UV

UV T ∇q V (q0 (µ) + UV qr ; µ)

Yf f

yes

yes

UV T ∇q V (q0 (µ) + UV qr ; µ)

Yf f

yes

yes

method

proposal 2

n
M
P
i=1

i
ξM
(µ)VT Mi V

n
M
P

i
ξM
(µ)VT Mi V+
i=1
T
βUV ∇qq V (q0 (µ) ; µ)UV

α

Table 6.1
Terms appearing in the equations of motion for various model-reduction techniques, including the two proposed structurepreserving methods.

6.1. Riemannian-metric approximation g̃r . The function gr : (vr , wr ; µ) 7→ vrT VT M (µ)Vwr is
defined by a low-dimensional symmetric positive-definite matrix VT M (µ)V. Thus, the task of approximating this matrix is consistent with problem (P1) of Section 2; we therefore propose computing an approximated
Riemannian metric g̃r : Rn × Rn × D → R as
g̃r (v1 , v2 ; µ) ≡v1T M̃ (µ) v2 ,

(6.3)

where M̃ (µ) is an n × n matrix that must be symmetric and positive definite. The first method (proposal
1 in Table 6.1) employs the reduced-basis sparsification technique, i.e., it approximates M̃ (µ) via Eq. (2.1).
The second method (proposal 2 in Table 6.1) employs matrix gappy POD and approximates this matrix by
Eq. (2.6). Procedures 1 (Section 2.1) and 2 (Section 2.2) provide the offline and online steps to implement
these approximations.
16

6.2. Potential-energy-function approximation Ṽr . Noting that only ∇qr Ṽr appears in the reducedorder equations of motion (see Eqs. (6.1)–(6.2)), problem (P2) of Section 3 applies to this scenario, and so we
approximate the potential energy according to the method described in that section. In particular, Eq. (3.1)
defines the approximated reduced potential energy. Further, we set the reference configuration to equilibrium
q̄ = q0 to avoid the limitations associated with other choices (see the discussion in Section 3). Procedure 3
of Section 3.2 describes the offline/online decomposition for implementing this approximation.
6.3. Dissipation-function approximation Fr . To maintain the Rayleigh-damping structure, we
simply approximate the damping matrix as a linear combination of the approximated mass matrix and
Hessian of the potential at equilibrium
β
α T
v M̃ (µ) v + vT ∇qr qr Ṽr (0; µ)v,
(6.4)
2
2
where α and β are the Rayleigh damping coefficients defined in Section 4.
6.4. External-force approximation f˜r . The following form of the approximated external force preserves structure, i.e., ensures it is derived from applying the Lagrange–D’Alembert principle to an approximated force f̃ limited to variations in the reduced configuration space Qr :
F̃r (v; µ) =

f˜r (qr , q˙r , t; µ) ≡ VT f̃ (q, q̇, t; µ) .

(6.5)

Thus, the task of generating this approximation can be reduced to computing f̃ (q, q̇, t; µ)—an approximation
to the vector-valued function f (q, q̇, t; µ). That is, we assign no special mathematical properties to f aside
from the fact that it is a vector. One way to accomplish this is by the DEIM/gappy POD approach described
in Section 5.2.2.
The error in this approximation can be bounded using a result derived from the error in the gappy POD
approximation of f (e.g., see [7, Appendix D]). We obtain



+

kf˜r − fr k2 ≤ k I − Wf PT Wf PT f k2 ≤ kR−1 k2 k I − Wf WfT f k2 ,
(6.6)
where Wf is an orthogonal basis used to represent the external force, and PT Wf = QR is the thin QR matrix
factorization. This result assumes that PT Wf has full rank. Thus, the accuracy of this approximation relies
both on the sampling matrix PT and how close f is to the range of Wf . To achieve accuracy, we compute
Wf via POD, which minimizes the average value of k I − Wf WfT f k22 over the training data.
6.4.1. Exactness conditions. Exactness conditions are similar to those described in Section 2.2.3 for
the matrix gappy POD approximation. In the general case where m < N , if f = 0, then the approximation
is exact, i.e., f˜r = fr . If instead f has at least one non-zero entry, then sufficient conditions for an exact
approximation are 1) f ∈ range (Wf ) and 2) PT Wf has full column rank. The first of these conditions
holds, for example, when Wf is computed via POD, the POD basis is not truncated, f is independent of q
and q̇, µ? ∈ {µi }, and if a snapshot of the external force was collected at the considered time instance. The
second of these can be enforced by the method for choosing P, which is beyond the scope of this paper. In
the full-sampling case where m = N , condition 2 holds automatically, so we only require condition 1 in this
case.
6.4.2. Implementation. Procedure 4 provides the offline and online steps for implementing the
external-force approximation.
7. Numerical experiments. Although the Galerkin and proposed reduced-order models have a theoretical advantage over the gappy POD and collocation reduced-order models in terms of preserving Lagrangian structure, it is unclear if this translates to improved numerical results in practice. This section
investigates this question by applying the model-reduction techniques to a practical problem: the clamped–
free truss structure shown in Figure 7.1.
3
We set the material properties to those of aluminum, i.e., density ρ = 2700 kg/m and elastic modulus
9
E = 62 × 10 Pa. The external force is composed of four components:
f (µ, t) =

4
X
i=1

17

ri (µ, t)ri ,

(7.1)

Procedure 4: External-force approximation via gappy POD
Offline stage
1
2
3
4

Collect snapshots of the external force Xf ≡ {f (q, q̇, t; µ) | µ ∈ {µi }, t ∈ Tsample (µ)}
Compute a POD basis Wf using Algorithm 1 with inputs Xf and ηf ∈ [0, 1].
Determine the sampling matrix P.

+
Compute the low-dimensional matrix Yf = VT Wf PT Wf .
Online stage (given µ? )

1
2

?
Compute m  N entries of the external force PT f (q, q̇, t;
 µT ).

Compute the low-dimensional matrix–vector product Yf P f (q, q̇, t; µ? ) .

width
length
force 1
force 2
force 3
force 4

height

y

1m
1m

1m

z

x

Fig. 7.1. Clamped–free parameterized truss structure

where ri ∈ RN , i = 1, . . . , 4 correspond to unit loads uniformly distributed across designated nodes and
ri : D × [0, T] → R, i = 1, . . . , 4 denote the component-force magnitudes. Figure 7.1 depicts the spatial
distribution of the forces, which lead to vectors ri , i = 1, . . . , 4 through the finite-element formulation
described below. The parameterized, time-dependent magnitudes of these forces are
(
γi (µ) sin (λi (µ) (t − T/4)) , t ≥ T/4
ri (µ, t) =
,
(7.2)
0, otherwise
where γi : D → R and λi : D → R, i = 1, . . . , 4 denote the maximum force magnitudes and force frequencies,
respectively. Similarly, the initial condition is composed of four components
q(0; µ) =

4
X

si (µ)si ,

(7.3)

i=1

where si is the steady-state displacement of the truss subjected to load ri γi (µ̄) with µ̄ = (0, . . . , 0) denoting the nominal point in parameter space. The equilibrium configuration is simply the undeformed truss
represented by q0 (µ) = 0; thus, the configuration space is centered at equilibrium.
16
The truss is parameterized by 16 parameters µ ≡ (µ1 , . . . , µ16 ) ∈ [−1, 1] that affect the geometry,
initial condition, and applied force as described in Table 7.1.
The problem is discretized by the finite-element method. The model consists of sixteen three-dimensional
bar elements per bay with three degrees of freedom per node; this results in 12 degrees of freedom per bay.
We consider a problem with 250 bays, which leads to N = 3 × 103 degrees of freedom in the full-order
model. The bar elements model geometric nonlinearity, which results in a high-order nonlinearity in the
potential-energy V .
18

length (m)
200 + 50µ1

bar
cross-sectional
area (m2 )
0.0025(1 + 0.5µ2 )

width (m)

height (m)

10(1 + µ3 )

10(1 + µ4 )

initial condition
max magnitude (N)
si
f (1 + 0.5µi+4 )
i

external-force
magnitude
γi
f (1 + 0.5µi+8 )
i

external-force
frequency
λi
3ω0 (1 + 0.5µi+12 )

Table 7.1
Effect of parameters on truss geometry, initial conditions, and applied forces. Here, f i , i = 1, . . . , 4 denote the nominal
force magnitudes (to be specified within each experiment) and ω0 denotes the lowest-magnitude eigenvalue of the structure at
the nominal point µ̄.

This discretization leads to a model that corresponds to a Lagrangian dynamical system, with configuration manifold Q = RN , Riemannian metric g(v, w; µ) = vT M (µ)w, nonlinear potential-energy function
V and dissipation function F (q̇; µ) ≡ 21 q̇T C (µ)q̇. Here, C (µ) = αM (µ) + β∇qq V (0; µ) corresponds to
Rayleigh damping. Here, α and β are chosen such that the damping ratio is a specified value ζ for the
uncoupled ODEs associated with the smallest two eigenvalues of the matrix pencil (M(µ̄), ∇qq V (0; µ̄)) [10].
To numerically solve the Lagrangian equations of motion in the time interval [0, T] with T = 25 seconds,
we employ the implicit midpoint rule (a symplectic integrator). This ensures that the numerical solution will
yield symplectic time-evolution maps in the conservative case. We employ a globalized Newton solver with a
More–Thuente linesearch [12] to solve the system of nonlinear algebraic equations arising at each time step.
Convergence of the Newton iterations is declared when the residual norm reaches 10−6 of its value computed
using a zero acceleration and the values of the displacement and velocity at the beginning of the timestep.
The linear system arising at each Newton iteration is solved directly.
The experiments compare the performance of four reduced-order models: Galerkin projection (Section
5.1), collocation (Section 5.2.1), and gappy POD (Section 5.2.2), and the proposed structure-preserving
methods. All reduced-order models (ROMs) employ the same POD reduced basis V, which is computed
by applying Algorithm 1 with snapshots of the configuration variables and an energy criterion ηq ∈ [0, 1]
specified within each experiment. The POD bases Wθ employed by the gappy POD approach (see Section
5.2.2) are generated in the same way. In all cases, snapshots are only collected for the first half of the
time interval at the training points; as a result, the second half of the time interval can be considered
predictive—even for the training set.
Reduced-order models with complexity reduction employ the same sampling matrix P, which is generated
using GNAT’s greedy sample-mesh algorithm [7]. These models are also implemented using the samplemesh concept [7]. To solve optimization problems (2.2), we use the Poblano toolbox for unconstrained
optimization [12]. The initial guess for each of these problems is chosen as PT PV. In practice, we always
found the constraints to be inactive at the unconstrained solution to (2.12); therefore, this reduces to a linear
least-squares problem that we solve directly.
To compare the performance of the reduced-order models, we will consider the response quantity of
interest to be the y-displacement of the bottom-left node of the end face of the truss in Figure 7.1; we denote
this (parameterized, time-dependent) quantity by y ∈ R. The reported errors will be a normalized 1-norm
(in time) of the error in this quantity:
P
|yROM (t; µ? ) − yHFM (t; µ? )|
t∈Tsample (µ? )

error =

.


|Tsample

(µ? )|

max

t∈Tsample (µ? )

yHFM

(t; µ? )

−

min

t∈Tsample (µ? )

yHFM

(7.4)

(t; µ? )

Here, yROM denotes the response computed by a reduced-order model, yHFM is the high-fidelity ‘truth’
response, and Tsample (µ? ) ⊂ [0, T] denotes the time instances selected by the time integrator for online point
µ? .5 In addition to the error in Eq. 7.4, we will compare the speedup achieved by the reduced-order models,
measured as the ratio of the reduced-order-model simulation time to the full-order-model simulation time.
All computations are carried out in Matlab on a Mac Pro with 2 × 2.93 GHz 6-Core Intel Xeon processors
and 64 GB of memory.
5 We employ this error measure because it is insensitive to shifts in the average value of the displacement, unlike other
measures such as the average 1-norm.

19

7.1. Conservative case. We first consider the conservative case characterized by zero damping ζ = 0
and no external forces µi = −2 for i = 9, . . . , 16. This scenario is particularly interesting, as the full-order
model corresponds to a conservative Lagrangian dynamical system characterized by energy conservation and
symplectic time-evolution maps, and our method can also be interpreted as preserving Hamiltonian structure
(see Appendix A). Because we numerically solve the equations of motion using the implicit midpoint rule,
which is a symplectic integrator, the numerical solution is also characterized by a sympletic time-evolution
map. This will also hold for reduced-order models that preserve Lagrangian structure, i.e., the Galerkin
reduced-order model and the two proposed techniques (see Table 6.1). Note also that the dynamics of
undamped, unforced structures are typically quite stiff, which often leaves reduced-order models prone to
instabilities. As a result, we are free to vary parameters µi , i = 1, . . . , 8, which affect only the geometry and
2
initial condition. We set the nominal forces that affect the initial condition to f 1 = f 2 = 2kg × 9.81m/s
2
and f 3 = f 4 = 0.4kg × 9.81m/s .
We first perform a timestep-verification study for the nominal point µ̄ characterized by µ̄i = 0, i = 1, . . . 8
to ensure we employ an appropriate timestep in the numerical experiments. Results are shown in Figure
7.2. A timestep size of ∆t = 0.008 seconds yields an observed convergence rate in the time-averaged tip
displacement of 1.98, which is close to the asymptotic rate of convergence of the implicit midpoint rule, and
an approximated error in the time-averaged tip displacement using Richardson extrapolation of 5.16 × 10−7 .
We can therefore declare this to be an appropriate timestep size for the numerical experiments. Further, we
note that the average number of Newton iterations per timestep is 3.15, so the geometric nonlinearity in the
potential-energy function is significant.

3

∆t = 0.002
∆t = 0.004
∆t = 0.008

10−6
relative error in time-averaged tip displacement

4

tip displacement (m)

2
1
0

10−7

−1
−2
−3

−4
0

5

15
10
time (seconds)

20

10−8 −3
10

25

∆t

10−2

Fig. 7.2. Conservative case: Timestep verification study at the nominal point. A timestep of ∆t = 0.008 seconds gives a
convergence rate of 1.98 and approximate relative error of 5.16 × 10−7 . Note that the three responses in the left figure nearly
overlap.

7.1.1. Fixed parameters. We now test the model-reduction techniques in a fixed-parameters scenario.
That is, we employ the nominal point in the parameter space for both the training and online points: {µi } = µ̄
and µ? = µ̄. Recall that we only collect snapshots for the first half of the time interval, so the second half
can be considered a predictive regime. Note that the two proposed structure-preserving methods are the
same for this case: they both exactly approximate the reduced mass matrix when the parameters are fixed.
The POD reduced basis V is generated using an energy criterion of ηq = 1 − 10−5 in Algorithm 1 of
Appendix B; this leads to a basis dimension of only n = 11  N . The gappy POD-based reduced-order
model employs an energy criterion of 1 (i.e., no truncation) for its reduced bases Wθ (see Section 5.2.2).
Figures 7.3 and 7.4 report results for the reduced-order models as the number of sample indices varies.6
6 In all response plots, a ‘flat line’ indicates that the nonlinear solver failed to converge after 500 Newton iterations at three
different time steps.

20

First, note that the Galerkin reduced-order model is accurate (relative error of 5.42%); however, it
yields a speedup of only 1.69. This is to be expected, as it preserves Lagrangian structure, but has no
complexity-reduction mechanism (see Section 5.1.1). In addition, the proposed reduced-order model—which
also preserves structure, yet has a complexity-reduction mechanism—yields a stable and accurate response
regardless of the number of sample points chosen. For example, 0.4% sampling yields a relative error of 7.3%
and a speedup of 207.0. Sampling 2% of the indices yields an error of 0.71% and a speedup of 34.5, and
sampling 5% of the indices leads to 0.48% error and a speedup of 15.7. Note that sampling beyond 5% does
not improve the method’s accuracy; however, it degrades the speedup, as it requires computing more entries
of the vector-valued functions.
The other complexity-reducing reduced-order models (gappy POD and collocation) are always unstable
except for collocation in the full-sampling case, when it is equivalent to Galerkin. This clearly highlights the
practical benefits of preserving structure in model reduction, as existing structure-destroying complexityreduction methods failed, even in the relatively simple scenario of fixed parameter values.7

2

0

−2

−4

0

10

20

4

tip y-displacement (m)

4

tip y-displacement (m)

tip y-displacement (m)

4

2

0

−2

−4

0

10

20

2

0

−2

−4

0

10

time (seconds)

time (seconds)

time (seconds)

(a) 0.4% sampling

(b) 2% sampling

(c) 100% sampling

20

Fig. 7.3. Conservative, fixed-parameters case: reduced-order model responses as a function of sampling percentage m/N ×
100%. Legend: full-order model (black), Galerkin ROM (dark blue), structure-preserving ROM (light blue), gappy POD ROM
(red), collocation ROM (green), end of training time interval (black circle).

7.1.2. Varying parameters. We now consider a fully predictive scenario with µ? 6∈ {µi }. We use p = 6
training points and determine {µi } using Latin hypercube sampling. The online points are subsequently
chosen randomly in the parameter space. Figure 7.5 depicts the tip displacement for the training points.
Note that the responses are significantly different from one another. The two proposed structure-preserving
reduced-order models will now be different from one another, as the parameters are varying, which means
that the parameterized mass matrix will be approximated differently by the two techniques (see Table 6.1).
The reduced-order models employ a POD reduced basis with a truncation energy criterion of ηq =
1 − 10−6 , which yields a basis dimension of n = 147  N . Again, the gappy POD-based reduced-order
model employs a truncation criterion of ηθ = 1 for its reduced bases. Figure 7.6 reports the tip displacements
generated by the reduced-order models for the three randomly chosen online points, and Figure 7.7 reports
the speedup and errors achieved by the reduced-order models as a function of the number of sample indices.
Again, note that the Galerkin reduced-order model is stable and accurate, as it generates relative errors
of 18.7%, 14.5%, and 9.16% at the three online points, respectively. However, it yields discouraging speedups
of 0.81 (i.e., the simulation was slower than for the full-order model), 1.61, and 1.32 at these points. The
proposed structure-preserving methods are always stable and quite accurate. They yield nearly the same
performance, although method two (which employs the matrix gappy POD approximation) generates lower
errors for online points with 4.9% sampling. From Figure 7.6, note that the high-frequency oscillations that
characterize the proposed methods’ responses are smoothed out when the sampling percentage reaches 20%.
7 We will show in Section 7.2 that introducing dissipation improves the performance of both collocation and gappy POD.
Note that gappy POD was also unstable for other attempted energy criteria of ηθ = 1 − 10−9 and ηθ = 1 − 10−8 .

21

103

Galerkin
gappy POD
collocation
structure pres

relative error

102
speedup

100

Galerkin
gappy POD
collocation
structure pres

10−1

101
100

10−2
10−1
0

20

40
60
80
sampling percentage

0

100

20

40
60
80
sampling percentage

100

Fig. 7.4. Conservative, fixed-parameters case: reduced-order model performance as a function of sampling percentage
m/N × 100%. Missing data points for gappy POD and collocation ROMs indicate unstable responses.
5

tip y-displacement (m)

4
3
2
1

µ12
µ3
µ4
µ5
µ6
µ

0

−1
−2
−3

−4
−5

5

10

15

20

25

time (seconds)
Fig. 7.5. Conservative, varying-parameters case: tip displacement for the training set {µi }.

In particular, proposed method 2 generates speedups of 15.9, 28.5, and 26.2 and relative errors of 11.6%,
13.0%, and 11.6% for 4.9% sampling. For 20% sampling, the method generates speedups of 4.84, 9.82, and
7.72, and relative errors of 1.51%, 5.83%, and 1.09%.
In this example, the gappy POD reduced-order model is unstable for all sampling percentages, and
the collocation reduced-order model is only stable for 100% sampling (at which point it is mathematically
equivalent to the Galerkin reduced-order model). This is not surprising, as these methods do not preserve
problem structure, nor do they guarantee energy conservation. This poor performance can be attributed to
the stiff dynamics that characterize the considered conservative Lagrangian dynamical system, which lead
to instabilities for both reduced-order models.
This example strongly showcases the practical importance of preserving Lagrangian structure: the proposed structure-preserving reduced-order models are the only models that yield both fast and accurate
results.
7.2. Non-conservative case. We now consider the non-conservative case in which the nonconservative dissipative and external forces are nonzero. That is, we set ζ = sin(5◦ ) and all parameters
2
µi , i = 1, . . . , 16 are free to vary. We again set the nominal forces to f¯1 = f¯2 = 2kg × 9.81m/s and
2
f¯3 = f¯4 = 0.4kg × 9.81m/s .
As before, we perform a timestep-verification study for the nominal point µ̄ characterized by µ̄i = 0, i =
22

4

5

0

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2
2

0

−2

0

−2

0

10

−4

20

0

time (seconds)

10

−5

20

0

time (seconds)

10

20

time (seconds)

(a) 4.9% sampling
4

5

0

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2
2

0

−2

0

−2

0

10

−4

20

0

time (seconds)

10

−5

20

0

time (seconds)

10

20

time (seconds)

(b) 20% sampling
4

5

0

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2
2

0

−2

0

−2

0

10

time (seconds)

20

−4

0

10

time (seconds)

20

−5

0

10

20

time (seconds)

(c) 100% sampling
Fig. 7.6. Conservative, varying-parameters case: reduced-order model responses as a function of sampling percentage
m/N × 100% for three randomly chosen online points. Legend: full-order model (black), Galerkin ROM (dark blue), structurepreserving ROM method 1 (magenta), structure-preserving ROM method 2 (light blue), gappy POD ROM (red), collocation
ROM (green), end of training time interval (black circle).

1, . . . 16 to discover an appropriate timestep. A timestep size of ∆t = 0.1 seconds leads to an approximated
error using Richardson extrapolation of 1.07 × 10−4 . We can therefore declare this to be an appropriate
timestep size for the numerical experiments. Further, we note that the average number of Newton iterations
per timestep is 2.56, so the nonlinearity remains significant.
7.2.1. Fixed parameters. We again test the different methods in the fixed-parameters case where
{µi } = µ̄ and µ? = µ̄. As above, we only collect snapshots for the first half of the time interval, and the
proposed structure-preserving methods yield the same results.
The POD reduced basis V is generated using an energy criterion of ηq = 1 − 10−5 , which leads to a
23

100

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

101

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

relative error

speedup

102

10−1

100

10−1
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(a) online point 1

100

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

101

100
0

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

relative error

speedup

102

10−1

20

40
60
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(b) online point 2

100

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

101

100
0

Galerkin
gappy POD
collocation
structure pres 1
structure pres 2

relative error

speedup

102

10−1

20

40
60
80
sampling percentage

100

10−2
0

20

40
60
80
sampling percentage

100

(c) online point 3
Fig. 7.7. Conservative, varying-parameters case: reduced-order model performance as a function of sampling percentage
m/N × 100%.
24

basis dimension of n = 6  N . The gappy POD-based reduced-order model employs an energy criterion of
1 for its reduced bases Wθ . Figures 7.8 and 7.9 report results for the reduced-order models as the number
of sample indices varies.
Again, the Galerkin reduced-order model is accurate, with a relative error of 1.57%, but produces a
speedup of only 1.33. The proposed structure-preserving method is always stable as expected. Its performance is dependent upon the sampling percentage, with (arguably) the best performance achieved for
2% sampling (6.28% error and 36.5 speedup). For 0.2% sampling, the method produces 16.1% error and a
speedup of 251; 20% sampling leads to 5.39% error and a speedup of 4.6.
The gappy POD reduced-order model is unstable for 0.2%, 2%, and 5% sampling, but stabilizes at 20%;
compared to the conservative case, this stability likely results from less stiff dynamics due to the presence
of damping. This yields its best performance of 1.53% error, but only a 4.1 speedup.8 The collocation
reduced-order model is stable only for full sampling, when it is equivalent to Galerkin.

0

−2

−4

0

10

20

2

tip y-displacement (m)

2

tip y-displacement (m)

tip y-displacement (m)

2

0

−2

−4

0

10

time (seconds)

−2

−4

20

time (seconds)

(a) 0.2% sampling

0

0

10

20

time (seconds)

(b) 5% sampling

(c) 20% sampling

Fig. 7.8. Non-conservative, fixed-parameters case: reduced-order model responses as a function of sampling percentage
m/N × 100%. Legend: full-order model (black), Galerkin ROM (dark blue), structure-preserving ROM (light blue), gappy POD
ROM (red), collocation ROM (green), end of training time interval (black circle).

100

Galerkin
gappy POD
collocation
structure preserving

102

Galerkin
gappy POD
collocation
structure preserving

relative error

speedup

103

10−1

101

100
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

Fig. 7.9. Non-conservative, fixed-parameters case: reduced-order model performance as a function of sampling percentage
m/N × 100%. Missing data points for gappy POD and collocation ROMs indicate unstable responses.

8 A truncation criterion of 1 yielded the best performance for Gappy POD. For η = 1 − 10−9 , Gappy POD was unstable for
θ
all sampling percentages. It was also unstable for all sampling percentages when it employed an energy criterion of ηθ = 1−10−8 .

25

7.2.2. Varying parameters. We now consider the parameter-varying case where µ? 6∈ {µi }. We again
employ p = 6 training points and determine {µi } using Latin hypercube sampling. We choose the online
points randomly in the parameter space. Figure 7.10(a) shows the tip displacement for the training points;
clearly, the responses are significantly different from one another.

4

2

10

µ12
µ3
µ4
µ5
µ6
µ

tip y-displacement (m)

tip y-displacement (m)

6

0

−2

0

−5

−10

−4

−6

5

µ12
µ3
µ4
µ5
µ6
µ

5

10

15

20

−15

25

time (seconds)

5

10

15

20

25

time (seconds)

(a) Section 7.2 experiments

(b) Section 7.3 experiments (higher nonlinearity)

Fig. 7.10. Non-conservative, parameter-varying case: tip displacement for the training set {µi } for two sets of experiments.

Because we are in a fully predictive scenario, the two proposed structure-preserving reduced-order models
again yield different results. All reduced-order models employ an energy criterion of ηq = 1 − 10−5 , which
leads to a basis dimension of n = 12. We employ ηθ = 1 for the gappy POD reduced-order model.
Figures 7.11 and 7.12 report the results for this predictive study at the online points. At all three points,
Galerkin is accurate (relative errors of 9.8%, 7.5%, and 13.5%), but does not yield significant speedups
(speedups of 1.2, 1.4, and 1.1). As is apparent from the plots, the two proposed structure-preserving methods
yield nearly the same performance. At 0.4% sampling, method 1 produces relative errors of 11.0%, 2.82%,
and 10.3% and speedups of 73.3, 96.3, and 82.3. At 2% sampling, method 1 yields relative errors of 10.9%,
4.38%, and 7.97% and speedups of 19.2, 21.6, and 16.8.
In this example, gappy POD does not stabilize until 40% sampling, at which point the speedup is less
than 1. Thus, gappy POD does not yield performance improvement for this problem. Collocation stabilizes
at 80% sampling, and also fails to generate any performance improvement.
7.3. Effect of nonlinearity. We now aim to characterize the dependence of problem nonlinearity on
the proposed methods’ performances. Recall from Section 3 that the potential-energy approximation is
computed by matching the gradient of the potential energy to first order about the equilibrium configuration q0 (µ? ). In the presence of stronger nonlinearity, we expect the configuration to deviate further from
equilibrium, which should degrade the accuracy of the approximation.
To numerically assess the effect of nonlinearity, we repeat the experiments from Section 7.2.2 using
the same training and online points, but we increase the nominal forces by a factor of 2.5 to f 1 = f 2 =
2
2
5kg × 9.81m/s and f 3 = f 4 = 1kg × 9.81m/s . We first perform a timestep-verification study for the
nominal point µ̄. As expected, a smaller timestep size of ∆t = 0.025 seconds is required, as it corresponds
to an approximated error using Richardson extrapolation of 3.62 × 10−4 .
Figure 7.10(b) displays the tip displacement for the training points. Note that the responses are similar to
those for the previous study (see Figure 7.10(a)), but have larger magnitudes and therefore imply a stronger
geometric nonlinearity. The reduced-order models employ a POD reduced basis of dimension n = 14, which
was obtained by an energy criterion of ηq = 1 − 10−5 ; gappy POD uses ηθ = 1 for its nonlinear-function
bases.
Figures 7.13 and 7.14 report the reduced-order models’ performances for this problem. As in the previous
case, Galerkin is accurate (relative errors of 8.3%, 3.0%, and 10.0% at the online points), but does not
26

4

2

2

0

−2

−4

0

10

0

−2

−4

20

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2

0

time (seconds)

10

0

−2

−4

20

0

time (seconds)

10

20

time (seconds)

(a) predictions at three randomly chosen online points, 0.4% sampling
4

2

0

−2

−4

0

10

0

−2

−4

20

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2
2

0

time (seconds)

10

0

−2

−4

20

0

time (seconds)

10

20

time (seconds)

(b) predictions at three randomly chosen online points, 40% sampling
4

2

0

−2

−4

0

10

20

time (seconds)

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

2
2

0

−2

−4

0

10

time (seconds)

20

0

−2

−4

0

10

20

time (seconds)

(c) predictions at three randomly chosen online points, 80% sampling
Fig. 7.11. Non-conservative, parameter-varying case: reduced-order model responses as a function of sampling percentage
m/N × 100%. Legend: full-order model (black), Galerkin ROM (dark blue), structure-preserving ROM method 1 (magenta),
structure-preserving ROM method 2 (light blue), gappy POD ROM (red), collocation ROM (green), end of training time
interval (black circle).

generate significant speedups (1.67, 1.71, and 1.0). The proposed structure-preserving techniques again yield
very similar results to each other; however, the errors are significantly larger than in the the experiments
from Section 7.2.2 characterized by a less severe nonlinearity. For 0.5% sampling, proposed method 1 yields
relative errors of 21.3%, 11.1%, and 15.9% at the online points and speedups of 116.4, 160, and 98.9. Thus,
increasing the nonlinearity in the problem does have a deleterious effect on the methods’ performances.
However, it is important to note that other complexity-reducing reduced-order models fail to generate
significant performance improvement on this more highly nonlinear problem. In particular collocation is
always unstable for a sampling percentage less than 60%, and gappy POD is always unstable when the
27

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

101

100

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

relative error

speedup

102

10−1

100

10−1
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(a) online point 1

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

speedup

102
101

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

10−1

100
10−1
0

100

relative error

103

20

40
60
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(b) online point 2

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

speedup

102
101

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

10−1

100
10−1
0

100

relative error

103

20

40
60
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(c) online point 3
Fig. 7.12. Non-conservative, parameter-varying case: reduced-order model performance as a function of sampling percentage m/N × 100%.

28

percentage is less than 80%. As a result, the best speedup obtained by either of the methods is only 2.77
(collocation for 60% sampling for online point 2).
5

5

0

−5

−10

0

10

20

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

5

0

−5

−10

0

time (seconds)

10

0

−5

0

20

10

20

time (seconds)

time (seconds)

(a) predictions at three randomly chosen online points, 0.5% sampling
5

5

0

−5

−10

0

10

20

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

5

0

−5

−10

0

time (seconds)

10

0

−5

0

20

10

20

time (seconds)

time (seconds)

(b) predictions at three randomly chosen online points, 60% sampling
5

5

0

−5

−10

0

10

20

time (seconds)

tip y-displacement (m)

tip y-displacement (m)

tip y-displacement (m)

5

0

−5

−10

0

10

time (seconds)

20

0

−5

0

10

20

time (seconds)

(c) predictions at three randomly chosen online points, 80% sampling
Fig. 7.13. Non-conservative, highly nonlinear parameter-varying case: reduced-order model responses as a function of
sampling percentage m/N × 100%. Legend: full-order model (black), Galerkin ROM (dark blue), structure-preserving ROM
method 1 (magenta), structure-preserving ROM method 2 (light blue), gappy POD ROM (red), collocation ROM (green), end
of training time interval (black circle).

8. Conclusions. This paper has presented an efficient structure-preserving model-reduction strategy
applicable to simple mechanical systems. The methodology directly approximates the quantities that define
the problem’s Lagrangian structure and subsequently derives the equations of motion, while ensuring low
online computational cost. The method is distinct from typical model-reduction methods for nonlinear ODEs;
these methods are typically based on collocation and DEIM/gappy POD techniques that approximate the
equations of motion and destroy Lagrangian structure. At the core of the methodology are the reduced29

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

102

100

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

relative error

speedup

103

10−1

101

100
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(a) online point 1

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

102

100

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

relative error

speedup

103

10−1

101

100
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(b) online point 2

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

102

100

Galerkin
gappy POD
collocation
structure preserving 1
structure preserving 2

relative error

speedup

103

10−1

101

100
0

20

60
40
80
sampling percentage

10−2
0

100

20

40
60
80
sampling percentage

100

(c) online point 3
Fig. 7.14. Non-conservative, highly nonlinear parameter-varying case: reduced-order model performance as a function of
sampling percentage m/N × 100%.

30

basis sparsification (RBS) and matrix gappy POD techniques for approximating parameterized reduced
matrices while preserving symmetry and positive definiteness; we also employed the former method to
preserve potential-energy structure.
Numerical experiments on a geometrically nonlinear parameterized truss structure highlight the method’s
benefits: preserving Lagrangian structure ensured the method always generated stable responses that were
often very accurate. Other model-reduction techniques were often unstable; achieving stability usually
required too many sample indices to lead to significant performance gains for those methods. The experiments
also showed that both RBS and matrix gappy POD led to nearly the same performance across a range of
experiments.
Future work includes devising a method to improve the method’s robustness in the presence of strong
nonlinearity (e.g., by non-local approximation of the potential-energy function), applying the method to a
truly large-scale problem, devising a technique-specific method for choosing the sample indices, and deriving
error bounds and error estimates that rigorously assess the accuracy of the method’s predictions. Finally,
the RBS and matrix gappy POD methods are relevant to a wider class of problems than model reduction
for Lagrangian systems; future work will investigate to their applicability to other scenarios.
Appendix A. Hamiltonian dynamics. When the Hamiltonian formulation of classical mechanics is
taken, the proposed structure-preserving reduced-order models also preserve problem structure. For simplicity, we consider conservative systems with no dissipation or applied external forces. The Hamiltonian
ingredients for conservative simple mechanical systems are then the same as those for Lagrangian dynamical
systems as described in Section 1:
• A differentiable configuration manifold Q, which we set to Q = RN .
• A parameterized Riemannian metric g(v, w; µ), which we set to g(v, w; µ) = vT M (µ)w, where
M (µ) denotes the N × N parameterized symmetric positive-definite mass matrix.
• A parameterized potential-energy function V : Q × D → R.
Again, the kinetic energy can be expressed as T (q̇; µ) = 21 g(q̇, q̇; µ) = 12 q̇T M (µ)q̇, and the Lagrangian
becomes L(q, q̇; µ) = 12 q̇T M (µ)q̇ − V (q; µ).
The conjugate momenta p : [0, T] × D → RN can then be derived as
∂L
(t; µ)
∂ q̇
= M (µ)q̇(t; µ).

p(t; µ) =

(A.1)
(A.2)
−1

In terms of the conjugate momenta, the kinetic energy then becomes Tq (p; µ) = pT M (µ) p. By definition,
the Hamiltonian H : RN × RN × D → R is the Legendre transformation of the Lagrangian function:
H(p, q; µ) = q̇T p − L(q, q̇, t)
1
−1
= pT M (µ) p + V (q; µ).
2

(A.3)
(A.4)

Equivalently, it is the total energy, or sum of the kinetic and potential energies for classical mechanical
systems. The equations of motion can then be obtained by applying Hamilton’s equations of motion
∂H
∂p
∂H
ṗ = −
,
∂q
q̇ =

(A.5)
(A.6)

which are equivalent to
q̇ = M (µ)q

(A.7)

M (µ)q̈ + ∇q V (q; µ) = 0

(A.8)

from the definition of the Hamiltonian and conjugate momenta. Note that these are equivalent to the
equations of motion for a simple mechanical system (4.5) in the conservative case derived using the Lagrangian
formalism.
31

A.1. Galerkin reduced-order model. To construct the structure-preserving Galerkin reduced-order
model for the Hamiltonian formalism, we follow this same recipe, but with the reduced ingredients defined
previously. First, we define the reduced configuration space Qr ∈ Rn with Qr ≡ {q̄(µ) + Vqr | qr ∈ Qr },
and subsequently the reduced Lagrangian from Eqs. (5.5)–(5.6)
Lr (qr , q˙r ; µ) ≡ L(q̄(µ) + Vqr , Vq˙r ; µ)
1
= q˙r T VT M (µ)Vq˙r − V (q̄(µ) + Vqr ; µ).
2

(A.9)
(A.10)

Then, the reduced conjugate momenta pr : [0, T] × D → Rn can be computed as
∂Lr
(t; µ)
∂ q˙r
= VT M (µ)Vq˙r (t; µ),

pr (t; µ) =

(A.11)
(A.12)

and the reduced Hamiltonian is, by definition,
Hr (pr , qr ; µ) = q˙r T pr − L(qr , q˙r , t)
−1
1
pr + V (q̄(µ) + Vqr ; µ).
= pr T VT M (µ)V
2

(A.13)
(A.14)

Applying Hamilton’s equations of motion then yields
∂Hr
∂pr
∂Hr
,
p˙r = −
∂qr
q˙r =

(A.15)
(A.16)

or equivalently
q˙r = VT M (µ)Vqr
T

T

V M (µ)Vq¨r + V ∇q V (q̄(µ) + Vqr ; µ) = 0

(A.17)
(A.18)

Again, these are equivalent to the Galerkin reduced-order equations of motion (5.11) obtained from the
Lagrangian perspective in the absence of non-conservative forces. Clearly, this formulation for a reducedorder model preserves Hamiltonian structure, as the equations of motion are consistent with the Hamiltonian
formalism.
A.2. Proposed structure-preserving methods. The proposed methods preserve Hamiltonian structure. The derivation follows that of the Galerkin reduced-order model in the previous section, but with the
reduced Lagrangian replaced by the approximated reduced Lagrangian from Eq. (A.10):
1
g̃r (q˙r , q˙r ; µ) − Ṽr (qr ; µ)
2
1
= q˙r T M̃ (µ) q˙r − Ṽr (qr ; µ),
2

L̃r (qr , q˙r ; µ) =

(A.19)
(A.20)

where M̃ (µ) is an n × n symmetric positive definite matrix generated by one of the methods presented in
Section 2, and Ṽr (qr ; µ) is an approximated potential approximated according to the method outlined in
Section 3.
The approximated reduced conjugate momenta are then
p˜r (t; µ) =

∂ L̃r
(t; µ)
∂ q˙r

= M̃ (µ) q˙r (t; µ),
32

(A.21)
(A.22)

and the approximated reduced Hamilton is
H̃r (p˜r , qr ; µ) = q˙r T p˜r − L(qr , q˙r , t)
1
−1
= p˜r T M̃ (µ) p˜r + Ṽr (qr ; µ).
2

(A.23)
(A.24)

Applying Hamilton’s equations yields the following equations of motion
q˙r = M̃ (µ) qr

(A.25)

M̃ (µ) q¨r + ∇qr Ṽr (qr ; µ) = 0,

(A.26)

which are computationally inexpensive to solve, and also retain Hamiltonian structure, as they are consistent
with the Hamiltonian formalism. Note that again these are equivalent to the equations of motion obtained
with the proposed strategy in the Lagrangian case with non-conservative forces set to zero.
Appendix B. Proper orthogonal decomposition. Algorithm 1 describes the method for computing
a proper orthogonal decomposition (POD) basis given a set of snapshots. The method amounts to computing
the singular value decomposition of the snapshot matrix; the left singular vectors define the POD basis.
Algorithm 1 Proper-orthogonal-decomposition basis computation (normalized snapshots)
x
Input: Set of snapshots X ≡ {xi }ni=1
⊂ RN , energy criterion η ∈ [0, 1]
Output: W (X , η)
1: Compute the thin singular value decomposition X = UΣVT , where X ≡ [x1 /kx1 k · · · xnx /kxnx k].
2: Choose dimension of truncated basis n = ne (η), where

ne (η) ≡ min i

(B.1)

i∈V(η)

V(η) ≡ {n ∈ {1, . . . , nx } |

n
X

σi2 /

i=1

3:

nx
X

σj2 ≥ η},

(B.2)

j=1

and Σ ≡ diag
 (σi ) with σ1 ≥ · · · ≥ σnx ≥ 0.

W (X , η) = u1 · · · un , where U ≡ u1 · · · unx .

Appendix C. Solving the matrix gappy POD optimization problem. This approach reformulates the constraints of problem (2.11) in terms of eigenvalues of the reduced matrix. That is, problem (2.11)
is reformulated as
nA
X

T

minimize
kP A (µ)P −
PT Ak Pxk k2F
x≡(x1 ,...,xnA )
k=1
subject to λ̃j (x) ≥  > 0,

(C.1)

j = 1, . . . , n.

Here, λ̃j (x), j = 1, . . . , n are the eigenvalues of the low-dimensional matrix

n
A
P

VT Ai Vxi and  denotes a

i=1

numerical threshold for defining a full-rank matrix. This problem can be numerically solved, e.g., using a
gradient-based algorithm.
The gradient of the quadratic objective function is obvious. The gradient of the constraint can be derived
by assuming distinct eigenvalues:
n

A
P
∂
VT Ak Vxk
∂ λ̃j
k=1
= ỹjT
ỹj
(C.2)
∂xi
∂xi

= ỹjT VT Ai V ỹj .
(C.3)
33

∂ λ̃

Here, ỹj is the eigenvector associated with eigenvalue λ̃j . This indicates that computing the gradient ∂xji is
inexpensive and requires the following steps:
n
A
P
VT Ak Vxk .
1. Compute the eigenvector ỹj ∈ Rn of the matrix
k=1

2. Compute the low-dimensional matrix–vector product w = VT Ai V ỹj .
3. Compute the low-dimensional vector–vector product ỹjT w.
We propose using the unconstrained solution to problem (C.1) as the initial guess. In practice, this
solution is often feasible, so it is typically unnecessary to handle the constraints directly. In the rare cases
where it is necessary to deal with multiple equal eigenvalues—or a number of nearby eigenvalues—the
methods presented by Andrew and Tan [3] can be used to produce a numerically stable gradient of the
constraint; this was not required in the numerical experiments reported in Section 7.
Appendix D. Proofs.
D.1. Proof of Theorem 2.1. The proof relies on a generalization of the well known Cauchy interlacing
thereom. To prove this generalization, we use a theorem (Theorem 4.3.10) from Ref. [16] that we restate
below.
(s)
(r)
Theorem D.1. Let two sequences of interlacing real numbers be given by (λi )ni=1 and (λi )m
i=1 as
(r)
(s)
(r)
(s)
described by inequality (2.5) when m = n + 1. Define Λ = diag(λi ) and Λ = diag(λi ). Then, there
exists a real number α ∈ R and a vector y ∈ Rn such that Λ(s) are the eigenvalues of the real symmetric
matrix
 (r)

Λ
y
B̂(bordered) ≡
.
yT α
The following corollary is a direct consequence of the above theorem.
Corollary D.2. Given B̂(s) ∈ SPD (m) and B̂(r) ∈ SPD (m − 1), where SPD (k) denotes the set of
k × k symmetric positive-definite matrices, whose eigenvalues interlace, then
∃ Um such that UTm B̂(s) Um = B̂(r)

UTm Um = I.

with

(D.1)

Proof. Using the above theorem, a matrix B̂(bordered) ∈ SPD (m) exists that shares the same eigenvalues
with B̂(s) . Let Q(bordered) , Q(s) , and Q(r) be the (square, orthogonal) matrices of eigenvectors for B̂(bordered) ,
B̂(s) , and B̂(r) , respectively. Then,


Q(bordered)

T

B̂(bordered) Q(bordered) = Λ(s) = Q(s) B̂(s) (Q(s) )T ,

which implies that
B̂(bordered) = Q(bordered) (Q(s) )T B̂(s) Q(s) (Q(bordered) )T .
From the definition of B̂(bordered) it also follows that
[I

0]B̂(bordered) [I

0]T = Λ(r) ,

where I is the (m − 1) × (m − 1) identity matrix and 0 is the zero column vector of length m − 1, and thus
h
i
B̂(r) = Q(r) Λ(r) (Q(r) )T = Q(r) [I 0]B̂(bordered) [I 0]T (Q(r) )T .
(D.2)
Combining the above, we can write
B̂(r) = Q(r) [I

0]Q(bordered) (Q(s) )T B̂(s) Q(s) (Q(bordered) )T [I


and so (D.1) is satisfied taking Um = Q(r) [I

0]Q(bordered) (Q(s) )T
34

T

.

0]T (Q(r) )T

The generalization of the Cauchy interlacing thereom now follows.
Theorem D.3. Given two matrices B̂(s) ∈ SPD (m) and B̂(r) ∈ SPD (n) (with m ≥ n), then
UT B̂(s) U = B̂(r)

∃ U such that

with

(r)

UT U = I

(D.3)

(s)

if and only if the eigenvalues λi , i = 1, . . . , n interlace the eigenvalues λi , i = 1, . . . , m defined as
(r)

= λi x̂i ,

(s)
B̂(s) x̂i

(s) (s)
λi x̂i ,

B̂(r) x̂i

(r) (r)

=

i = 1, . . . , n

(D.4)

i = 1, . . . , m.

(D.5)

The definition of interlacing is given by inequality (2.5).
Proof. It is well known that if B̂(s) ∈ SPD (m) is given along with an orthogonal m × n matrix U (with
m ≥ n), then the eigenvalues of UT B̂(s) U must interlace those of B̂(s) . This is referred to as the Cauchy
interlacing theorem (e.g., see [20]).
The converse of the Cauchy interlacing theorem is less widely known. The case m = n follows trivially
using an eigenvalue decomposition. The case m = n + 1 corresponds to the above corollary. The proof is
completed by generalizing the corollary to the m > n + 1 case. This follows from an inductive argument
where one considers a projection that reduces the matrix dimension of B̂(s) by one. According to the
above corollary, we have a great deal of flexibility in choosing this lower dimensional matrix if its eigenvalues
interlace those of the higher dimension matrix. We then choose a lower dimensional matrix whose eigenvalues
not only interlace those of B̂(s) but whose eigenvalues are also interlaced by those of B̂(r) . That is,
(s)

λi

(s)

≤ µi ≤ λi+1

(r)

and µi ≤ λi

≤ µi+m−n−1 ,

where µi denotes the ith smallest eigenvalue of the intermediate matrix. Rewriting this we obtain the
following intervals for the eigenvalues µi :
(
(s)
(r)
max(λi , λi−m+n+1 ) i ≥ m − n
µi ≥
(s)
λi
i<m−n
and
(
µi ≤

(s)

(r)

min(λi+1 , λi ) i ≤ n
.
(s)
λi+1
i>n

Using the interlacing property for B̂(s) and B̂(r) , one can verify that the intervals for the µi are nonempty.
(s)
(s)
(r)
(s)
(r)
(r)
That is, λi+1 ≥ λi and for those i such that λi−m+n+1 is defined, we have λi+1 ≥ λi−m+n+1 , λi ≥
(s)

(r)

(r)

λi , and λi ≥ λi−m+n+1 Thus, there exists an orthogonal matrix Um such that the (m−1) × (m−1)
matrix UTm B̂(s) Um has eigenvalues that are interlaced by those of B̂(r) . We repeat this procedure each time
reducing the matrix dimension by one until the final reduction where we take the lower dimension matrix to
be equal to B̂(r) . This implies that there exists a set of projection matrices such that UT B̂(s) U is equal to
B̂(r) where U = Um Um−1 ...Un+1 .
Equipped with the generalized Cauchy interlacing theorem, we now prove the exactness condition for
the A (µ) term which is restated in slightly simplified notation.
Theorem D.4. Let A (µ) have the form
A (µ) = h1 (µ) A1 + h2 (µ) A2

(D.6)

where A1 ∈ SPD (N ), A2 ∈ SPD (N ), and h1 , h2 : D → R.
Then,
∃ UA such that

UA T PT A (µ) P UA = VT A (µ) V,
35

∀µ ∈ D

(D.7)

if and only if the eigenvalues of the general matrix pencil
(r)

= λi D(r) xi ,

(s)

= λi D(s) xi ,

B(r) xi

(r)

(r)

i = 1, . . . , n

(D.8)

(s)

(s)

i = 1, . . . , m

(D.9)

interlace the eigenvalues of
B(s) xi
where
D(r)

=

B(r)

=




V T A1 V ,
 T

V A2 V ,

and

D(s)

=

B(s)

=

 T

P A1 P ,
 T

P A2 P .

The definition of interlacing is given by
(s)

λi

(r)

≤ λi

(s)

≤ λi+m−n

f or i = 1, ..., n

(D.10)

where the eigenvalues are indexed in order of increasing magnitude.
Proof. Clearly (D.7) can only hold for any µ ∈ D and any functions h1 (µ) and h2 (µ) if and only if
UA T D(s) UA = D(r) and UA T B(s) UA = B(r) .

(D.11)

Using a carefully chosen linear transformation, it follows that proving the theorem is equivalent to proving
the following:
∃ U such that UT B̂(s) U = B̂(r)
(r)

if and only if the eigenvalues λi
in Eqs. (D.8)–(D.9) also satisfy

with

UT U = I

(D.12)

(s)

interlace the eigenvalues of λi , where the eigenvalues previously defined
B̂(r) x̂(r) = λ(r) x̂(r)

(D.13)

B̂(s) x̂(s) = λ(s) x̂(s) .

The linear transformation relies on Cholesky factorizations given by D(s) = L(s) (L(s) )T and D(r) =
L(r) (L(r) )T . These factors lead to the following definitions
B̂(s)

=

(L(s) )−1 B(s) (L(s) )−T ,

x̂(s)

=

(L(s) )T x(s) ,

B̂(r)

=

(L(r) )−1 B(r) (L(r) )−T , x̂(r)

=

(L(r) )T x(r) , and

U

=

(L(s) )T UA (L(r) )−T ,

which can be used in Eqs. (D.8), (D.9) and (D.11) to obtain Eqs. (D.12) and (D.13). The proof is completed
by recognizing that Eqs. (D.12) and (D.13) correspond to the generalized Cauchy interlacing thereorm.
D.2. Proof of Theorem 2.2. Proof. If condition 2 holds, then the unconstrained solution to problem
(2.11) is
+
ξA (µ? ) = P̄T Wa P̄T v (A (µ? )) .
(D.14)
If condition 1 holds, then the vectorized matrix can be expressed as
v (A (µ? )) = Wa z (µ? ) ,

(D.15)

or equivalently
A (µ? ) =

nA
X
i=1

36

z i (µ? ) Ai ,

(D.16)


where z ≡ z 1 , . . . , z nA . Substituting Eq. (D.15) into Eq. (D.14) gives ξA (µ? ) = z (µ? ) and so Eq. (2.6)
yields
?

Ã (µ ) =

nA
X

z i (µ? ) VT Ai V.

(D.17)

i=1

Comparing Eqs. (D.17) and (D.16) gives the exactness result: Ã (µ? ) = VT A (µ? )VT .
T

g
g
g
D.3. Proof of Lemma 1. Proof. The equation U
V UV = I simply states that the m columns of UV
g
are orthogonal and so any orthogonal matrix U
V satisfies the first part of (3.8). Thus, solvability amounts
T
T
g
g T
to finding an orthogonal matrix U
V such that Ṽ c̃ = UV P c̃. For a solution to exist, however, it is
T

T
T
g
obviously necessary that ||ṼT c̃||2 = ||U
V P c̃||2 . If the vector P c̃ lies within the span of the basis defined
T
T
T
T
g
g T
g T
by the columns of U
V , then UV P c̃ preserves its 2-norm and so ||UV P c̃||2 = ||P c̃||2 . That is,
T
g
application of U
corresponds to a rotation of PT c̃ about the origin and so length is preserved. If instead
V
T

T
T
g
the vector PT c̃ lies only partially within the span of the orthogonal basis, then ||U
V P c̃||2 < ||P c̃||2 .
T
g
That is, application of U
corresponds to a rotation of the component of PT c̃ lying within the span of the
V
orthogonal basis. This implies that a necessary condition for a solution to (3.8) is that

||ṼT c̃||2 ≤ ||PT c̃||2 .

(D.18)

Case 1: m = n
g
PT c̃ must lie within the range of U
V (as it is a full rank square matrix) and so it is necessary to have
n×n
g
equality in (D.18) when m = n. One possible U
V in this case is obtained by first defining a Q1 ∈ R
n×n
T
T
and a Q2 ∈ R
such that the first row of Q1 is α1 c̃ Ṽ with α1 = 1/kc̃ Ṽk2 . Likewise, the first row of
Q2 is taken as α1 c̃T P; note that α1 also normalizes this row because we assume kṼT c̃k2 = kPT c̃k2 . All
remaining rows are chosen so that both Q1 and Q2 are orthogonal matrices. This gives
Q1 ṼT c̃ = ||ṼT c̃||2 e1 = ||PT c̃||2 e1 = Q2 PT c̃,
where e1 is the first canonical unit vector of length n (first element is one and all other n − 1 components
T
g
g = QT Q2 . Thus, equality in (D.18) is also
are zero). A suitable U
V that solves (3.8) is then given by UV
1
sufficient when m = n.
Case 2: m > n
g
g
The matrix U
V is now rectangular. One possible UV is obtained by defining Q1 as before while instead
defining an m × m orthogonal matrix Qfull with the first row again set to α2 c̃T P with α2 = 1/kPT c̃k2 . This
gives
Q1 ṼT c̃ = ||ṼT c̃||2 e1 and Qfull PT c̃ = ||PT c̃||2 ẽ1
where ẽ1 ∈ Rm×1 is the first canonical unit vector of length m. If ||PT c̃||2 = ||ṼT c̃||2 , then a suitable
Ṽ solving (3.8) is given by taking Q2 to be the first n rows of Qfull (as Q1 ṼT c̃ = Q2 PT c̃) and taking
T
g
U
= QT1 Q2 . If ||PT c̃||2 > ||ṼT c̃||2 , then we define a vector y as an arbitrary linear combination of the
V
last m − n rows of Qfull such that y has unit norm. The first row of Q2 is then taken as
(Q2 )1 = α3

PT c̃
+
||PT c̃||2

q

1 − α32 y,

where (Q)k denotes the kth row of a matrix Q and α3 = ||ṼT c̃||2 /||PT c̃||2 . The remaining rows of Q2 are
simply (Qfull )k for k = 2, . . . , n. It is easy to verify that Q2 is again orthogonal and that Q2 PT c̃ = ||ṼT c̃||2 e1 .
37

T

g
Thus, Q2 PT c̃ = ||ṼT c̃||2 e1 and U
= QT1 Q2 is a possible solution implying that ||PT c̃||2 ≥ ||ṼT c̃||2 is a
V
necessary and sufficient condition when m > n.
Appendix E. The authors acknowledge Julien Cortial for both insightful discussions and for providing
the original nonlinear-truss code that was modified to generate the numerical results. The authors also
acknowledge Clancey Rowley for useful comments received at the 2013 SIAM Conference on Computational
Science and Engineering.
This research was supported in part by an appointment to the Sandia National Laboratories Truman
Fellowship in National Security Science and Engineering, sponsored by Sandia Corporation (a wholly owned
subsidiary of Lockheed Martin Corporation) as Operator of Sandia National Laboratories under its U.S.
Department of Energy Contract No. DE-AC04-94AL85000. The authors also acknowledge support by the
Department of Energy Office of Advanced Scientific Computing Research under contract 10-014804.
REFERENCES
[1] S. Adhikari, Damping models for structural vibration, PhD thesis, Cambridge University, September 2000.
[2] S.S. An, T. Kim, and D.L. James, Optimizing cubature for efficient integration of subspace deformations, ACM Transactions on Graphics (TOG), 27 (2008), p. 165.
[3] Alan L. Andrew and Roger C. E. Tan, Computation of derivatives of repeated eigenvalues and the corresponding
eigenvectors of symmetric matrix pencils, SIAM J. Matrix Anal. Appl., 20 (1998), pp. 78–100.
[4] P. Astrid, S. Weiland, K. Willcox, and T. Backx, Missing point estimation in models described by proper orthogonal
decomposition, IEEE Transactions on Automatic Control, 53 (2008), pp. 2237–2251.
[5] C. Beattie and S. Gugercin, Structure-preserving model reduction for nonlinear port-hamiltonian systems, in Decision
and Control and European Control Conference (CDC-ECC), 2011 50th IEEE Conference on, IEEE, 2011, pp. 6564–
6569.
[6] K. Carlberg, C. Bou-Mosleh, and C. Farhat, Efficient non-linear model reduction via a least-squares Petrov–Galerkin
projection and compressive tensor approximations, International Journal for Numerical Methods in Engineering, 86
(2011), pp. 155–181.
[7] K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem, The GNAT method for nonlinear model reduction: effective
implementation and application to computational fluid dynamics and turbulent flows, Journal of Computational
Physics, 242 (2013), pp. 623–647.
[8] K. Carlberg, R. Tuminaro, and P. Boggs, Efficient structure-preserving model reduction for nonlinear mechanical
systems with application to structural dynamics, in AIAA Paper 2012-1969, 53rd AIAA/ASME/ASCE/AHS/ASC
Structures, Structural Dynamics and Materials Conference, Honolulu, Hawaii, April 23–26 2012.
[9] S. Chaturantabut and D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal
on Scientific Computing, 32 (2010), pp. 2737–2764.
[10] I. Chowdhury and S.P. Dasgupta, Computation of Rayleigh damping coefficients for large systems, The Electronic
Journal of Geotechnical Engineering, 8 (2003).
[11] M. Drohmann, B. Haasdonk, and M. Ohlberger, Reduced basis approximation for nonlinear parametrized evolution
equations based on empirical operator interpolation, SIAM Journal on Scientific Computing, 34 (2012), pp. A937–
A969.
[12] D.M. Dunlavy, T.G. Kolda, and E. Acar, Poblano v 1.0: A Matlab toolbox for gradient-based optimization, Sandia
National Laboratories, Albuquerque, NM and Livermore, CA, Tech. Rep. SAND, 1422 (2010).
[13] R. Everson and L. Sirovich, Karhunen–Loève procedure for gappy data, Journal of the Optical Society of America A,
12 (1995), pp. 1657–1664.
[14] D. Galbally, K. Fidkowski, K. Willcox, and O. Ghattas, Non-linear model reduction for uncertainty quantification in
large-scale inverse problems, International Journal for Numerical Methods in Engineering, 81 (2009), pp. 1581–1608.
[15] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary
differential equations, vol. 31, Springer Verlag, 2006.
[16] R.A. Horn and C.R. Johnson, Matrix Analysis, Cambridge University Press, 1990.
[17] K. Ito and SS Ravindran, A reduced basis method for control problems governed by PDEs, Control and estimation of
distributed parameter systems, (1998), pp. 153–168.
[18] S. Lall, P. Krysl, and J.E. Marsden, Structure-preserving model reduction for mechanical systems, Physica D: Nonlinear Phenomena, 184 (2003), pp. 304–318.
[19] J.E. Marsden and M. West, Discrete mechanics and variational integrators, Acta Numerica, 10 (2001), p. 514.
[20] Beresford N Parlett, The symmetric eigenvalue problem, vol. 7, SIAM, 1980.
[21] C. Prud’homme, D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici, Reliable realtime solution of parameterized partial differential equations: Reduced-basis output bound methods, Journal of Fluids
Engineering, 124 (2002), pp. 70–80.
[22] C. Prud’Homme, Dimitrios V Rovas, Karen Veroy, Luc Machiels, Yvon Maday, Anthony T Patera, Gabriel
Turinici, et al., Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound
methods, Journal of Fluids Engineering, 124 (2001), pp. 70–80.
38

[23] D. Ryckelynck, A priori hyperreduction method: an adaptive approach, Journal of Computational Physics, 202 (2005),
pp. 346–366.
[24] T. Tonn, Reduced-Basis Method (RBM) for Non-Affine Elliptic Parameterized PDEs, PhD thesis, Ulm University, December 2011.
[25] B. Vandereycken, Low-rank matrix completion by Riemannian optimization, SIAM Journal on Optimization, 23 (2013),
pp. 1214—1236.
[26] D. Wirtz, D. C. Sorensen, and B. Haasdonk, A-posteriori error estimation for DEIM reduced nonlinear dynamical
systems, Preprint Series, Stuttgart Research Centre for Simulation Technology, (2012).

39

