Serkan Gugercin

Zlatko Drmac, Serkan Gugercin

This paper introduces a new framework for constructing the Discrete Empirical Interpolation Method DEIM projection operator. The interpolation node selection procedure is formulated using the QR factorization with column pivoting, and it enjoys a sharper error bound for the DEIM projection error. Furthermore, for a subspace $\mathcal{U}$ given as the range of an orthonormal $U$, the DEIM projection does not change if $U$ is replaced by $U Ω$ with arbitrary unitary matrix $Ω$. In a large-scale setting, the new approach allows modifications that use only randomly sampled rows of $U$, but with the potential of producing good approximations with corresponding probabilistic error bounds. Another salient feature of the new framework is that robust and efficient software implementation is easily developed, based on readily available high performance linear algebra packages.

Sean Reiter, Igor Pontes Duff, Ion Victor Gosea et al.

In this work, we consider the $H_2$ optimal model reduction of dynamical systems that are linear in the state equation and up to quadratic nonlinearity in the output equation. As our primary theoretical contributions, we derive gradients of the squared $H_2$ system error with respect to the reduced model quantities and, from the stationary points of these gradients, introduce Gramian-based first-order necessary conditions for the $H_2$ optimal approximation of a linear quadratic output (LQO) system. The resulting $H_2$ optimality framework neatly generalizes the analogous Gramian-based optimality framework for purely linear systems. Computationally, we show how to enforce the necessary optimality conditions using Petrov-Galerkin projection; the corresponding projection matrices are obtained from a pair of Sylvester equations. Based on this result, we propose an iteratively corrected algorithm for the $H_2$ model reduction of LQO systems, which we refer to as LQO-TSIA (linear quadratic output two-sided iteration algorithm). Numerical examples are included to illustrate the effectiveness of the proposed computational method against other existing approaches.

Peter Benner, Pawan Goyal, Serkan Gugercin

We investigate the optimal model reduction problem for large-scale quadratic-bilinear (QB) control systems. Our contributions are threefold. First, we discuss the variational analysis and the Volterra series formulation for QB systems. We then define the $\mathcal H_2$-norm for a QB system based on the kernels of the underlying Volterra series and also propose a truncated $\mathcal H_2$-norm. Next, we derive first-order necessary conditions for an optimal approximation, where optimality is measured in term of the truncated $\mathcal H_2$-norm of the error system. We then propose an iterative model reduction algorithm, which upon convergence yields a reduced-order system that approximately satisfies the newly derived optimality conditions. We also discuss an efficient computation of the reduced Hessian, using the special Kronecker structure of the Hessian of the system. We illustrate the efficiency of the proposed method by means of several numerical examples resulting from semi-discretized nonlinear partial differential equations and show its competitiveness with the existing model reduction schemes for QB systems such as moment-matching methods and balanced truncation.

Serkan Gugercin, Tatjana Stykel, Sarah Wyatt

In this paper, we investigate interpolatory projection framework for model reduction of descriptor systems. With a simple numerical example, we first illustrate that employing subspace conditions from the standard state space settings to descriptor systems generically leads to unbounded H2 or H-infinity errors due to the mismatch of the polynomial parts of the full and reduced-order transfer functions. We then develop modified interpolatory subspace conditions based on the deflating subspaces that guarantee a bounded error. For the special cases of index-1 and index-2 descriptor systems, we also show how to avoid computing these deflating subspaces explicitly while still enforcing interpolation. The question of how to choose interpolation points optimally naturally arises as in the standard state space setting. We answer this question in the framework of the H2-norm by extending the Iterative Rational Krylov Algorithm (IRKA) to descriptor systems. Several numerical examples are used to illustrate the theoretical discussion.

Christopher A. Beattie, Serkan Gugercin, Volker Mehrmann

We consider the model reduction problem for linear time-invariant dynamical systems having nonzero (but otherwise indeterminate) initial conditions. Building upon the observation that the full system response is decomposable as a superposition of the response map for an unforced system having nontrivial initial conditions and the response map for a forced system having null initial conditions, we develop a new approach that involves reducing these component responses independently and then combining the reduced responses into an aggregate reduced system response. This approach allows greater flexibility and offers better approximation properties than other comparable methods.

