Cross-Gramian-Based Dominant Subspaces

A standard approach for model reduction of linear input-output systems is balanced truncation, which is based on the controllability and observability properties of the underlying system. The related dominant subspace projection model reduction method similarly utilizes these system properties, yet instead of balancing, the associated subspaces are directly conjoined. In this work we extend the dominant subspace approach by computation via the cross Gramian for linear systems, and describe an a-priori error indicator for this method. Furthermore, efficient computation is discussed alongside numerical examples illustrating these findings.


Introduction
Input-output systems map an input function to an output function via a dynamical system. The input excites or perturbs the state of the dynamical system and the output is some transformation of the state. Typically, these input and output functions are low-dimensional while the intermediate dynamical system is high(er)-dimensional. In applications from natural sciences and engineering, the dimensionality of the dynamical system may render the numerical computation of outputs from inputs excessively expensive or at least demanding. Model reduction addresses this computational challenge by algorithms that provide surrogate systems, which approximate the input-output mapping of the original system with a low(er)-dimensional intermediate dynamical system. Practically, the trajectory of the dynamical system's state is constrained to a subspace of the original system's state-space, for example by using truncated projections. A standard approach for projection-based model reduction of inputoutput systems is balanced truncation [25], which transforms the state-space unitarily to a representation that is sorted (balanced) in terms of the input's P r e p r i n t effect on the state (controllability) as well as the state's effect on the output (observability) and discards (truncates) the least important states according to this measure. Instead of balancing, this work investigates a dominant subspaces approach [28], that conjoins instead of balances, the most controllable and most observable subspaces into a projection. This unbalanced model reduction method may yield less accurate reduced-order systems, yet allows a computationally advantageous formulation while also preserving stability and providing an error quantification. The method proposed in this work, combines the dominant subspace projection model reduction (DSPMR) method [28] with the cross Gramian (matrix) [12], which encodes controllability and observability information of an underlying input-output system. For this cross-Gramian-based dominant subspace method, an a-priori error indicator is developed, and the numerical issues arising in the wake of large-scale systems are addressed, specifically by utilizing the hierarchical approximate proper orthogonal decomposition (HAPOD) [17]. Compared to other cross Gramian and SVD model reduction techniques such as [20], the proposed method does not need multiple decompositions, but a single HAPOD. The considered class of input-output systems are generalized linear (time-invariant) systems 1 , mapping input u : R → R M via the state x : R → R N -a solution to an ordinary differential equation -to the output y : R → R Q : with a system matrix A ∈ R N ×N , an input matrix B ∈ R N ×M , an output matrix C ∈ R Q×N and a mass matrix E ∈ R N ×N . In the scope of this work, we assume E to be non-singular as well as the matrix pencil (A, E) to be asymptotically stable, meaning the eigenvalues of the associated generalized eigenproblem lie in the open left half-plane. This type of system arises, for example, in spatial discretizations of partial differential equations using the finite element method. In Section 2 the cross Gramian for generalized linear systems is introduced, followed by Section 3, briefly describing projection-based model reduction, and extending the dominant subspace projection method to the cross Gramian together with an error indicator. The proposed model reduction technique is then tested numerically in Section 4 and a summary is given in Section 5.

