Zlatko Drmač

Zlatko Drmač, Igor Mezić, Ryan Mohr

The goals and contributions of this paper are twofold. It provides a new computational tool for data driven Koopman spectral analysis by taking up the formidable challenge to develop a numerically robust algorithm by following the natural formulation via the Krylov decomposition with the Frobenius companion matrix, and by using its eigenvectors explicitly -- these are defined as the inverse of the notoriously ill-conditioned Vandermonde matrix. The key step to curb ill-conditioning is the discrete Fourier transform of the snapshots; in the new representation, the Vandermonde matrix is transformed into a generalized Cauchy matrix, which then allows accurate computation by specially tailored algorithms of numerical linear algebra. The second goal is to shed light on the connection between the formulas for optimal reconstruction weights when reconstructing snapshots using subsets of the computed Koopman modes. It is shown how using a certain weaker form of generalized inverses leads to explicit reconstruction formulas that match the abstract results from Koopman spectral theory, in particular the Generalized Laplace Analysis.

Zlatko Drmač, Arvind K. Saibaba

New contributions are offered to the theory and practice of the Discrete Empirical Interpolation Method (DEIM). These include a detailed characterization of the canonical structure; a substantial tightening of the error bound for the DEIM oblique projection, based on index selection via a strong rank revealing QR factorization; and an extension of the DEIM approximation to weighted inner products defined by a real symmetric positive-definite matrix $W$. The weighted DEIM ($W$-DEIM) can be deployed in the more general framework where the POD Galerkin projection is formulated in a discretization of a suitable energy inner product such that the Galerkin projection preserves important physical properties such as e.g. stability. Also, a special case of $W$-DEIM is introduced, which is DGEIM, a discrete version of the Generalized Empirical Interpolation Method that allows generalization of the interpolation via a dictionary of linear functionals.

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.

Zlatko Drmač, Ela Đimoti

The Dynamic Mode Decomposition (DMD) and the more general Extended DMD (EDMD) are powerful tools for computational analysis of dynamical systems in data-driven scenarios. They are built on the theoretical foundation of the Koopman composition operator and can be considered as numerical methods for data snapshot-based extraction of spectral information of the composition operator associated with the dynamics, spectral analysis of the structure of the dynamics, and for forecasting. In high fidelity numerical simulations, the state space is high dimensional and efficient numerical methods leverage the fact that the actual dynamics evolves on manifolds of much smaller dimension. This motivates computing low rank approximations in a streaming fashion and the DMD matrix is adaptively updated with newly received data. In this way, large number of high dimensional snapshots can be processed very efficiently. Low dimensional representation also requires fast updating for online applications. This paper revisits the pioneering works of Hemati, Williams and Rowley (Physics of Fluids, 2014), and Zhang, Rowley, Deem and Cattafesta (SIAM Journal on Applied Dynamical Systems, 2019) on the streaming DMD and proposes improvements in functionality (using residual bounds, Exact DMD vectors), computational efficiency (more efficient algorithm with smaller memory footprint) and numerical robustness (smaller condition numbers and better forecasting skill).

Ryan Mohr, Maria Fonoberova, Zlatko Drmač et al.

Hierarchical support vector regression (HSVR) models a function from data as a linear combination of SVR models at a range of scales, starting at a coarse scale and moving to finer scales as the hierarchy continues. In the original formulation of HSVR, there were no rules for choosing the depth of the model. In this paper, we observe in a number of models a phase transition in the training error -- the error remains relatively constant as layers are added, until a critical scale is passed, at which point the training error drops close to zero and remains nearly constant for added layers. We introduce a method to predict this critical scale a priori with the prediction based on the support of either a Fourier transform of the data or the Dynamic Mode Decomposition (DMD) spectrum. This allows us to determine the required number of layers prior to training any models.

Iva Manojlović, Maria Fonoberova, Ryan Mohr et al.

