Peter Benner

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.

Maximilian Behr, Peter Benner, Jan Heiland

The differential Sylvester equation and its symmetric version, the differential Lyapunov equation, appear in different fields of applied mathematics like control theory, system theory, and model order reduction. The few available straight-forward numerical approaches if applied to large-scale systems come with prohibitively large storage requirements. This shortage motivates us to summarize and explore existing solution formulas for these equations. We develop a unifying approach based on the spectral theorem for normal operators like the Sylvester operator $\mathcal S(X)=AX+XB$ and derive a formula for its norm using an induced operator norm based on the spectrum of $A$ and $B$. In view of numerical approximations, we propose an algorithm that identifies a suitable Krylov subspace using Taylor series and use a projection to approximate the solution. Numerical results for large-scale differential Lyapunov equations are presented in the last sections.

Peter Benner, Zvonimir Bujanović, Patrick Kürschner et al.

This paper introduces a new algorithm for solving large-scale continuous-time algebraic Riccati equations (CARE). The advantage of the new algorithm is in its immediate and efficient low-rank formulation, which is a generalization of the Cholesky-factored variant of the Lyapunov ADI method. We discuss important implementation aspects of the algorithm, such as reducing the use of complex arithmetic and shift selection strategies. We show that there is a very tight relation between the new algorithm and three other algorithms for CARE previously known in the literature -- all of these seemingly different methods in fact produce exactly the same iterates when used with the same parameters: they are algorithmically different descriptions of the same approximation sequence to the Riccati solution.

Peter Benner, Christian Himpe, Tim Mitchell

The identification of reduced-order models from high-dimensional data is a challenging task, and even more so if the identified system should not only be suitable for a certain data set, but generally approximate the input-output behavior of the data source. In this work, we consider the input-output dynamic mode decomposition method for system identification. We compare excitation approaches for the data-driven identification process and describe an optimization-based stabilization strategy for the identified systems.

Peter Benner, Sergey Dolgov, Venera Khoromskaia et al.

In this paper, we study and implement the structural iterative eigensolvers for the large-scale eigenvalue problem in the Bethe-Salpeter equation (BSE) based on the reduced basis approach via low-rank factorizations in generating matrices, introduced in the previous paper. The approach reduces numerical costs down to $\mathcal{O}(N_b^2)$ in the size of atomic orbitals basis set, $N_b$, instead of practically intractable $\mathcal{O}(N_b^6)$ complexity scaling for the direct diagonalization of the BSE matrix. As an alternative to rank approximation of the static screen interaction part of the BSE matrix, we propose to restrict it to a small active sub-block, with a size balancing the storage for rank-structured representations of other matrix blocks. We demonstrate that the enhanced reduced-block approximation exhibits higher precision within the controlled numerical cost, providing as well a distinct two-sided error estimate for the BSE eigenvalues. It is shown that further reduction of the asymptotic computational cost is possible due to ALS-type iteration in block tensor train (TT) format applied to the quantized-TT (QTT) tensor representation of both long eigenvectors and rank-structured matrix blocks. The QTT-rank of these entities possesses almost the same magnitude as the number of occupied orbitals in the molecular systems, $N_o$, hence the overall asymptotic complexity for solving the BSE problem can be estimated by $\mathcal{O}(\log(N_o) N_o^{2})$. We confirm numerically a considerable decrease in computational time for the presented iterative approach applied to various compact and chain-type molecules, while supporting sufficient accuracy.

Pawan Goyal, Igor Pontes Duff, Peter Benner

Scientific machine learning for inferring dynamical systems combines data-driven modeling, physics-based modeling, and empirical knowledge. It plays an essential role in engineering design and digital twinning. In this work, we primarily focus on an operator inference methodology that builds dynamical models, preferably in low-dimension, with a prior hypothesis on the model structure, often determined by known physics or given by experts. Then, for inference, we aim to learn the operators of a model by setting up an appropriate optimization problem. One of the critical properties of dynamical systems is stability. However, this property is not guaranteed by the inferred models. In this work, we propose inference formulations to learn quadratic models, which are stable by design. Precisely, we discuss the parameterization of quadratic systems that are locally and globally stable. Moreover, for quadratic systems with no stable point yet bounded (e.g., chaotic Lorenz model), we discuss how to parameterize such bounded behaviors in the learning process. Using those parameterizations, we set up inference problems, which are then solved using a gradient-based optimization method. Furthermore, to avoid numerical derivatives and still learn continuous systems, we make use of an integral form of differential equations. We present several numerical examples, illustrating the preservation of stability and discussing its comparison with the existing state-of-the-art approach to infer operators. By means of numerical examples, we also demonstrate how the proposed methods are employed to discover governing equations and energy-preserving models.

Harshit Kapadia, Lihong Feng, Peter Benner

When repeated evaluations for varying parameter configurations of a high-fidelity physical model are required, surrogate modeling techniques based on model order reduction are desired. In absence of the governing equations describing the dynamics, we need to construct the parametric reduced-order surrogate model in a non-intrusive fashion. In this setting, the usual residual-based error estimate for optimal parameter sampling associated with the reduced basis method is not directly available. Our work provides a non-intrusive optimality criterion to efficiently populate the parameter snapshots, thereby, enabling us to effectively construct a parametric surrogate model. We consider separate parameter-specific proper orthogonal decomposition (POD) subspaces and propose an active-learning-driven surrogate model using kernel-based shallow neural networks, abbreviated as ActLearn-POD-KSNN surrogate model. To demonstrate the validity of our proposed ideas, we present numerical experiments using two physical models, namely Burgers' equation and shallow water equations. Both the models have mixed -- convective and diffusive -- effects within their respective parameter domains, with each of them dominating in certain regions. The proposed ActLearn-POD-KSNN surrogate model efficiently predicts the solution at new parameter locations, even for a setting with multiple interacting shock profiles.