Generalized Cross Gramian
In this section, the cross Gramian matrix, introduced in [12], is briefly reviewed from the point of view of generalized linear time-invariant (LTI) systems (1).
Fundamental to system-theoretic model reduction are the controllability and observability operators [1], which are given for (1) by the generalized controllability P r e p r i n t operator C : L 2 → R N and the generalized observability operator O : R N → L 2 : The (generalized) cross Gramian 2 is then defined as a composition of the generalized controllability and observability operators: and jointly quantifies controllability and observability of square systems -systems with the same number of inputs and outputs M = Q. As for linear, square systems with E = I, the generalized cross Gramian is a solution to a matrix equation; in this case a Sylvester-type equation: which can be shown using integration-by-parts of (2): Besides the cross Gramian, the (generalized) controllability Gramian W C := CC * and (generalized) observability Gramian W O := O * O are defined accordingly [31,35]. For systems with a symmetric Hankel operator H := OC, H = H * , the (generalized) cross Gramian has the property: Hence for symmetric systems, either, W X or {W C , W O } can be used interchangeably, if controllability and observability are to be concurrently evaluated. For non-symmetric and especially non-square systems, an approximation to the cross Gramian is defined, based on the column-wise partitioning of the input matrix B and row-wise partitioning of the output matrix C: c q , the non-symmetric cross Gramian [19] for (1) is defined as: which is the cross Gramian of the average system (A,B,C).

Model Reduction
One of the main numerical applications of the cross Gramian is model (order) reduction, which aims to determine lower order surrogate systems for (1), with respect to the state-space dimension N = dim(x(t)). The reduced order model (ROM) with x r : R → R n , n N , has a reduced system matrix A r ∈ R n×n , a reduced input matrix B r ∈ R n×M , a reduced output matrix C r ∈ R Q×n and a reduced mass matrix E r ∈ R n×n , such that the reduced system's outputỹ : R → R Q approximates the full order model's output: in a suitable norm. Following, the projection-based dominant subspace model reduction method is extended to exploit the cross Gramian for computation.

Projection-Based Model Reduction
A commonplace approach to construct reduced order models is mapping the state-space trajectory x(t) to a lower dimensional subspace, using a reduction operator V 1 : R N → R n and a lifting operator U 1 : R n → R N : In the case of (generalized) linear systems (1), the operators U 1 ∈ R N ×n and V 1 ∈ R n×N , can be directly applied to the system components A, B, C and E to obtain the reduced quantities: Hence, the aim is the computation of suitable reducing and lifting operators U 1 , V 1 , which are typically assumed to be bi-orthogonal The dominant subspaces method considered in this work is additionally orthogonal V 1 := U 1 , thus, the reduction process is a Galerkin projection, which is stability preserving, if the symmetric part of the system matrix A is negative definite, and the mass matrix E positive definite [8], This is a generalization of the stability preservation for systems with E = I, mentioned in [28]. If a system does not fulfill (5), a stabilization procedure, see for example [4,Sec. 4], can be applied to the reduced order model. P r e p r i n t

Dominant Subspaces
The Dominant Subspace Projection Model Reduction (DSPMR) is introduced in [28]. The idea behind DSPMR is, instead of balancing controllability and observability Gramians, to combine the associated principal subspaces obtained from approximate system Gramians. This yields a simple model reduction algorithm which is based upon low-rank factors of the controllability and observability Gramians. In [28], a low-rank Cholesky (LR Chol) factor is used, while [21] utilizes singular vectors of a truncated singular value decomposition (tSVD), The controllability and observability subspaces encoded in the matrix factors are now conjoined and orthogonalized, by either a SVD ( [28]) or a (rank-revealing) QR-decomposition ( [21,22]). Either, the left singular vectors U , or the Q factor, can be taken as Galerkin projections, respectively: An extension to the DSPMR method is also proposed in [28], called Refined Dominant Subspace Projection Model Reduction. The eponymous refinement is given by weighting factors ω C , ω O > 0 for the controllability and observability subspace bases respectively. The weighting factors are selected as the Frobenius norm of the respective low-rank factors, Obviously, this is only sensible for the Cholesky factor variant, as the norm of the (orthonormal) singular vectors is one. A similar idea for combining weighted subspaces is also used in the cotangent lift method from [27].

Cross-Gramian-Based Dominant Subspaces
Instead of the controllability and observability Gramians, also the cross Gramian can be used to obtain a dominant subspace projection. An SVD of the cross Gramian, produces left and right singular vectors aggregated in matrices U X and V X , which induce subspaces associated to controllability (U X ) and observability (V X ) of the underlying system (A, B, C).