We consider the training process of a neural network as a dynamical system acting on the high-dimensional weight space. Each epoch is an application of the map induced by the optimization algorithm and the loss function. Using this induced map, we can apply observables on the weight space and measure their evolution. The evolution of the observables are given by the Koopman operator associated with the induced dynamical system. We use the spectrum and modes of the Koopman operator to realize the above objectives. Our methods can help to, a priori, determine the network depth; determine if we have a bad initialization of the network weights, allowing a restart before training too long; speeding up the training time. Additionally, our methods help enable noise rejection and improve robustness. We show how the Koopman spectrum can be used to determine the number of layers required for the architecture. Additionally, we show how we can elucidate the convergence versus non-convergence of the training process by monitoring the spectrum, in particular, how the existence of eigenvalues clustering around 1 determines when to terminate the learning process. We also show how using Koopman modes we can selectively prune the network to speed up the training procedure. Finally, we show that incorporating loss functions based on negative Sobolev norms can allow for the reconstruction of a multi-scale signal polluted by very large amounts of noise.

Zlatko Drmač, Benjamin Peherstorfer

Loewner rational interpolation provides a versatile tool to learn low-dimensional dynamical-system models from frequency-response measurements. This work investigates the robustness of the Loewner approach to noise. The key finding is that if the measurements are polluted with Gaussian noise, then the error due to noise grows at most linearly with the standard deviation with high probability under certain conditions. The analysis gives insights into making the Loewner approach robust against noise via linear transformations and judicious selections of measurements. Numerical results demonstrate the linear growth of the error on benchmark examples.

Zlatko Drmač, Ivana Šain Glibić

This paper presents a new method for computing all eigenvalues and eigenvectors of quadratic matrix pencil. It is an upgrade of the quadeig algorithm by Hammarling, Munro and Tisseur, which attempts to reveal and remove by deflation certain number of zero and infinite eigenvalues before QZ iterations. Proposed modifications of the quadeig framework are designed to enhance backward stability and to make the process of deflating infinite and zero eigenvalues more numerically robust. In particular, careful preprocessing allows scaling invariant/component-wise backward error and thus better condition number. Further, using an upper triangular version of the Kronecker canonical form enables deflating additional infinite eigenvalues, in addition to those inferred from the rank of leading coefficient matrix. Theoretical analysis and empirical evidence from thorough testing of the software implementation confirm superior numerical performances of the proposed method.

Nela Bosner, Zvonimir Bujanović, Zlatko Drmač

This paper proposes a combination of a hybrid CPU--GPU and a pure GPU software implementation of a direct algorithm for solving shifted linear systems $(A - σI)X = B$ with large number of complex shifts $σ$ and multiple right-hand sides. Such problems often appear e.g. in control theory when evaluating the transfer function, or as a part of an algorithm performing interpolatory model reduction, as well as when computing pseudospectra and structured pseudospectra, or solving large linear systems of ordinary differential equations. The proposed algorithm first jointly reduces the general full $n\times n$ matrix $A$ and the $n\times m$ full right-hand side matrix $B$ to the controller Hessenberg canonical form that facilitates efficient solution: $A$ is transformed to a so-called $m$-Hessenberg form and $B$ is made upper-triangular. This is implemented as blocked highly parallel CPU--GPU hybrid algorithm; individual blocks are reduced by the CPU, and the necessary updates of the rest of the matrix are split among the cores of the CPU and the GPU. To enhance parallelization, the reduction and the updates are overlapped. In the next phase, the reduced $m$-Hessenberg--triangular systems are solved entirely on the GPU, with shifts divided into batches. The benefits of such load distribution are demonstrated by numerical experiments. In particular, we show that our proposed implementation provides an excellent basis for efficient implementations of computational methods in systems and control theory, from evaluation of transfer function to the interpolatory model reduction.

Zlatko Drmač, Igor Mezić, Ryan Mohr

The Dynamic Mode Decomposition (DMD) is a tool of trade in computational data driven analysis of fluid flows. More generally, it is a computational device for Koopman spectral analysis of nonlinear dynamical systems, with a plethora of applications in applied sciences and engineering. Its exceptional performance triggered developments of several modifications that make the DMD an attractive method in data driven framework. This work offers further improvements of the DMD to make it more reliable, and to enhance its functionality. In particular, data driven formula for the residuals allows selection of the Ritz pairs, thus providing more precise spectral information of the underlying Koopman operator, and the well-known technique of refining the Ritz vectors is adapted to data driven scenarios. Further, the DMD is formulated in a more general setting of weighted inner product spaces, and the consequences for numerical computation are discussed in detail. Numerical experiments are used to illustrate the advantages of the proposed method, designated as DDMD_RRR (Refined Rayleigh Ritz Data Driven Modal Decomposition).