Sean Reiter, Mark Embree, Serkan Gugercin et al.

This work investigates the reduction of phasor measurement unit (PMU) data through low-rank matrix approximations. To reconstruct a PMU data matrix from fewer measurements, we propose the framework of interpolatory matrix decompositions (IDs). In contrast to methods relying on principal component analysis or singular value decomposition, IDs recover the complete data matrix using only a few of its rows (PMU datastreams) and/or a few of its columns (snapshots in time). This row-/column-based compression enables real-time monitoring of power transmission systems using measurements from a smaller subset of pilot datastreams, thereby minimizing communication bandwidth. The ID perspective gives a rigorous error bound on the quality of the data compression. We propose selecting the pilot measurements used in an ID via the discrete empirical interpolation method (DEIM), a greedy algorithm that aims to control the error bound. This bound yields a computable estimate of the reconstruction error during online operations. A violation of this estimate suggests a change in the system's operating conditions and thus serves as a tool for fault detection. Following a disturbance, DEIM can be used to localize the event source across all buses with high accuracy. Numerical tests on synthetic PMU data demonstrate DEIM's excellent performance in data compression and validate the proposed DEIM-based fault-detection and localization method.

Sean Reiter, Ion Victor Gosea, Igor Pontes Duff et al.

This paper addresses the $\mathcal{H}_2$-optimal approximation of linear dynamical systems with quadratic-output functions, also known as linear quadratic-output systems. Our major contributions are threefold. First, we derive interpolatory first-order optimality conditions for the linear quadratic-output $\mathcal{H}_2$ minimization problem. These conditions correspond to the mixed-multipoint tangential interpolation of the full-order linear- and quadratic-output transfer functions, and generalize the Meier-Luenberger optimality framework for the $\mathcal{H}_2$-optimal model reduction of linear time-invariant systems. Second, given the optimal interpolation data, we show how to enforce the interpolatory optimality conditions explicitly by Petrov-Galerkin projection of the full-order model. Third, to find the optimal interpolation data, we build on this projection framework and propose a generalization of the iterative rational Krylov algorithm for the $\mathcal{H}_2$-optimal model reduction of linear quadratic-output systems, called LQO-IRKA. Upon convergence, LQO-IRKA produces reduced linear quadratic-output systems that satisfy the interpolatory optimality conditions. The method only requires solving shifted linear systems and matrix-vector products, thus making it suitable for large-scale problems. Numerical examples are included to illustrate the effectiveness of the proposed method.

Christopher Beattie, Serkan Gugercin

The last two decades have seen major developments in interpolatory methods for model reduction of large-scale linear dynamical systems. Advances of note include the ability to produce (locally) optimal reduced models at modest cost; refined methods for deriving interpolatory reduced models directly from input/output measurements; and extensions for the reduction of parametrized systems. This chapter offers a survey of interpolatory model reduction methods starting from basic principles and ranging up through recent developments that include weighted model reduction and structure-preserving methods based on generalized coprime representations. Our discussion is supported by an assortment of numerical examples.

Garret Flagg, Christopher Beattie, Serkan Gugercin

The Iterative Rational Krylov Algorithm (IRKA) of [8] is an interpolatory model reduction approach to the optimal $\mathcal{H}_2$ approximation problem. Even though the method has been illustrated to show rapid convergence in various examples, a proof of convergence has not been provided yet. In this note, we show that in the case of state-space symmetric systems, IRKA is a locally convergent fixed point iteration to a local minimum of the underlying $\mathcal{H}_2$ approximation problem.

Philipp Schulze, Benjamin Unger, Christopher Beattie et al.

We present a framework for constructing structured realizations of linear dynamical systems having transfer functions of the form $C(\sum_{k=1}^K h_k(s)A_k)^{-1}B$ where $h_1,h_2,\ldots,h_K$ are prescribed functions that specify the surmised structure of the model. Our construction is data-driven in the sense that an interpolant is derived entirely from measurements of a transfer function. Our approach extends the Loewner realization framework to more general system structure that includes second-order (and higher) systems as well as systems with internal delays. Numerical examples demonstrate the advantages of this approach.