P r e p r i n t
In [34,Sec. 4.3], it is noted, that the sole use of either, U X or V X , as a Galerkin projection, will largely omit observability or controllability information respectively. Hence, both subspaces should be incorporated in the reducing and lifting operator. Balanced truncation, for example, determines a suitable Petrov-Galerkin projection 3 , where U 1 = V 1 , by simultaneous diagonalization of the controllability and observability Gramians.
For the proposed variant of the dominant subspace method (for an algorithmic description see Section 3.4.1), the left and right singular vectors are conjoined as before, but also scaled column-wise by the singular vector's associated singular values: So, instead of normalizing the controllability and observability subspaces (as a whole), as in refined DSPMR, based on the common controllability-observability measure, the singular values of the cross Gramian, the vectors spanning the compound subspace are scaled individually. Here explicitly a (rank-revealing) SVD is used, instead of a QR decomposition, as the singular values D CO will be used for an error indicator. An advantage of the cross-Gramian-based dominant subspace projection method is this common measure of minimality, the singular values σ i = D CO,ii associated jointly to the "controllability" and "observability" subspaces.

Error Indicator
In this section an error indicator for the cross-Gramian-based dominant subspace method is developed. Previous works, such as [33,39,29,38], already introduced error bounds for the Hardy H 2 -norm. Here, an H 2 -error indicator of simple structure using time-domain quantities is proposed, which is loosely related to the simplified balanced gains approach from [11]. The H 2 -norm is particularly interesting, since an error estimation has relevance for the frequency-domain and the time-domain [37, Ch. 2], since it also describes the energy (L 2 -norm) of the system's impulse response. Before this error indicator is derived, a straightforward property of the matrix exponential is presented.

Lemma 1
Given matrices A ∈ R N ×N and U ∈ R N ×n , n ≤ N , the following holds: Proof. The proof is a trivial consequence on the associativity of the matrix product.
3 Balanced truncation yields a Galerkin projection for state-space symmetric systems, A = A , B = C , E = E [9]. P r e p r i n t Next, the error indicator is constructed, which is derived from the L 2 -norm of the impulse response error system, for which we assume for ease of exposition E = I: We consider only SISO (Single-Input-Single-Output) systems for this error indicator and impulse inputs u(t) = δ(t).
First, the H 2 -norm of the error system, in impulse response form, is transformed in a manner so that Lemma 1 can be applied. Note, that the error system of a SISO system is also a SISO system with a scalar and thus symmetric impulse response: applying the definition of the reduced quantities (4), and subsequently the result of Lemma 1, gives: The next step is approximating the matrix exponentials e AU1U 1 t and e U1U 1 At by the homogeneous system's solution operator, which allows to factor the previous representation to: Forming the cross Gramian by moving the projection error terms (I −U 1 U 1 ) out of the integral, applying von Neumann's trace inequality [24], and replacing the cross Gramian by its SVD (6) yields:

P r e p r i n t
We note that U 1 is the orthogonalized concatenation of U X and V X , and BC is of rank one due to the SISO nature of the system, hence, The last upper estimate results from the relation of the spectral norm and the Frobenius norm ( A 2 ≤ A F ). This derivation yields the following error indicator: The L 2 impulse response model reduction error for a cross-Gramian-based dominant subspaces reduced order model is approximated by: One might assume that the spectral norm representation of the error indicator would be more convenient, yet in Section 4 we will show the advantage of the final Frobenius norm form.

Remark 1
Using the Cauchy-Schwarz inequality as in [14], this impulse response error indicator can be extended to squarely integrable inputs u ∈ L 2 .