Peter Benner, Akwum Onwunta, Martin Stoll

This paper aims at the efficient numerical solution of stochastic eigenvalue problems. Such problems often lead to prohibitively high dimensional systems with tensor product structure when discretized with the stochastic Galerkin method. Here, we exploit this inherent tensor product structure to develop a globalized low-rank inexact Newton method with which we tackle the stochastic eigenproblem. We illustrate the effectiveness of our solver with numerical experiments.

Yue Qiu, Sara Grundel, Martin Stoll et al.

We study the modeling and simulation of gas pipeline networks, with a focus on fast numerical methods for the simulation of transient dynamics. The obtained mathematical model of the underlying network is represented by a nonlinear differential algebraic equation (DAE). With our modeling, we reduce the number of algebraic constraints, which correspond to the $(2,2)$ block in our semi-explicit DAE model, to the order of junction nodes in the network, where a junction node couples at least three pipelines. We can furthermore ensure that the $(1, 1)$ block of all system matrices including the Jacobian is block lower triangular by using a specific ordering of the pipes of the network. We then exploit this structure to propose an efficient preconditioner for the fast simulation of the network. We test our numerical methods on benchmark problems of (well-)known gas networks and the numerical results show the efficiency of our methods.

Pawan Goyal, Peter Benner

The engineering design process often relies on mathematical modeling that can describe the underlying dynamic behavior. In this work, we present a data-driven methodology for modeling the dynamics of nonlinear systems. To simplify this task, we aim to identify a coordinate transformation that allows us to represent the dynamics of nonlinear systems using a common, simple model structure. The advantage of a common simple model is that customized design tools developed for it can be applied to study a large variety of nonlinear systems. The simplest common model -- one can think of -- is linear, but linear systems often fall short in accurately capturing the complex dynamics of nonlinear systems. In this work, we propose using quadratic systems as the common structure, inspired by the lifting principle. According to this principle, smooth nonlinear systems can be expressed as quadratic systems in suitable coordinates without approximation errors. However, finding these coordinates solely from data is challenging. Here, we leverage deep learning to identify such lifted coordinates using only data, enabling a quadratic dynamical system to describe the system's dynamics. Additionally, we discuss the asymptotic stability of these quadratic dynamical systems. We illustrate the approach using data collected from various numerical examples, demonstrating its superior performance with the existing well-known techniques.

Xin Du, Peter Benner

This paper discusses model order reduction of LTI systems over limited frequency intervals within the framework of balanced truncation. Two new \emph{frequency-dependent balanced truncation} methods were developed, one is \emph{SF-type frequency-dependent balanced truncation} to copy with the cases that only a single dominating point of the operating frequency interval is pre-known, the other is \emph{interval-type frequency-dependent balanced truncation} to deal with the cases that both of the upper and lower bound of frequency interval are known \emph{a priori}. SF-type error bound and interval-type error bound are derived for the first time to estimate the desired approximation error over pre-specified frequency interval. We show that the new methods generally lead to good in-band approximation performance, at the same time, provide accurate error bounds under certain conditions. Examples are included for illustration.

Süleyman Yildiz, Pawan Goyal, Thomas Bendokat et al.

We present a framework for learning Hamiltonian systems using data. This work is based on a lifting hypothesis, which posits that nonlinear Hamiltonian systems can be written as nonlinear systems with cubic Hamiltonians. By leveraging this, we obtain quadratic dynamics that are Hamiltonian in a transformed coordinate system. To that end, for given generalized position and momentum data, we propose a methodology to learn quadratic dynamical systems, enforcing the Hamiltonian structure in combination with a weakly-enforced symplectic auto-encoder. The obtained Hamiltonian structure exhibits long-term stability of the system, while the cubic Hamiltonian function provides relatively low model complexity. For low-dimensional data, we determine a higher-dimensional transformed coordinate system, whereas for high-dimensional data, we find a lower-dimensional coordinate system with the desired properties. We demonstrate the proposed methodology by means of both low-dimensional and high-dimensional nonlinear Hamiltonian systems.

Peter Benner, Christian Himpe

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.

Peter Benner, Norman Lang

We consider the application of implicit and linearly implicit (Rosenbrock-type) peer methods to matrix-valued ordinary differential equations. In particular the differential Riccati equation (DRE) is investigated. For the Rosenbrock-type schemes, a reformulation capable of avoiding a number of Jacobian applications is developed that, in the autonomous case, reduces the computational complexity of the algorithms. Dealing with large-scale problems, an efficient implementation based on low-rank symmetric indefinite factorizations is presented. The performance of both peer approaches up to order 4 is compared to existing implicit time integration schemes for matrix-valued differential equations.

Petar Mlinarić, Takayuki Ishizaki, Aranya Chakrabortty et al.

We study nonlinear power systems consisting of generators, generator buses, and non-generator buses. First, looking at a generator and its bus' variables jointly, we introduce a synchronization concept for a pair of such joint generators and buses. We show that this concept is related to graph symmetry. Next, we extend, in two ways, the synchronization from a pair to a partition of all generators in the networks and show that they are related to either graph symmetry or equitable partitions. Finally, we show how an exact reduced model can be obtained by aggregating the generators and associated buses in the network when the original system is synchronized with respect to a partition, provided that the initial condition respects the partition. Additionally, the aggregation-based reduced model is again a power system.

Cleophas Kweyu, Lihong Feng, Matthias Stein et al.