Kapil Ahuja, Eric de Sturler, Serkan Gugercin et al.

Science and engineering problems frequently require solving a sequence of dual linear systems. Besides having to store only few Lanczos vectors, using the BiConjugate Gradient method (BiCG) to solve dual linear systems has advantages for specific applications. For example, using BiCG to solve the dual linear systems arising in interpolatory model reduction provides a backward error formulation in the model reduction framework. Using BiCG to evaluate bilinear forms -- for example, in quantum Monte Carlo (QMC) methods for electronic structure calculations -- leads to a quadratic error bound. Since our focus is on sequences of dual linear systems, we introduce recycling BiCG, a BiCG method that recycles two Krylov subspaces from one pair of dual linear systems to the next pair. The derivation of recycling BiCG also builds the foundation for developing recycling variants of other bi-Lanczos based methods, such as CGS, BiCGSTAB, QMR, and TFQMR. We develop an augmented bi-Lanczos algorithm and a modified two-term recurrence to include recycling in the iteration. The recycle spaces are approximate left and right invariant subspaces corresponding to the eigenvalues closest to the origin. These recycle spaces are found by solving a small generalized eigenvalue problem alongside the dual linear systems being solved in the sequence. We test our algorithm in two application areas. First, we solve a discretized partial differential equation (PDE) of convection-diffusion type. Such a problem provides well-known test cases that are easy to test and analyze further. Second, we use recycling BiCG in the Iterative Rational Krylov Algorithm (IRKA) for interpolatory model reduction. IRKA requires solving sequences of slowly changing dual linear systems. We show up to 70% savings in iterations, and also demonstrate that for a model reduction problem BiCG takes (about) 50% more time than recycling BiCG.

Zlatko Drmac, Serkan Gugercin, Christopher Beattie

Vector Fitting (VF) is a popular method of constructing rational approximants that provides a least squares fit to frequency response measurements. In an earlier work, we provided an analysis of VF for scalar-valued rational functions and established a connection with optimal $H_2$ approximation. We build on this work and extend the previous framework to include the construction of effective rational approximations to matrix-valued functions, a problem which presents significant challenges that do not appear in the scalar case. Transfer functions associated with multi-input/multi-output (MIMO) dynamical systems typify the class of functions that we consider here. Others have also considered extensions of VF to matrix-valued functions and related numerical implementations are readily available. However to our knowledge, a detailed analysis of numerical issues that arise does not yet exist. We offer such an analysis including critical implementation details here. One important issue that arises for VF on matrix-valued functions that has remained largely unaddressed is the control of the McMillan degree of the resulting rational approximant; the McMillan degree can grow very high in the case of large input/output dimensions. We introduce two new mechanisms for controlling the McMillan degree of the final approximant, one based on alternating least-squares minimization and one based on ancillary system-theoretic reduction methods. Motivated in part by our earlier work on the scalar VF problem as well as by recent innovations for computing optimal $H_2$ approximation, we establish a connection with optimal $H_2$ approximation, and are able to improve significantly the fidelity of VF through numerical quadrature, with virtually no increase in cost or complexity. We provide several numerical examples to support the theoretical discussion and proposed algorithms.

Serkan Gugercin, Rostyslav V. Polyuga, Christopher Beattie et al.

Port-Hamiltonian systems result from port-based network modeling of physical systems and are an important example of passive state-space systems. In this paper, we develop the framework for model reduction of large-scale multi-input/multi-output port-Hamiltonian systems via tangential rational interpolation. The resulting reduced-order model not only is a rational tangential interpolant but also retains the port-Hamiltonian structure; hence is passive. This reduction methodology is described in both energy and co-energy system coordinates. We also introduce an $\mathcal{H}_2$-inspired algorithm for effectively choosing the interpolation points and tangential directions. The algorithm leads a reduced port-Hamiltonian model that satisfies a subset of $\mathcal{H}_2$-optimality conditions. We present several numerical examples that illustrate the effectiveness of the proposed method showing that it outperforms other existing techniques in both quality and numerical efficiency.

Benjamin Unger, Serkan Gugercin