Algorithmic and Low-Rank Computation
The computation of the proposed cross-Gramian-based dominant subspace projection, as well as the classic dominant subspace projection consists of two phases: First, the computation of the system Gramians, either the cross Gramian, or the controllability and observability Gramians. And second, the assembly of the reducing (and lifting) operator. For large-scale systems, the computation of dense system Gramians, which are of dimension N × N , may be infeasible or at least inefficient. To this end, lowrank representations of the Gramians can be computed, for the cross Gramian, in example by the implicitly restarted Arnoldi algorithm [34], the factorized iteration [3], a factored ADI [5] or a low-rank empirical cross Gramian [18]. Overall, the cross-Gramian-based dominant subspace algorithm is summarized by: 1. Compute (low-rank) cross Gramian.
2.a Compute SVD of the cross Gramian.

2.b Weight left (controllability) singular vectors with singular values.
2.c Weight right (observability) singular vectors with singular values.

2.d Compute POD of conjoined and weighted left and right singular vectors.
P r e p r i n t The resulting POD modes (left singular vectors) represent the Galerkin projection. In principle, a similar procedure can be conducted using controllability and observability Gramians, yet it is not immediately clear if the SVD of the (weighted) conjoined singular vectors yields in general an equally useful measure.

Numerical Results
Following, two numerical examples are presented to illustrate the previous findings. These numerical experiments are conducted using MATLAB 2018a [23].
The system Gramians needed for the dominant subspace methods, the controllability and observability Gramian for (refined) DSPMR, and the cross Gramian for the cross-Gramian-based dominant subspaces, are computed as empirical Gramians [16] using emgr -empirical Gramian framework in version 5.5 [15]. All simulated trajectories for the construction of these empirical dominant subspaces are computed using the implicit Euler method. Before we present the numerical results, the practical computation of the cross-Gramian-based dominant subspaces is discussed.

Fused Computation
Even for moderately sized systems, the computation of the (cross) Gramian's singular vectors may be a computationally challenging task 4 . To compute the dominant subspace projections from the cross Gramian, or the controllability and observability Gramians, the hierarchical approximate proper orthogonal decomposition (HAPOD) [17] is used. Particularly, the HAPOD allows the computation of the dominant subspace directly, which means the Gramian SVD(s) and the concatenation SVD in one. The HAPOD enables a swift computation of left singular vectors of arbitrary partitioned data sets, based on a selected projection error (on the input data) ε > 0 and a tree hierarchy with the data (Gramian) partitions as leafs. The tree hierarchy utilized for the experiments in this work is given by a combination of special topologies discussed in [17], the incremental HAPOD (maximally unbalanced binary tree) and the distributed HAPOD (star). Two incremental HAPODs are performed for the Gramian partitions respectively and subsequently a distributed HAPOD of the resulting singular vectors from both sub-trees yields the dominant subspace projection. Fig. 1 illustrates the overall HAPOD tree.
Since the HAPOD computes only left singular vectors, but the right singular vectors of the cross Gramian are also needed, the HAPOD of the cross Gramian (left singular vectors) and the transposed cross Gramian (right singular vectors) is computed. In the following numerical examples, the (full-order) empirical linear cross Gramian [16, Sec. 3.1.3] is used, as in-memory storage of the Gramian(s) is possible. For settings, where only parts of the cross Gramian can be kept in memory, the low-rank empirical cross Gramian [18] for the left singular vectors, and the low-rank empirical cross Gramian of the adjoint system for the right P r e p r i n t singular vectors, can be utilized, since the cross Gramian of the adjoint system is equal to the system's transposed cross Gramian: Since a projection-error-driven POD method is used, the following error bound holds 5 , for a given mean projection error ε > 0: This means the error indicator (7) can be bounded using the prescribed cross Gramian's projection error, thus making it an a-priori error indicator.
P r e p r i n t