The Poisson-Boltzmann equation (PBE) is a nonlinear elliptic PDE that arises in biomolecular modeling and is a fundamental tool for structural biology. It is used to calculate electrostatic potentials around an ensemble of fixed charges immersed in an ionic solution. Efficient numerical computation of the PBE yields a high number of degrees of freedom in the resultant algebraic system of equations. Coupled with the fact that in most cases the PBE requires to be solved multiple times for a large number of system configurations, this poses great computational challenges to conventional numerical techniques. To accelerate such computations, we here present the reduced basis method (RBM) which greatly reduces this computational complexity by constructing a reduced order model of typically low dimension. In this study, we employ a simple version of the PBE for proof of concept and discretize the linearized PBE (LPBE) with a centered finite difference scheme. The resultant linear system is solved by the aggregation-based algebraic multigrid method at different samples of ionic strength on a three-dimensional Cartesian grid. The discretized LPBE, which we call the high-fidelity full order model (FOM), yields solution as accurate as other LPBE solvers. We then apply the RBM to FOM. The discrete empirical interpolation method (DEIM) is applied to the Dirichlet boundary conditions which are nonaffine with the parameter, to reduce the complexity of the reduced order model (ROM). From the numerical results, we notice that the RBM reduces the model order from $\mathcal{N} = 2\times 10^{6}$ to $N = 6$ at an accuracy of $10^{-9}$ and reduces computational time by a factor of approximately $7,600$. DEIM, on the other hand, is also used in the offline-online phase of solving the ROM for different values of parameters which provides a speed-up of $20$ for a single iteration of the greedy algorithm.

Peter Benner, Sergey Dolgov, Akwum Onwunta et al.

Many problems in computational science and engineering are simultaneously characterized by the following challenging issues: uncertainty, nonlinearity, nonstationarity and high dimensionality. Existing numerical techniques for such models would typically require considerable computational and storage resources. This is the case, for instance, for an optimization problem governed by time-dependent Navier-Stokes equations with uncertain inputs. In particular, the stochastic Galerkin finite element method often leads to a prohibitively high dimensional saddle-point system with tensor product structure. In this paper, we approximate the solution by the low-rank Tensor Train decomposition, and present a numerically efficient algorithm to solve the optimality equations directly in the low-rank representation. We show that the solution of the vorticity minimization problem with a distributed control admits a representation with ranks that depend modestly on model and discretization parameters even for high Reynolds numbers. For lower Reynolds numbers this is also the case for a boundary control. This opens the way for a reduced-order modeling of the stochastic optimal flow control with a moderate cost at all stages.

Peter Benner, Yue Qiu, Martin Stoll

We consider the problem of estimating the uncertainty in statistical inverse problems using Bayesian inference. When the probability density of the noise and the prior are Gaussian, the solution of such a statistical inverse problem is also Gaussian. Therefore, the underlying solution is characterized by the mean and covariance matrix of the posterior probability density. However, the covariance matrix of the posterior probability density is full and large. Hence, the computation of such a matrix is impossible for large dimensional parameter spaces. It is shown that for many ill-posed problems, the Hessian matrix of the data misfit part has low numerical rank and it is therefore possible to perform a low-rank approach to approximate the posterior covariance matrix. For such a low-rank approximation, one needs to solve a forward partial differential equation (PDE) and the adjoint PDE in both space and time. This in turn gives $\mathcal{O}(n_x n_t)$ complexity for both, computation and storage, where $n_x$ is the dimension of the spatial domain and $n_t$ is the dimension of the time domain. Such computations and storage demand are infeasible for large problems. To overcome this obstacle, we develop a new approach that utilizes a recently developed low-rank in time algorithm together with the low-rank Hessian method. We reduce both the computational complexity and storage requirement from $\mathcal{O}(n_x n_t)$ to $\mathcal{O}(n_x + n_t)$. We use numerical experiments to illustrate the advantages of our approach.

Igor Pontes Duff, Sara Grundel, Peter Benner

In this paper, we propose new algebraic Gramians for continuous-time linear switched systems, which satisfy generalized Lyapunov equations. The main contribution of this work is twofold. First, we show that the ranges of those Gramians encode the reachability and observability spaces of a linear switched system. As a consequence, a simple Gramian-based criterion for reachability and observability is established. Second, a balancing-based model order reduction technique is proposed and, under some sufficient conditions, stability preservation and an error bound are shown. Finally, the efficiency of the proposed method is illustrated by means of numerical examples.

Peter Benner, Venera Khoromskaia, Boris N. Khoromskij

We introduce and analyze the new range-separated (RS) canonical/Tucker tensor format which aims for numerical modeling of the 3D long-range interaction potentials in multi-particle systems. The main idea of the RS tensor format is the independent grid-based low-rank representation of the localized and global parts in the target tensor which allows the efficient numerical approximation of $N$-particle interaction potentials. The single-particle reference potential like $1/\|x\|$ is split into a sum of localized and long-range low-rank canonical tensors represented on a fine 3D $n\times n\times n$ Cartesian grid. The smoothed long-range contribution to the total potential sum is represented on the 3D grid in $O(n)$ storage via the low-rank canonical/Tucker tensor. We prove that the Tucker rank parameters depend only logarithmically on the number of particles $N$ and the grid-size $n$. Agglomeration of the short range part in the sum is reduced to an independent treatment of $N$ localized terms with almost disjoint effective supports, calculated in $O(N)$ operations. Thus, the cumulated sum of short range clusters is parametrized by a single low-rank canonical reference tensor with a local support, accomplished by a list of particle coordinates and their charges. The RS canonical/Tucker tensor representations reduce the cost of multi-linear algebraic operations on the 3D potential sums arising in modeling of multi-dimensional data by radial basis functions, say, in computation of the electrostatic potential of a protein, in 3D integration and convolution transforms, computation of gradients, forces and the interaction energy of a many-particle systems, and in low parametric fitting of multi-dimensional scattered data by reducing all of them to 1D calculations.

Pawan Goyal, Peter Benner