Kolmogorov $n$-widths and Hankel singular values are two commonly used concepts in model reduction. Here we show that for the special case of linear time-invariant dynamical (LTI) systems, these two concepts are directly connected. More specifically, the greedy search applied to the Hankel operator of an LTI system resembles the minimizing subspace for the Kolmogorov n-width and the Kolmogorov $n$-width of an LTI system equals its $(n+1)st$ Hankel singular value once the subspaces are appropriately defined. We also establish a lower bound for the Kolmorogov $n$-width for parametric LTI systems and illustrate that the method of active subspaces can be viewed as the dual concept to the minimizing subspace for the Kolmogorov $n$-width.

Tobias Breiten, Christopher Beattie, Serkan Gugercin

This paper develops an interpolatory framework for weighted-$\mathcal{H}_2$ model reduction of MIMO dynamical systems. A new representation of the weighted-$\mathcal{H}_2$ inner products in MIMO settings is introduced and used to derive associated first-order necessary conditions satisfied by optimal weighted-$\mathcal{H}_2$ reduced-order models. Equivalence of these new interpolatory conditions with earlier Riccati-based conditions given by Halevi is also shown. An examination of realizations for equivalent weighted-$\mathcal{H}_2$ systems leads then to an algorithm that remains tractable for large state-space dimension. Several numerical examples illustrate the effectiveness of this approach and its competitiveness with Frequency Weighted Balanced Truncation and an earlier interpolatory approach, the Weighted Iterative Rational Krylov Algorithm.

Garret Flagg, Christopher Beattie, Serkan Gugercin

We introduce an interpolation framework for H-infinity model reduction founded on ideas originating in optimal-H2 interpolatory model reduction, realization theory, and complex Chebyshev approximation. By employing a Loewner "data-driven" framework within each optimization cycle, large-scale H-infinity norm calculations can be completely avoided. Thus, we are able to formulate a method that remains effective in large-scale settings with the main cost dominated by sparse linear solves. Several numerical examples illustrate that our approach will produce high fidelity reduced models consistently exhibiting better H-infinity performance than those produced by balanced truncation; these models often are as good as (and occasionally better than) those models produced by optimal Hankel norm approximation. In all cases, these reduced models are produced at far lower cost than is possible either with balanced truncation or optimal Hankel norm approximation.

Christopher A. Beattie, Serkan Gugercin, Sarah Wyatt

We investigate the use of inexact solves for interpolatory model reduction and consider associated perturbation effects on the underlying model reduction problem. We give bounds on system perturbations induced by inexact solves and relate this to termination criteria for iterative solution methods. We show that when a Petrov-Galerkin framework is employed for the inexact solves, the associated reduced order model is an exact interpolatory model for a nearby full-order system; thus demonstrating backward stability. We also give evidence that for $\h2$-optimal interpolation points, interpolatory model reduction is robust with respect to perturbations due to inexact solves. Finally, we demonstrate the effecitveness of direct use of inexact solves in optimal ${\mathcal H}_2$ approximation. The result is an effective model reduction strategy that is applicable in realistically large-scale settings.

Branimir Anic, Christopher A. Beattie, Serkan Gugercin et al.

This paper introduces an interpolation framework for the weighted-H2 model reduction problem. We obtain a new representation of the weighted-H2 norm of SISO systems that provides new interpolatory first order necessary conditions for an optimal reduced-order model. The H2 norm representation also provides an error expression that motivates a new weighted-H2 model reduction algorithm. Several numerical examples illustrate the effectiveness of the proposed approach.

Saifon Chaturantabut, Chris Beattie, Serkan Gugercin

This paper presents a structure-preserving model reduction approach applicable to large-scale, nonlinear port-Hamiltonian systems. Structure preservation in the reduction step ensures the retention of port-Hamiltonian structure which, in turn, assures the stability and passivity of the reduced model. Our analysis provides a priori error bounds for both state variables and outputs. Three techniques are considered for constructing bases needed for the reduction: one that utilizes proper orthogonal decompositions; one that utilizes $\mathcal{H}_2/\mathcal{H}_{\infty}$-derived optimized bases; and one that is a mixture of the two. The complexity of evaluating the reduced nonlinear term is managed efficiently using a modification of the discrete empirical interpolation method (DEIM) that also preserves port-Hamiltonian structure. The efficiency and accuracy of this model reduction framework are illustrated with two examples: a nonlinear ladder network and a tethered Toda lattice.