FOM Benchmark
The first numerical example compares the cross-Gramian-based dominant subspace method with the classic refined dominant subspace method as well as (empirical) balanced truncation 6 for the "FOM" example in [28], which is also part of the SLICOT Benchmark Collection [10]. This linear SISO system (with E = I) of the structure:ẋ is of order N = 1006, and the system components are given by: The empirical Gramians are constructed using impulse input u(t) = δ(t) and the reduced systems are also tested with this input, to evaluate the error indicator. In Fig. 2, the (empirical) cross-Gramian-based dominant subspaces method, the (empirical) refined dominant subspaces method, (empirical) balanced truncation and the error indicator from Section 3.4 are compared, for a given projection error ε ∈ {10 −3 , . . . , 10 −12 }. The same projection error is selected for the controllability and observability Gramians used by the refined DSPMR and low-rank empirical balanced truncation. In Fig. 2a the prescribed projection error of the respective Gramians is plotted against the resulting L 2 model reduction error. For a given projection error the refined DSPMR method produces the lowest model reduction error, and lowrank balanced truncation the largest, while the proposed cross-Gramian-based dominant subspace method is in between. The error indicator overestimates the error for larger and underestimates for smaller projection errors. Note, that the error indicator just scales the projection error by a constant, hence it appears as a line in the log-log plot. Fig. 2b depicts the resulting reduced order of the tested methods against the model reduction error. Balanced truncation produces the smallest, and DSPMR the largest reduced models, again the cross-Gramian-based method is in between. These results follow intuitions that DSPMR produces the most accurate, but largest subspaces, while balanced truncation may have a smaller, but less accurate subspaces. Hence, the cross-Gramian-based dominant subspace method appears as a compromise. The error indicator is rather conservative, which is due to its simple structure.
P r e p r i n t  P r e p r i n t

Convection Benchmark
The second numerical example evaluates the convection benchmark [36, Convection] 7 from the Oberwolfach Benchmark Collection [26]. This is a twodimensional computational fluid dynamics application of thermal flow modelled by a convection-diffusion partial differential equation: with the solution temperature T (x, t), the thermal conductivity κ, the fluid speed v of fixed direction, and the heat generation rateq. The model is discretized in space using the finite element method, yielding a generalized linear system (1) of order N = 9669 ≈ 10 4 , a single input M = 1 and five outputs Q = 5. For a more detailed description of this benchmark see [26] and references therein. This model is tested in two variants: First, in a diffusion-dominant setting with zero flow speed v = 0, and second, in a convection-dominant setting with a flow speed v = 0.5. Due to the MIMO nature of the system we use the average system (see (3)) for the error indicator computation. This set of experiments is organized in the same manner as Section 4.2, but conducted for the prescribed projection errors ε ∈ {10 −2 , . . . , 10 −8 }. As indicated in Section 3.4, the average system (averaged over outputs)C := Q q=1 c q is used for the computation of the error indicator. The resulting reduced order models are tested with impulse input u(t) = δ(t).

Diffusion-Dominant Variant
The experimental results of the diffusion-dominant variant (v = 0) is depicted in Fig. 3. In Fig. 3a the prescribed projection error for the (empirical) system Gramians versus the resulting L 2 model reduction error is plotted. As in Section 4.2, the DSPMR method produces the reduced order models with the lowest model reduction error. Reduced systems from (empirical) balanced truncation and the (empirical) cross-Gramian-based dominant subspace method result in similar errors, while the error indicator behaves behaves like a conservative upper bound to the cross-Gramian-based model reduction error. Fig. 3b shows the model reduction error for the reduced orders resulting from the prescribed projection error. Balanced truncation achieves the smallest and DSPMR the largest reduced models, the cross-Gramian-based dominant subspace method reduced order model dimension lies in between, and the error indicator shows a similar behavior as the latter.

Convection-Dominant Variant
The experimental results of the convection-dominant variant (v = 0.5) is presented in Fig. 4. Overall, the plots Fig. 4a and Fig. 4b are similar to the diffusion-dominant variant, yet the model reduction errors are higher and flatten out for the cross-Gramian-based method for ε < 10 −4 . This means that the model reduction error does not follow the error indicator, which is equal for diffusion-and convection-dominant benchmark variants.