Measurement noise is an integral part while collecting data of a physical process. Thus, noise removal is necessary to draw conclusions from these data, and it often becomes essential to construct dynamical models using these data. We discuss a methodology to learn differential equation(s) using noisy and irregular sampled measurements. In our methodology, the main innovation can be seen in the integration of deep neural networks with the neural ordinary differential equations (ODEs) approach. Precisely, we aim at learning a neural network that provides (approximately) an implicit representation of the data and an additional neural network that models the vector fields of the dependent variables. We combine these two networks by constraining using neural ODEs. The proposed framework to learn a model describing the vector field is highly effective under noisy measurements. The approach can handle scenarios where dependent variables are not available at the same temporal grid. Moreover, a particular structure, e.g., second-order with respect to time, can easily be incorporated. We demonstrate the effectiveness of the proposed method for learning models using data obtained from various differential equations and present a comparison with the neural ODE method that does not make any special treatment to noise.

Pawan Goyal, Süleyman Yıldız, Peter Benner

Discovering a suitable coordinate transformation for nonlinear systems enables the construction of simpler models, facilitating prediction, control, and optimization for complex nonlinear systems. To that end, Koopman operator theory offers a framework for global linearization for nonlinear systems, thereby allowing the usage of linear tools for design studies. In this work, we focus on the identification of global linearized embeddings for canonical nonlinear Hamiltonian systems through a symplectic transformation. While this task is often challenging, we leverage the power of deep learning to discover the desired embeddings. Furthermore, to overcome the shortcomings of Koopman operators for systems with continuous spectra, we apply the lifting principle and learn global cubicized embeddings. Additionally, a key emphasis is paid to enforce the bounded stability for the dynamics of the discovered embeddings. We demonstrate the capabilities of deep learning in acquiring compact symplectic coordinate transformation and the corresponding simple dynamical models, fostering data-driven learning of nonlinear canonical Hamiltonian systems, even those with continuous spectra.

Pawan Goyal, Igor Pontes Duff, Peter Benner

Machine-learning technologies for learning dynamical systems from data play an important role in engineering design. This research focuses on learning continuous linear models from data. Stability, a key feature of dynamic systems, is especially important in design tasks such as prediction and control. Thus, there is a need to develop methodologies that provide stability guarantees. To that end, we leverage the parameterization of stable matrices proposed in [Gillis/Sharma, Automatica, 2017] to realize the desired models. Furthermore, to avoid the estimation of derivative information to learn continuous systems, we formulate the inference problem in an integral form. We also discuss a few extensions, including those related to control systems. Numerical experiments show that the combination of a stable matrix parameterization and an integral form of differential equations allows us to learn stable systems without requiring derivative information, which can be challenging to obtain in situations with noisy or limited data.

Peter Benner, Pawan Goyal