Garret M. Flagg, Serkan Gugercin

The ADI iteration is closely related to the rational Krylov projection methods for constructing low rank approximations to the solution of Sylvester equation. In this paper we show that the ADI and rational Krylov approximations are in fact equivalent when a special choice of shifts are employed in both methods. We will call these shifts pseudo H2-optimal shifts. These shifts are also optimal in the sense that for the Lyapunov equation, they yield a residual which is orthogonal to the rational Krylov projection subspace. Via several examples, we show that the pseudo H2-optimal shifts consistently yield nearly optimal low rank approximations to the solutions of the Lyapunov equations.

Andrea Carracedo Rodriguez, Serkan Gugercin, Jeff Borggaard

Interpolatory projection methods for model reduction of nonparametric linear dynamical systems have been successfully extended to nonparametric bilinear dynamical systems. However, this is not the case for parametric bilinear systems. In this work, we aim to close this gap by providing a natural extension of interpolatory projections to model reduction of parametric bilinear dynamical systems. We introduce necessary conditions that the projection subspaces must satisfy to obtain parametric tangential interpolation of each subsystem transfer function. These conditions also guarantee that the parameter sensitivities (Jacobian) of each subsystem transfer function is matched tangentially by those of the corresponding reduced order model transfer function. Similarly, we obtain conditions for interpolating the parameter Hessian of the transfer function by including extra vectors in the projection subspaces. As in the parametric linear case, the basis construction for two-sided projections does not require computing the Jacobian or the Hessian.

Meghan O'Connell, Misha E. Kilmer, Eric de Sturler et al.

In nonlinear imaging problems whose forward model is described by a partial differential equation (PDE), the main computational bottleneck in solving the inverse problem is the need to solve many large-scale discretized PDEs at each step of the optimization process. In the context of absorption imaging in diffuse optical tomography, one approach to addressing this bottleneck proposed recently (de Sturler, et al, 2015) reformulates the viewing of the forward problem as a differential algebraic system, and then employs model order reduction (MOR). However, the construction of the reduced model requires the solution of several full order problems (i.e. the full discretized PDE for multiple right-hand sides) to generate a candidate global basis. This step is then followed by a rank-revealing factorization of the matrix containing the candidate basis in order to compress the basis to a size suitable for constructing the reduced transfer function. The present paper addresses the costs associated with the global basis approximation in two ways. First, we use the structure of the matrix to rewrite the full order transfer function, and corresponding derivatives, such that the full order systems to be solved are symmetric (positive definite in the zero frequency case). Then we apply MOR to the new formulation of the problem. Second, we give an approach to computing the global basis approximation dynamically as the full order systems are solved. In this phase, only the incrementally new, relevant information is added to the existing global basis, and redundant information is not computed. This new approach is achieved by an inner-outer Krylov recycling approach which has potential use in other applications as well. We show the value of the new approach to approximate global basis computation on two DOT absorption image reconstruction problems.

Cankat Tilki, Tobias Breiten, Serkan Gugercin

We develop the interpolatory $\mathcal{H}_2$ optimal model reduction framework for linear control systems posed on infinite dimensional state, input and output spaces. Specifically, we consider linear systems formulated as controlled abstract Cauchy problems on a Banach space and approximate them via Petrov-Galerkin projection onto finite dimensional trial and test subspaces. We show that the resulting reduced order transfer function interpolates the original at prescribed points, and we characterize precisely how the projection subspaces must be constructed to enforce this interpolation. Building on this, we develop a data-driven realization framework -- an infinite dimensional analogue of the Loewner approach -- that recovers the system behavior directly from input-output data without requiring access to the underlying operators. Finally, we derive $\mathcal{H}_2$ optimality conditions for the reduced model and show that the classical interpolatory characterization persists in this infinite dimensional setting: first-order optimality requires Hermite interpolation of the transfer function at the mirror images of the reduced model's poles. Taken together, these results establish that the interpolatory $\mathcal{H}_2$ optimal model reduction theory extends naturally and completely to infinite dimensional linear control systems with infinite dimensional input and output spaces.

Klajdi Sinani, Serkan Gugercin

In this paper we establish the interpolatory model reduction framework for optimal approximation of MIMO dynamical systems with respect to the $\mathcal{H}_2$ norm over a finite-time horizon, denoted as the $\mathcal{H}_2(t_f)$ norm. Using the underlying inner product space, we derive the interpolatory first-order necessary optimality conditions for approximation in the $\mathcal{H}_2(t_f)$ norm. Then, we develop an algorithm, which yields a locally optimal reduced model that satisfies the established interpolation-based optimality conditions. We test the algorithm on various numerical examples to illustrate its performance.

Alexander Grimm, Christopher Beattie, Zlatko Drmač et al.

Given a set of response observations for a parametrized dynamical system, we seek a parametrized dynamical model that will yield uniformly small response error over a range of parameter values yet has low order. Frequently, access to internal system dynamics or equivalently, to realizations of the original system is either not possible or not practical; only response observations over a range of parameter settings might be known. Respecting these typical operational constraints, we propose a two phase approach that first encodes the response data into a high fidelity intermediate model of modest order, followed then by a compression stage that serves to eliminate redundancy in the intermediate model. For the first phase, we extend non-parametric least-squares fitting approaches so as to accommodate parameterized systems. This results in a (discrete) least-squares problem formulated with respect to both frequency and parameter that identifies "local" system response features. The second phase uses an $\mathbf{\mathcal{H}}_2$-optimal model reduction strategy accommodating the specialized parametric structure of the intermediate model obtained in the first phase. The final compressed model inherits the parametric dependence of the intermediate model and maintains the high fidelity of the intermediate model, while generally having dramatically smaller system order. We provide a variety of numerical examples demonstrating our approach.

Cankat Tilki, Serkan Gugercin

We present an in-depth analysis of the Koopman semigroup via wavelet transform. Towards this goal, we start by introducing the wavelet-based observables and show that they are eigenfunctions of the Koopman semigroup when this semigroup is considered over the Banach space of continuous functions on a compact forward-invariant set endowed with the supremum norm. We then construct closed-form expressions of the action of the Koopman semigroup and its resolvent in terms of these observables. To approximate the action of Koopman semigroup numerically, we combine Extended Dynamic Mode Decomposition (EDMD) with the proposed wavelet-based observables leading to the Wavelet Dynamic Mode Decomposition via Continuous Wavelet Transform (cWDMD) algorithm. We validate our theoretical results on two numerical examples.

Michael S. Ackermann, Ion Victor Gosea, Serkan Gugercin et al.

The data-driven modeling of dynamical systems has become an essential tool for the construction of accurate computational models from real-world data. In this process, the inherent differential structures underlying the considered physical phenomena are often neglected making the reinterpretation of the learned models in a physically meaningful sense very challenging. In this work, we present three data-driven modeling approaches for the construction of dynamical systems with second-order differential structure directly from frequency domain data. Based on the second-order structured barycentric form, we extend the well-known Adaptive Antoulas-Anderson algorithm to the case of second-order systems. Depending on the available computational resources, we propose variations of the proposed method that prioritize either higher computation speed or greater modeling accuracy, and we present a theoretical analysis for the expected accuracy and performance of the proposed methods. Three numerical examples demonstrate the effectiveness of our new structured approaches in comparison to classical unstructured data-driven modeling.

Caleb C. Magruder, Serkan Gugercin, Christopher A. Beattie

Linear time-periodic (LTP) dynamical systems frequently appear in the modeling of phenomena related to fluid dynamics, electronic circuits, and structural mechanics via linearization centered around known periodic orbits of nonlinear models. Such LTP systems can reach orders that make repeated simulation or other necessary analysis prohibitive, motivating the need for model reduction. We develop here an algorithmic framework for constructing reduced models that retains the linear time-periodic structure of the original LTP system. Our approach generalizes optimal approaches that have been established previously for linear time-invariant (LTI) model reduction problems. We employ an extension of the usual H2 Hardy space defined for the LTI setting to time-periodic systems and within this broader framework develop an a posteriori error bound expressible in terms of related LTI systems. Optimization of this bound motivates our algorithm. We illustrate the success of our method on two numerical examples.