We discuss model reduction for a particular class of quadratic-bilinear (QB) descriptor systems. The main goal of this article is to extend the recently studied interpolation-based optimal model reduction framework for QBODEs [Benner et al. '16] to a class of descriptor systems in an efficient and reliable way. Recently, it has been shown in the case of linear or bilinear systems that a direct extension of interpolation-based model reduction techniques to descriptor systems, without any modifications, may lead to poor reduced-order systems. Therefore, for the analysis, we aim at transforming the considered QB descriptor system into an equivalent QBODE system by means of projectors for which standard model reduction techniques for QBODEs can be employed, including aforementioned interpolation scheme. Subsequently, we discuss related computational issues, thus resulting in a modified algorithm that allows us to construct \emph{near}--optimal reduced-order systems without explicitly computing the projectors used in the analysis. The efficiency of the proposed algorithm is illustrated by means of a numerical example, obtained via semi-discretization of the Navier-Stokes equations.

Xin Du, Peter Benner

Balanced truncation is one of the most common model order reduction schemes. In this paper, we study finite-frequency model order reduction (FF-MOR) problems of linear continuous-time systems within the framework of balanced truncation method. Firstly, we construct a family of parameterized frequency-dependent (PFD) mappings which generate discrete-time PFD mapped systems and continuous-time PFD mapped systems of the given continuous-time system. The relationships between the maximum singular value of the given system over pre-specified frequency ranges and the maximum singular value of the PFD mapped systems over entire frequency range are established. By exploiting the properties of the discrete-time PFD mapped systems, a new parameterized frequency-dependent balanced truncation (PFDBT) method providing finite-frequency type error bound with respect to the maximum singular value of the approximation error systems are developed. Examples are included for illustration.

Peter Benner, Boris N. Khoromskij, Venera Khoromskaia et al.

We introduce and analyze a mesh-free two-level hybrid Chebyshev-Tucker tensor representation for approximating multivariate functions, which combines tensor-product Chebyshev interpolation with the low-rank Tucker decomposition of the tensor of Chebyshev coefficients. This construction allows to avoid the expensive rank-structured grid-based approximation of function-related tensors on large spatial grids, while benefiting from the Tucker decomposition of the moderate-sized core tensor of Chebyshev coefficients. Thus, we can compute the nearly optimal Tucker decomposition of the 3D function with controllable accuracy $\varepsilon >0$ without discretizing the function on a full fine grid in the domain, but only using its values at a small set of Chebyshev nodes computed either from the explicit analytic expression of the target function or from its data-sparse representation in a rank-structured tensor format with moderate rank parameter. Finally, we can represent the function in the algebraic Tucker format with optimal $\varepsilon$-rank on an arbitrarily large 3D tensor grid in the computational domain by discretizing the Chebyshev polynomials on that grid. The rank parameters of the nonlinear Tucker-ALS decomposition of the coefficient tensor can be much smaller than the polynomial degrees of the initial Chebyshev linear interpolation in the function independent polynomial basis set. It is shown that our techniques can be gainfully applied to the long-range part of the singular electrostatic potential of multi-particle systems represented on a fine grid in the range-separated (RS) tensor format. We provide error and complexity estimates and demonstrate the computational efficiency of the proposed techniques on challenging examples, including the collective electrostatic potential for large bio-molecular systems and lattice-type compounds.

Peter Benner, Venera Khoromskaia, Boris N. Khoromskij et al.

In this paper, we introduce a new interpolation scheme to approximate the density of states (DOS) for a class of rank-structured matrices with application to the Tamm-Dancoff approximation (TDA) of the Bethe-Salpeter equation (BSE). The presented approach for approximating the DOS is based on two main techniques. First, we propose an economical method for calculating the traces of parametric matrix resolvents at interpolation points by taking advantage of the block-diagonal plus low-rank matrix structure described in [6, 3] for the BSE/TDA problem. Second, we show that a regularized or smoothed DOS discretized on a fine grid of size $N$ can be accurately represented by a low rank quantized tensor train (QTT) tensor that can be determined through a least squares fitting procedure. The latter provides good approximation properties for strictly oscillating DOS functions with multiple gaps, and requires asymptotically much fewer ($O(\log N)$) functional calls compared with the full grid size $N$. This approach allows us to overcome the computational difficulties of the traditional schemes by avoiding both the need of stochastic sampling and interpolation by problem independent functions like polynomials etc. Numerical tests indicate that the QTT approach yields accurate recovery of DOS associated with problems that contain relatively large spectral gaps. The QTT tensor rank only weakly depends on the size of a molecular system which paves the way for treating large-scale spectral problems.

Tobias Damm, Peter Benner, Jan Hauth

The stochastic $H^\infty$-norm is defined as the $L^2$-induced norm of the input-output operator of a stochastic linear system. Like the deterministic $H^\infty$-norm it is characterised by a version of the bounded real lemma, but without a frequency domain description or a Hamiltonian condition. Therefore, we base its computation on a parametrised algebraic Riccati-type matrix equation.

Kirandeep Kour, Sergey Dolgov, Peter Benner et al.

High-dimensional data in the form of tensors are challenging for kernel classification methods. To both reduce the computational complexity and extract informative features, kernels based on low-rank tensor decompositions have been proposed. However, what decisive features of the tensors are exploited by these kernels is often unclear. In this paper we propose a novel kernel that is based on the Tucker decomposition. For this kernel the Tucker factors are computed based on re-weighting of the Tucker matrices with tuneable powers of singular values from the HOSVD decomposition. This provides a mechanism to balance the contribution of the Tucker core and factors of the data. We benchmark support tensor machines with this new kernel on several datasets. First we generate synthetic data where two classes differ in either Tucker factors or core, and compare our novel and previously existing kernels. We show robustness of the new kernel with respect to both classification scenarios. We further test the new method on real-world datasets. The proposed kernel has demonstrated a higher test accuracy than the state-of-the-art tensor train multi-way multi-level kernel, and a significantly lower computational time.

Ali Forootani, Pawan Goyal, Peter Benner

The discovery of governing equations from data has been an active field of research for decades. One widely used methodology for this purpose is sparse regression for nonlinear dynamics, known as SINDy. Despite several attempts, noisy and scarce data still pose a severe challenge to the success of the SINDy approach. In this work, we discuss a robust method to discover nonlinear governing equations from noisy and scarce data. To do this, we make use of neural networks to learn an implicit representation based on measurement data so that not only it produces the output in the vicinity of the measurements but also the time-evolution of output can be described by a dynamical system. Additionally, we learn such a dynamic system in the spirit of the SINDy framework. Leveraging the implicit representation using neural networks, we obtain the derivative information -- required for SINDy -- using an automatic differentiation tool. To enhance the robustness of our methodology, we further incorporate an integral condition on the output of the implicit networks. Furthermore, we extend our methodology to handle data collected from multiple initial conditions. We demonstrate the efficiency of the proposed methodology to discover governing equations under noisy and scarce data regimes by means of several examples and compare its performance with existing methods.

Mohammad S. Khorrami, Pawan Goyal, Jaber R. Mianroodi et al.

The purpose of the current work is the development of a so-called physics-encoded Fourier neural operator (PeFNO) for surrogate modeling of the quasi-static equilibrium stress field in solids. Rather than accounting for constraints from physics in the loss function as done in the (now standard) physics-informed approach, the physics-encoded approach incorporates or "encodes" such constraints directly into the network or operator architecture. As a result, in contrast to the physics-informed approach in which only training is physically constrained, both training and output are physically constrained in the physics-encoded approach. For the current constraint of divergence-free stress, a novel encoding approach based on a stress potential is proposed. As a "proof-of-concept" example application of the proposed PeFNO, a heterogeneous polycrystalline material consisting of isotropic elastic grains subject to uniaxial extension is considered. Stress field data for training are obtained from the numerical solution of a corresponding boundary-value problem for quasi-static mechanical equilibrium. This data is also employed to train an analogous physics-guided FNO (PgFNO) and physics-informed FNO (PiFNO) for comparison. As confirmed by this comparison and as expected on the basis of their differences, the output of the trained PeFNO is significantly more accurate in satisfying mechanical equilibrium than the output of either the trained PgFNO or the trained PiFNO.

Peter Benner, Xiaobo Liu

We develop a mixed-precision iterative refinement framework for solving low-rank Lyapunov matrix equations $AX + XA^T + W =0$, where $W=LL^T$ or $W=LSL^T$. Via rounding error analysis of the algorithms we derive sufficient conditions for the attainable normwise residuals in different precision settings and show how the algorithmic parameters should be chosen. These conditions are independent of the choice of inner solver, provided that the prescribed residual accuracy is attained in the inner solves. Using the sign-function Newton iteration as the solver, we demonstrate that reduced precisions, such as half precision with unit roundoff $u_s$, can be used efficiently for Lyapunov equations with condition numbers of order $1/u_s$ without compromising the attainable solution quality. This provides an algorithmic framework towards exploiting native low-precision hardware to accelerate Lyapunov solvers without sacrificing accuracy.

Süleyman Yıldız, Pawan Goyal, Peter Benner

This paper presents an energy-preserving machine learning method for inferring reduced-order models (ROMs) by exploiting the multi-symplectic form of partial differential equations (PDEs). The vast majority of energy-preserving reduced-order methods use symplectic Galerkin projection to construct reduced-order Hamiltonian models by projecting the full models onto a symplectic subspace. However, symplectic projection requires the existence of fully discrete operators, and in many cases, such as black-box PDE solvers, these operators are inaccessible. In this work, we propose an energy-preserving machine learning method that can infer the dynamics of the given PDE using data only, so that the proposed framework does not depend on the fully discrete operators. In this context, the proposed method is non-intrusive. The proposed method is grey box in the sense that it requires only some basic knowledge of the multi-symplectic model at the partial differential equation level. We prove that the proposed method satisfies spatially discrete local energy conservation and preserves the multi-symplectic conservation laws. We test our method on the linear wave equation, the Korteweg-de Vries equation, and the Zakharov-Kuznetsov equation. We test the generalization of our learned models by testing them far outside the training time interval.

Celine Reddig, Pawan Goyal, Igor Pontes Duff et al.

Model reduction is an active research field to construct low-dimensional surrogate models of high fidelity to accelerate engineering design cycles. In this work, we investigate model reduction for linear structured systems using dominant reachable and observable subspaces. When the training set $-$ containing all possible interpolation points $-$ is large, then these subspaces can be determined by solving many large-scale linear systems. However, for high-fidelity models, this easily becomes computationally intractable. To circumvent this issue, in this work, we propose an active sampling strategy to sample only a few points from the given training set, which can allow us to estimate those subspaces accurately. To this end, we formulate the identification of the subspaces as the solution of the generalized Sylvester equations, guiding us to select the most relevant samples from the training set to achieve our goals. Consequently, we construct solutions of the matrix equations in low-rank forms, which encode subspace information. We extensively discuss computational aspects and efficient usage of the low-rank factors in the process of obtaining reduced-order models. We illustrate the proposed active sampling scheme to obtain reduced-order models via dominant reachable and observable subspaces and present its comparison with the method where all the points from the training set are taken into account. It is shown that the active sample strategy can provide us $17$x speed-up without sacrificing any noticeable accuracy.

Peter Benner, Volker Mehrmann, Anshul Prajapati et al.

We study linear time-invariant dissipative Hamiltonian differential-algebraic systems. We characterize when the systems are robustly asymptotically stable and derive exact conditions and bounds when this property is lost under structure-preserving perturbations.

Mohammad S. Khorrami, Pawan Goyal, Soroush Motahari et al.

The purpose of the current work is the development of an approach to account for quasi-static mechanical equilibrium in empirical (i.e., data-based) models for the stress field employing neural approximations (NAs), which include neural networks (NNs) and neural operators (NOs), in particular Fourier NOs (FNOs). Rather than including such constraints from physics in the loss function as done in the (now standard) physics-informed approach, the current approach incorporates or "encodes" such constraints directly into the architecture of the NA. As a result, both NA training and output are physically constrained in the physics-encoded approach, in contrast to the physics-informed approach, in which only training is physically constrained. For the current constraint of divergence-free stress, a novel encoding approach based on a stress potential is proposed. As a "proof-of-concept" example application of the current approach, a physics-encoded FNO (PeFNO) is developed for a heterogeneous polycrystalline material consisting of isotropic elastic grains and subject to uniaxial extension. Stress field data for this purpose are obtained from the numerical solution of corresponding boundary-value problems for quasi-static mechanical equilibrium. For comparison with the PeFNO, this data is also employed to develop an analogous physics-guided FNO (PgFNO) and physics-informed FNO (PiFNO). As expected theoretically, and confirmed by this computational comparison, for comparable accuracy of the stress field itself as compared to the data, the stress field output by the trained and tested PeFNO is significantly more accurate in satisfying mechanical equilibrium than the output of either the PgFNO or the PiFNO.

Ali Forootani, Harshit Kapadia, Sridhar Chellappa et al.

The sparse identification of nonlinear dynamical systems (SINDy) is a data-driven technique employed for uncovering and representing the fundamental dynamics of intricate systems based on observational data. However, a primary obstacle in the discovery of models for nonlinear partial differential equations (PDEs) lies in addressing the challenges posed by the curse of dimensionality and large datasets. Consequently, the strategic selection of the most informative samples within a given dataset plays a crucial role in reducing computational costs and enhancing the effectiveness of SINDy-based algorithms. To this aim, we employ a greedy sampling approach to the snapshot matrix of a PDE to obtain its valuable samples, which are suitable to train a deep neural network (DNN) in a SINDy framework. SINDy based algorithms often consist of a data collection unit, constructing a dictionary of basis functions, computing the time derivative, and solving a sparse identification problem which ends to regularised least squares minimization. In this paper, we extend the results of a SINDy based deep learning model discovery (DeePyMoD) approach by integrating greedy sampling technique in its data collection unit and new sparsity promoting algorithms in the least squares minimization unit. In this regard we introduce the greedy sampling neural network in sparse identification of nonlinear partial differential equations (GN-SINDy) which blends a greedy sampling method, the DNN, and the SINDy algorithm. In the implementation phase, to show the effectiveness of GN-SINDy, we compare its results with DeePyMoD by using a Python package that is prepared for this purpose on numerous PDE discovery

Igor Pontes Duff, Pawan Goyal, Peter Benner

This work primarily focuses on an operator inference methodology aimed at constructing low-dimensional dynamical models based on a priori hypotheses about their structure, often informed by established physics or expert insights. Stability is a fundamental attribute of dynamical systems, yet it is not always assured in models derived through inference. Our main objective is to develop a method that facilitates the inference of quadratic control dynamical systems with inherent stability guarantees. To this aim, we investigate the stability characteristics of control systems with energy-preserving nonlinearities, thereby identifying conditions under which such systems are bounded-input bounded-state stable. These insights are subsequently applied to the learning process, yielding inferred models that are inherently stable by design. The efficacy of our proposed framework is demonstrated through a couple of numerical examples.

Konrad Janik, Peter Benner

Measurement data is often sampled irregularly i.e. not on equidistant time grids. This is also true for Hamiltonian systems. However, existing machine learning methods, which learn symplectic integrators, such as SympNets [20] and HénonNets [4] still require training data generated by fixed step sizes. To learn time-adaptive symplectic integrators, an extension to SympNets, which we call TSympNets, was introduced in [20]. We adapt the architecture of TSympNets and extend them to non-autonomous Hamiltonian systems. So far the approximation qualities of TSympNets were unknown. We close this gap by providing a universal approximation theorem for separable Hamiltonian systems and show that it is not possible to extend it to non-separable Hamiltonian systems. To investigate these theoretical approximation capabilities, we perform different numerical experiments. Furthermore we fix a mistake in a proof of a substantial theorem [25, Theorem 2] for the approximation of symplectic maps in general, but specifically for symplectic machine learning methods.

Shuwen Sun, Lihong Feng, Peter Benner

Numerically solving a large parametric nonlinear dynamical system is challenging due to its high complexity and the high computational costs. In recent years, machine-learning-aided surrogates are being actively researched. However, many methods fail in accurately generalizing in the entire time interval $[0, T]$, when the training data is available only in a training time interval $[0, T_0]$, with $T_0<T$. To improve the extrapolation capabilities of the surrogate models in the entire time domain, we propose a new deep learning framework, where kernel dynamic mode decomposition (KDMD) is employed to evolve the dynamics of the latent space generated by the encoder part of a convolutional autoencoder (CAE). After adding the KDMD-decoder-extrapolated data into the original data set, we train the CAE along with a feed-forward deep neural network using the augmented data. The trained network can predict future states outside the training time interval at any out-of-training parameter samples. The proposed method is tested on two numerical examples: a FitzHugh-Nagumo model and a model of incompressible flow past a cylinder. Numerical results show accurate and fast prediction performance in both the time and the parameter domain.

Ion Victor Gosea, Luisa Peterson, Pawan Goyal et al.

In this work, we address the challenge of efficiently modeling dynamical systems in process engineering. We use reduced-order model learning, specifically operator inference. This is a non-intrusive, data-driven method for learning dynamical systems from time-domain data. The application in our study is carbon dioxide methanation, an important reaction within the Power-to-X framework, to demonstrate its potential. The numerical results show the ability of the reduced-order models constructed with operator inference to provide a reduced yet accurate surrogate solution. This represents an important milestone towards the implementation of fast and reliable digital twin architectures.

Peter Benner, Pascal den Boef, Patrick Kürschner et al.

We investigate the acceleration of stationary iterations for multi-term Sylvester equation by means of reduced rank extrapolation (RRE). Theoretical convergence results and implementations are provided for both small and large-scale problems. For the large-scale problems, an inexact non-stationary iteration is discussed, which makes use of low-rank matrix approximations. Numerical experiments illustrate the potential of the RRE acceleration which often leads to a substantial gain in convergence speed and therefore reducing the consumption of storage and computing time.

Konrad Janik, Peter Benner

Measurement data is often sampled irregularly, i.e., not on equidistant time grids. This is also true for Hamiltonian systems. However, existing machine learning methods, which learn symplectic integrators, such as SympNets [1] and HénonNets [2] still require training data generated by fixed step sizes. To learn time-adaptive symplectic integrators, an extension to SympNets called TSympNets is introduced in [3]. The aim of this work is to do a similar extension for HénonNets. We propose a novel neural network architecture called T-HénonNets, which is symplectic by design and can handle adaptive time steps. We also extend the T-HénonNet architecture to non-autonomous Hamiltonian systems. Additionally, we provide universal approximation theorems for both new architectures for separable Hamiltonian systems and discuss why it is difficult to handle non-separable Hamiltonian systems with the proposed methods. To investigate these theoretical approximation capabilities, we perform different numerical experiments.

Süleyman Yıldız, Konrad Janik, Peter Benner

We propose a new symplectic convolutional neural network (CNN) architecture by leveraging symplectic neural networks, proper symplectic decomposition, and tensor techniques. Specifically, we first introduce a mathematically equivalent form of the convolution layer and then, using symplectic neural networks, we demonstrate a way to parameterize the layers of the CNN to ensure that the convolution layer remains symplectic. To construct a complete autoencoder, we introduce a symplectic pooling layer. We demonstrate the performance of the proposed neural network on three examples: the wave equation, the nonlinear Schrödinger (NLS) equation, and the sine-Gordon equation. The numerical results indicate that the symplectic CNN outperforms the linear symplectic autoencoder obtained via proper symplectic decomposition.

Harshit Kapadia, Peter Benner, Lihong Feng

In situations where the solution of a high-fidelity dynamical system needs to be evaluated repeatedly, over a vast pool of parametric configurations and in absence of access to the underlying governing equations, data-driven model reduction techniques are preferable. We propose a novel active learning approach to build a parametric data-driven reduced-order model (ROM) by greedily picking the most important parameter samples from the parameter domain. As a result, during the ROM construction phase, the number of high-fidelity solutions dynamically grow in a principled fashion. The high-fidelity solution snapshots are expressed in several parameter-specific linear subspaces, with the help of proper orthogonal decomposition (POD), and the relative distance between these subspaces is used as a guiding mechanism to perform active learning. For successfully achieving this, we provide a distance measure to evaluate the similarity between pairs of linear subspaces with different dimensions, and also show that this distance measure is a metric. The usability of the proposed subspace-distance-enabled active learning (SDE-AL) framework is demonstrated by augmenting two existing non-intrusive reduced-order modeling approaches, and providing their active-learning-driven (ActLearn) extensions, namely, SDE-ActLearn-POD-KSNN, and SDE-ActLearn-POD-NN. Furthermore, we report positive results for two parametric physical models, highlighting the efficiency of the proposed SDE-AL approach.

Shuwen Sun, Lihong Feng, Peter Benner

We explore the promising performance of a transformer model in predicting outputs of parametric dynamical systems with external time-varying input signals. The outputs of such systems vary not only with physical parameters but also with external time-varying input signals. Accurately catching the dynamics of such systems is challenging. We have adapted and extended an existing transformer model for single output prediction to a multiple-output transformer that is able to predict multiple output responses of these systems. The multiple-output transformer generalizes the interpretability of the original transformer. The generalized interpretable attention weight matrix explores not only the temporal correlations in the sequence, but also the interactions between the multiple outputs, providing explanation for the spatial correlation in the output domain. This multiple-output transformer accurately predicts the sequence of multiple outputs, regardless of the nonlinearity of the system and the dimensionality of the parameter space.

Pawan Goyal, Peter Benner

Learning dynamical models from data plays a vital role in engineering design, optimization, and predictions. Building models describing dynamics of complex processes (e.g., weather dynamics, or reactive flows) using empirical knowledge or first principles are onerous or infeasible. Moreover, these models are high-dimensional but spatially correlated. It is, however, observed that the dynamics of high-fidelity models often evolve in low-dimensional manifolds. Furthermore, it is also known that for sufficiently smooth vector fields defining the nonlinear dynamics, a quadratic model can describe it accurately in an appropriate coordinate system, conferring to the McCormick relaxation idea in nonconvex optimization. Here, we aim at finding a low-dimensional embedding of high-fidelity dynamical data, ensuring a simple quadratic model to explain its dynamics. To that aim, this work leverages deep learning to identify low-dimensional quadratic embeddings for high-fidelity dynamical systems. Precisely, we identify the embedding of data using an autoencoder to have the desired property of the embedding. We also embed a Runge-Kutta method to avoid the time-derivative computations, which is often a challenge. We illustrate the ability of the approach by a couple of examples, arising in describing flow dynamics and the oscillatory tubular reactor model.

Pawan Goyal, Peter Benner

Measurement noise is an integral part while collecting data of a physical process. Thus, noise removal is a necessary step to draw conclusions from these data, and it often becomes quite essential to construct dynamical models using these data. We discuss a methodology to learn differential equation(s) using noisy and sparsely sampled measurements. In our methodology, the main innovation can be seen in of integration of deep neural networks with a classical numerical integration method. Precisely, we aim at learning a neural network that implicitly represents the data and an additional neural network that models the vector fields of the dependent variables. We combine these two networks by enforcing the constraint that the data at the next time-steps can be given by following a numerical integration scheme such as the fourth-order Runge-Kutta scheme. The proposed framework to learn a model predicting the vector field is highly effective under noisy measurements. The approach can handle scenarios where dependent variables are not available at the same temporal grid. We demonstrate the effectiveness of the proposed method to learning models using data obtained from various differential equations. The proposed approach provides a promising methodology to learn dynamic models, where the first-principle understanding remains opaque.

Pawan Goyal, Peter Benner

Discovering dynamical models to describe underlying dynamical behavior is essential to draw decisive conclusions and engineering studies, e.g., optimizing a process. Experimental data availability notwithstanding has increased significantly, but interpretable and explainable models in science and engineering yet remain incomprehensible. In this work, we blend machine learning and dictionary-based learning with numerical analysis tools to discover governing differential equations from noisy and sparsely-sampled measurement data. We utilize the fact that given a dictionary containing huge candidate nonlinear functions, dynamical models can often be described by a few appropriately chosen candidates. As a result, we obtain interpretable and parsimonious models which are prone to generalize better beyond the sampling regime. Additionally, we integrate a numerical integration framework with dictionary learning that yields differential equations without requiring or approximating derivative information at any stage. Hence, it is utterly effective in corrupted and sparsely-sampled data. We discuss its extension to governing equations, containing rational nonlinearities that typically appear in biological networks. Moreover, we generalized the method to governing equations that are subject to parameter variations and externally controlled inputs. We demonstrate the efficiency of the method to discover a number of diverse differential equations using noisy measurements, including a model describing neural dynamics, chaotic Lorenz model, Michaelis-Menten Kinetics, and a parameterized Hopf normal form.

Pawan Goyal, Peter Benner

Mathematical modeling is an essential step, for example, to analyze the transient behavior of a dynamical process and to perform engineering studies such as optimization and control. With the help of first-principles and expert knowledge, a dynamic model can be built, but for complex dynamic processes, appearing, e.g., in biology, chemical plants, neuroscience, financial markets, this often remains an onerous task. Hence, data-driven modeling of the dynamics process becomes an attractive choice and is supported by the rapid advancement in sensor and measurement technology. A data-driven approach, namely operator inference framework, models a dynamic process, where a particular structure of the nonlinear term is assumed. In this work, we suggest combining the operator inference with certain deep neural network approaches to infer the unknown nonlinear dynamics of the system. The approach uses recent advancements in deep learning and possible prior knowledge of the process if possible. We also briefly discuss several extensions and advantages of the proposed methodology. We demonstrate that the proposed methodology accomplishes the desired tasks for dynamics processes encountered in neural dynamics and the glycolytic oscillator.