Preserving the DAE structure in the Loewner model reduction and identification framework

We propose an extension of the Loewner framework to descriptor linear systems that preserves the DAE (differential algebraic equation) structure of the underlying system. More precisely, by means of post-processing the data, the behavior at infinity is matched. As it turns out, the conventional procedure constructs a reduced model by directly compressing the data and hence losing information at infinity. By transforming the matrix pencil composed of the E and A matrices into a generalized block diagonal form, we can separate the descriptor system into two subsystems; one corresponding to the polynomial part and the other to the strictly proper part of the transfer function. Different algorithms are implemented to transform the matrix pencil into block diagonal form. Furthermore, a data-driven splitting of the descriptor system can be achieved in the Loewner framework. Hence, the coefficients of the polynomial part can be estimated directly from data. Several numerical examples are presented to illustrate the theoretical discussion.

faithfully reproduce the behavior or response of the original system. Such highdimensional models usually appear from spatial discretization of partial differential equations (PDEs). The denser the grid used, the higher the number of equations that characterize the behavior. Reduced-order systems are described by significantly lower number of equations. Such systems could then be efficiently used as surrogates, replacing the large-scale systems in such tasks as simulation, analysis, or control. For details on different MOR techniques, we refer the reader to the books [1,4,12] and to the surveys [3,7,10].
In the current work, we will turn our attention to interpolatory model reduction methods. In this setup, ROMs (reduced-order models) are constructed by means of interpolating the transfer function of the original system at selected points.
The Loewner framework, as introduced in [34], will be the main focus of this study. It is a data-driven method that directly constructs a ROM using only measured data (samples of the transfer function). Additionally, by compressing the (usually large) data set, it extracts the dominant features. For a recent survey on the Loewner framework for linear systems, see [6]. For the nonlinear case, this method requires appropriate definition of transfer functions, as described in some of its extensions, i.e., to bilinear [5], quadratic-bilinear [23], and switched systems [24].
We propose an extension of the Loewner framework to a specific class of linear systems, i.e., descriptor linear systems. The new method can preserve the underlying DAE structure by means of post-processing of the data (estimating the polynomial part of the transfer function).
Whenever the dynamics of the driving system is described by not only differential but also algebraic equations, descriptor systems come into play. Such systems are described by DAEs and are used to model many different real-world applications, such as the ones appearing in electrical networks/circuits, biomedical and chemical processes, or computational fluid dynamics. We refer the reader to the books [22,33] for more details. Additionally, the following surveys address particular topics related to systems described by linear DAEs: partial stabilization, in [9], description of the solution in [39], controllability in [16], and numerical linear algebra methods in [11]. The concept of partial realization is covered in [13], while topics on simulation and control in [36].
Modeling of complex physical and technical processes such as control of laminar/turbulent fluid flow, very large integrated circuit simulation for chip design or components coupling in mechanical systems, typically leads to descriptor systems of considerably large dimension. This generates bottlenecks such as storage requirements and expensive computations when dealing with this class of systems. This is where reduced models play an important role. Consequently, in recent years, many such methods for reducing conventional linear systems (with no algebraic constraints) have been extended for reduction of systems described by DAEs. In particular, balanced truncation has been extended to descriptor systems in [26,35,38] while interpolation-based methods were considered in [25,37]. For more details on MOR methods applied to DAEs, we refer the readers to the extensive survey in [14].
The paper is organized in the following way: in Section 2, we give a short overview on descriptor systems, introduce the Weierstrass and Kronecker canonical forms together with some relaxed versions thereof. Additionally, we discuss how to explicitly compute the polynomial coefficients of the transfer function. In Section 3, we present several methods for splitting the strictly proper and polynomial part of the transfer function, for both systems involving regular and singular pencils. Section 4 introduces the classical Loewner framework and afterwards, the main proposed procedure that is based on it. The difference is that, by means of post-processing of the data, we are able to match the behavior at infinity of the underlying system. In Section 5, several numerical experiments are presented (three of them chosen from benchmark examples commonly used for MOR applications), while a summary of the findings and the conclusion are stated in Section 6.

Setup
We study linear time-invariant continuous-time descriptor systems described by the following equations:

Σ :
Eẋ(t) = Ax(t) + Bu(t), where A, E ∈ R n×n , C ∈ R p×n , B ∈ R n×m , and x(t) ∈ R n , u(t) ∈ R m , y(t) ∈ R p are the internal variable, input and output vectors, respectively. The order of system Σ in (1) is given by the number of internal variables n. If the E matrix is the identity matrix, then the system in (1) is a standard state space system. Otherwise, it is a descriptor system or generalized state space system. The goal is to approximate the original descriptor system Σ in (1) by a reducedorder linear time-invariant continuous-time descriptor systemΣ Σ : Êx (t) =Â(t) +Bu(t), whereÂ,Ê ∈ R r×r ,Ĉ ∈ R p×r ,B ∈ R r×m , andx(t) ∈ R r is the reduced internal variable, with r n. The control input u(t) is the same for both systems in (1) and (2). In general, it is required that the approximation error be small.
The input-output behavior of the descriptor system Σ is usually described by its transfer function H. Assuming that the pencil λE−A is regular, i.e., det(λE−A) = 0 for λ ∈ C, then one can write the transfer function H as follows: If the pencil λE − A is singular (not regular), i.e., det(λE − A) = 0 for all λ ∈ C, one needs to replace the conventional inverse of the matrix sE − A in (3) with a pseudo-inverse. In this case, the transfer function H can be written where X # ∈ C u×v is the Moore-Penrose inverse of the matrix X ∈ C v×u with u, v ∈ N (for more details, see [2]).
For simplicity of exposition, in what follows, we will treat solely the single-input, single-output (SISO) case with m = p = 1. The multi-input, multi-output case is technically more involved but it is based on the same ideas.
The transfer function shown in (3) can be expressed as a rational function in the variable s ∈ C where a j , b i ∈ R, z is the number of zeros and q is the number of poles. When z > q, the transfer function H(s) is improper. When z ≤ q, H(s) is proper and when z < q, H(s) is strictly proper. If z ≥ q, the transfer function H(s) in (3) can be decomposed into a strictly proper part and a polynomial part, as where are the strictly proper part and polynomial part, respectively. In (7), the polynomial coefficients of H poly (s) are denoted with p k , for k ∈ {0, 1, . . . , z − q}. In what follows, we discuss how to determine the polynomial part H poly (s) of the descriptor system Σ and how to explicitly compute the coefficients p k .

The Weierstrass canonical forms
Let Σ be a descriptor system as in (1) and assume that the pencil λE − A is regular. In this case, the pencil λE − A can be reduced to the Weierstrass canonical form, defined next.

Definition 1
The regular matrix pencil λE − A can be written in the Weierstrass canonical form [35] given by where the matrices S, T ∈ C n×n are invertible, J f ∈ C n f ×n f and N ∈ C n ∞ ×n ∞ are matrices in Jordan canonical form and N is a nilpotent matrix. The positive integers n f and n ∞ are the dimensions of the deflating subspaces of the pencil corresponding to the finite and infinite eigenvalues, respectively. Furthermore, I n f ∈ R n f ×n f and I n ∞ ∈ R n ∞ ×n ∞ are identity matrices. Let ν be the nilpotency index of N, i.e., the smallest positive integer ν such that N ν = 0. The positive integer ν is sometimes referred to as the index of the pencil λE − A and also the Hessenberg index of the descriptor system Σ.
In [15], a so-called quasi-Weierstrass canonical form is introduced. It is weaker than the classical form although it also contains relevant information.
Next, introduce a generalized Weierstrass canonical form of the matrix pencil λE− A. Here, the transformed E and A are indeed block diagonal, but the composing blocks need not have any special structure.

Definition 2
The regular matrix pencil λE − A can be written in the generalized quasi-Weierstrass canonical form (or generalized block diagonal form) given by Proof In (7), since H ∞ (s) = H poly (s) = z−q k=0 p k s k , it directly follows that the kth derivative of H ∞ (s) can be written as By substituting s = 0 in (15), it automatically follows . Using the definition of H ∞ in equation (12), the kth derivative of H ∞ is d k So if k > 0, then and

The right matrix pencil disk function
In this section, we assume that the original system is a discrete-time descriptor system computed, for example, by means of applying finite difference approximation of the derivative of the continuous internal variable in (1). Let λE − A be a regular matrix pencil with no eigenvalues on the unit circle. The Weierstrass canonical form of λE− A is given by where J 0 ∈ C n k ×n k and J ∞ ∈ C (n−n k )×(n−n k ) .

Definition 3
The right matrix pencil disk function is defined The matrices P 0 and P ∞ define skew projections onto the right deflating subspace corresponding to the eigenvalues of λE − A inside the unit circle and respectively outside the unit circle. The disk function provides a tool for spectral decomposition along the unit circle. Splittings for other curves can be computed by applying a suitable conformal mapping to the pencil λE − A.
The method for the computation of the right matrix disk function can be found in [9,41]. The steps of this method are summarized in Algorithm 3 from [9] and in Algorithm 8 from [41].

The Kronecker canonical forms
Let Σ be a continuous-time descriptor system, as in (1), and assume that the pencil λE − A is singular, i.e., det(λE − A) = 0 for all λ ∈ C. In this case, the pencil λE − A can be reduced to the Kronecker canonical form (KCF).

Definition 4
The matrix pencil λE − A can be written in the Kronecker canonical form [40] where P, Q ∈ C n×n are invertible matrices. The building blocks present in (19) (by following the description in [40]) are as follows: As mentioned in [8,40], to compute the Kronecker canonical form, it is recommended to first compute a quasi-triangular form (also called the generalized Schur form) in order to avoid numerical issues: This form can be obtained using the unitary transformation matrices P, Q ∈ C n×n . Moreover, the block diagonal structure is as follows: 1. The square regular pencil λE f − A f contains the finite eigenvalues of λE − A.
2. The square regular pencil λE ∞ −A ∞ contains the infinite eigenvalues of λE−A.
3. The singular pencils λE η − A η and λE − A contain the Kronecker row and column structure, respectively.

Separating the strictly proper and polynomial part of the transfer function
In Section 2, we presented the theoretical aspects of decomposing the matrix pencil λE−A into two subsystems, corresponding to the strictly proper part and polynomial part of the transfer function H(s). Next, we present different numerical algorithms for computing the matrices U and V. These are used in order to construct a quasi-Weierstrass form. Four such algorithms are studied in this section. First, the Jordan chain (JC) method finds the transformation matrices by calculating the generalized eigenvectors corresponding to the infinite eigenvalues. Then, the ADTF1 (additive decomposition of a transfer function no. 1) and ADTF2 (additive decomposition of a transfer function no. 2) methods from [30] use the QZ decomposition to separate the descriptor system. Finally, a recent procedure based on the matrix disk function in [41] is mentioned.

Jordan chain method
This method requires computing matrices U ∞ ∈ C n×n ∞ and V ∞ ∈ C n×n ∞ so that their columns span the same subspace as the left and, respectively, right generalized eigenvectors corresponding to the eigenvalue at infinity of the pencil λE − A (or, equivalently, of eigenvalues at 0 of the pencil λA − E) where J ∈ C n ∞ ×n ∞ is a Jordan block with zero eigenvalues. Then, the matrices U f , V f ∈ C n×n f can be found by computing the null space of matrices V * ∞ A * and, respectively, of U * ∞ A: Then, construct the transformation matrices U = U f U ∞ and V = V f V ∞ by putting together the matrices in (21) and (22).

Additive decomposition of a transfer function no. 1
This approach was introduced in [30] Section 4.1. The first step is to apply the QZ (known also as generalized Schur) decomposition to the matrix pencil λE−A. This is done by imposing that the finite eigenvalues of this pencil are located in the upper-left corner while the remaining infinite eigenvalues are in the lower-right corner: Next, construct matrices T and S in the following way The transformed matrices E 1 = TQEZS, A 1 = TQAZS are thus block diagonal. Numerical methods for solving coupled Sylvester equations as in (24) were proposed in [27,29,31]. The left and right transformation matrices can be explicitly constructed as follows U = (TQ) * , V = ZS.

Additive decomposition of a transfer function no. 2
Again, as in Section 3.1.2, the QZ factorization is applied to the matrix pencil λE−A with two different orderings of the eigenvalues (as in Section 4.2 in [30]). This is done by imposing that the finite and infinite eigenvalues of the pencil are located in the upper-left corner and in the lower-right corner, respectively: The transformation matrices are

Additive decomposition of a transfer function no. 3
The decompositions covered in Sections 3.1.2 and 3.1.3 rely on computing one or two QZ factorizations (as described in [30]). Recently, a new method was proposed in [41] that overcomes the usage of such factorizations and, instead, uses the disk function method (briefly mentioned in Section 2.2.3) to construct projection matrices for subspace extraction.
The complete procedure is stated in Algorithm 10 from [41]. The first step is to apply the disk function method to the pencil λ(αA)−E, where the real positive scalar α is chosen so that α < max |λ| : Here, α is used as eigenvalue scaling such that the infinite eigenvalues become zeros and everything else is pushed outside the unit circle. Hence, the disk function method produces skew projections onto the corresponding deflating subspaces.

Singular pencils
The computation of the Kronecker canonical form (KCF) received considerable attention in the second half of the previous century. This follows from the importance of the KCF in the area of systems and control theory, i.e., realization theory and controllability. In what follows, we mention only a few of the contributions.
1. A generalization of the staircase algorithm to singular pencils using unitary equivalence transformations is proposed in [40].
2. The so-called AB algorithm for computing the right singular structure and the Jordan structure of the zero eigenvalue for a singular pencil is introduced in [32]. 3. An algorithm based on repeated generalized singular value decompositions (or more precisely, cosine-sine decompositions of partitioned orthonormal matrices), in short referred to as RGSVD, is presented in [28]. 4. A faster algorithm which requires O(m 2 n) operations when the original pencil is m × n is proposed in [8]. 5. Robust algorithms for implementing the generalized Schur decomposition of an arbitrary pencil A − λB are presented in [19]. 6. A procedure based on the so-called Wong sequences is proposed in order to derive a quasi-Kronecker form in [17].
The procedures mentioned above are quite involved and not straightforward to implement in a practical setup.
The toolbox GUPTRI (Generalized Upper TRIangular form) presented in [20] can be used to compute the generalized Schur decomposition of an arbitrary pencil A − λB (regular or singular), hence revealing the Kronecker structure of the pencil. The algorithms are implemented in Fortran (also in Matlab) and are available as an open source toolbox (see [21]).
Most of the numerical experiments presented in Section 5 are based on the following block diagonalization procedure (in Section 3.2.1). This is also used as a direct numerical tool for estimating the polynomial coefficients.

The proposed procedure
We are using a numerical procedure inspired from the decomposition that was briefly mentioned in Section 3.1.3. Originally proposed in [30], this procedure uses unitary transformations to block diagonalize a regular matrix pencil based on two re-orderings of the pencil's eigenvalues.
For the first re-ordering, the blocks that have eigenvalues belonging to a set Γ are placed in the lower-right corner of the transformed upper triangular matrices. For the second re-ordering, the order is reversed (the blocks that have eigenvalues belonging to the set Γ are placed in the upper-left corner).
The procedure is based on performing two QZ (or generalized Schur) decompositions in complex arithmetic for a regular matrix pencil (A, E).
We are dealing with singular pairs (A, E) that have eigenvalues at ∞. So instead of analyzing the spectrum of the pencil (A, E), we inspect the spectrum of the pencil (E, A). Hence, one quantity of interest is represented by the eigenvalues at 0 (instead of the eigenvalues at infinity). From a numerical point of view, this leads to a more stable procedure. The set of interest, denoted with Γ , contains the eigenvalues at 0.
Nevertheless, it is not straightforward to accurately identify the set Γ (as it will be made clear by the numerical examples). Since the pencil has singularities, the true zeros might get mixed with such singularities in the classification process. Hence, a reliable choice of k (that denotes the number of zeros) is sometimes not straightforward to find.
In what follows, we present the steps of the proposed procedure. 1). Initialization step: start with system matrices (C, E, A, B) that form a realization of a linear descriptor system as in (1). Then, let C 1 = C, E 1 = E, A 1 = A and B 1 = B.
2). Perform a QZ (generalized Schur) decomposition of (E 1 , A 1 ): so that the elements of the vector γ d is a small positive real number; in practical applications, it is chosen as the eps value in Matlab-the machine precision ≈ 2.2204 · 10 −16 . The purpose of using τ is to avoid division by 0. Note that, in (28) Next, for k ∈ N, we split the decomposition into two parts corresponding to blocks of dimension n − k, and respectively, k (with eigenvalues in the set Γ ), as follows: so that, this time, the elements of the vector γ a 1 γ a 2 . . . γ a n ∈ R + ∪ {0} are sorted in an ascending order, i.e., γ a 1 γ a 2 · · · γ a n , where γ a i = |α a i | |β a i |+τ . As before, α a i ∈ C and β a i ∈ C are the (i, i) entries of matrices E a 1 , and of A a 1 . The positive constant τ is chosen as before. Note that, in (30), matrices E a 1 , A a 1 are upper triangular while Q a , Z a are unitary.
Next, use the same value k ∈ N to split the decomposition into two parts corresponding to blocks of dimension k (with eigenvalues in the set Γ ), and respectively n − k, as follows:

3). Proceed by putting together the transformation matrices
It follows that the matrix pencil λA 1 − E 1 can be block diagonalized as Put together the transformed model as follows:

4).
Split the transformed linear model into two blocks: Use the second blocks to compute the polynomial coefficients: whereÃ # 2 represents the Moore-Penrose inverse ofÃ 2 .
Remark 1 Note that the procedure described above represents a direct extension of the ADTF2 algorithm from [30] to the case of singular pencils.

Remark 2
The values α d i and β d i that appear in the second step of the proposed procedure are indeed of relevance. The reason is that they indicate the blocks corresponding to singularities. The proposed procedure could be altered so that the criterion for constructing the set Γ changes. For example, instead of computing the values γ d i , one can solely analyze the values α d i and β d i .
Example 1 We choose matrices (C, E, A, B) that form a non-minimal linear descriptor realization of the transfer function H(s) = 2s + 3 + 4 s−5 , i.e., Note that the pencil λE − A is singular. Hence, the transfer function H of the system Additionally, note that rank(E) = 2, while rank(A) = 3. Hence, the pencil (E, A) has only one nonzero finite eigenvalue at 1 5 , two eigenvalues at 0 and one singularity. We proceed with the proposed splitting method by first computing the QZ factorization corresponding to the descending ordering. For example, the values α d i , β d i , and γ d i are presented in Table 1. Also, compute the QZ factorization corresponding to the ascending ordering. Next, for k = 2, we identify the matrices corresponding to the block diagonalized 0Ã 2 in the fourth step: By applying the formulas in (34), we correctly identify p 0 = 3, p 1 = 2, and p = 0 for 2. Next, we vary the value of k between 1 and 4 and collect the estimated polynomial coefficients for these four values. The results are presented in Table 2. Note that, for both k ∈ {2, 3}, the original coefficients are recovered.

The classical Loewner approach
In many situations, a model for the descriptor system Σ is not explicitly available; the system matrix realization in (1) can not be directly obtained. Instead of the system matrices A, B, C, and E, a black box is provided that produces input-output measurements. These situations include VLSI modeling from chips or real-time simulation of multi-body dynamics with constraints. The Loewner framework is used to construct a reduced-order system by means of the measured data. In this context, data consist of samples of the underlying transfer function H(s) at selected interpolation points (in some cases on the imaginary axis, when approximating the frequency response of the original system). A Loewner matrix pencil λL − L s is composed of the Loewner matrix L and of the shifted Loewner matrix L s . The Loewner framework originates from rational approximation theory (see [6,34]) and can be used to recover or reduce the original dynamical system by means of interpolating the transfer function. For simplicity, assume that the data set contains an even number of measurements denoted with 2N. The next step is to partition the data measurements set into two disjoint data subsets, i.e., the left subset and the right subset as described in the following Left set : (μ i , v i ), i = 1, . . . , N, Right set: (λ j , w j ), j = 1, . . . , N.
In (36), consider as given N distinct left interpolation points Note that the sets of left and right interpolation points are assumed distinct, i.e., μ i = λ j , ∀i, j = 1, . . . , N.
We seek a linear descriptor systemΣ such that the associated transfer func-tionĤ(s) is an interpolant to the original one, H(s). More precisely, the following interpolation conditions should then hold: Next, arrange the measured data into matrix format as follows: The Loewner matrix L ∈ C N×N and the shifted Loewner matrix L s ∈ C N×N are defined as follows: and are the solutions of the following Sylvester equation where L T = R = 1 1 · · · 1 ∈ C 1×N . After rearranging the data into matrix format, we construct with basically no computational cost a Loewner linear descriptor model described by the following realization Σ L = (W, −L, −L s , V).
We distinguish two cases: 1. The right amount of data is available-in this case, one can exactly recover the original system; 2. A redundant amount of data is provided-in this case, one can approximate the original system.
In the first case, the matrix pencil λL−L s is regular. Then, the transfer function of the Loewner model isĤ(s) = −W(sL−L s ) −1 V. This situation is encountered when the number of measurements is just enough to recover the original system (as explained in [34], one would require 2n measurements to recover a system of order n).
In the second case, the pencil λL − L s is singular. Then, the definition of the transfer function is adapted as in [2], i.e.,Ĥ(s) = −W(sL − L s ) # V .
In real-world applications, the reason behind the singularity of the matrix pencil λL − L s is that data contain, in most cases, redundant or inaccurate measurements.
We focus on the second case, which is more relevant in practice. In this case, since the Loewner pencil is singular, we can use the singular value decomposition (SVD) to compress the data. The truncation order r can be chosen by inspecting the singular value decay of the data. Deciding on a specific value of r is usually done by making a trade-off between reduced-order model dimension and accuracy of fit. Consider the (short) SVDs of augmented Loewner matrices where The singular values are nonnegative real scalars that give information on the numerical rank of matrices and are stored on the main diagonal of matrices S 1 and S 2 as N ]) and The matrices Y, X ∈ C N×r are obtained by selecting the first r columns of the matrices Y 1 and X 2 . The reduced Loewner system is constructed by projecting the matrices L, L s , V, W with Y * and X to the left and respectively to the right aŝ The projected Loewner model is given byΣ L = (Ŵ, −L, −L s ,V). By means of compression, the pencil λL −L s does not have any infinite eigenvalues. Hence, Σ L is a strictly proper system (its transfer function does not have a polynomial part). Hence, we can say that, in the modeling stage, the behavior at infinity of the underlying system is generically lost. We will address this issue in the next sections.
One can construct transformation matrices U ∞ and V ∞ as described in Section 3 that are used to block diagonalize the pencil λL − L s . The recovered polynomial coefficients can be computed similarly as introduced in Section 2.2.2, by replacing the inversion operator with the Moore-Penrose inversion operator, as follows

The main procedure
In the case of redundant data, significantly large amount of measurements are used in the modeling step. In this case, the matrix pencil λL − L s is singular. Instead of compressing the raw Loewner model as described in (42), we can alternatively first perform a post-processing procedure to the data in order to extract the required information.
The main procedure which adapts the Loewner framework to descriptor systems by preserving the DAE structure is given in Algorithm 1.
the Loewner matrices corresponding to the data in (45) constructed as in (39): Next, setting s = λ j and multiplying equation (48) with the j th unit vector e j on the right, we write: A similar reasoning is used to show that the equality in (46) also holds for the left interpolation points μ i .

Remark 3 A trust interval I (as used in
Step 3 of Algorithm 1) is an interval such that for any k 1 , k 2 ∈ [a, b], the polynomial coefficient estimates corresponding to k 1 , i.e., p (k 1 ) , and the ones corresponding to k 2 , i.e., p (k 2 ) have similar values, i.e., where ζ > 0 is a small tolerance value. As an example, in Fig. 7, one can choose the confidence interval to be I = [25,55], based on the values of the estimated coefficients.

Remark 4 In
Step 4, k is the estimated value denoting the dimension of the subsystem corresponding to the polynomial part. It is to be noted that this subsystem is not necessarily in a minimal realization format. In practice, e.g., the MNA example in Section 5.3.1, this value is chosen to be k = 17, although a minimal realization would only be two dimensional.
Remark 5 Descriptor systems appear in practical applications mostly with index ν ∈ {1, 2, 3} (the polynomial part is up to a quadratic term). Hence, in such cases, the value h ∈ N would be less or equal than 2. In all numerical examples in Section 5, we displayed only the first three estimated coefficients, i.e., p 0 , p 1 , and p 2 . For the MNA example, by choosing a tolerance value of ρ = 10 −10 , we included in the calculations only the first 2 polynomial coefficients; hence, h = 1 was considered (since p 2 < ρ).
Remark 6 The Loewner system Σ L corresponding to the post-processed data in (45) is strictly proper provided that the polynomial coefficients are perfectly estimated (exactly the same as the true coefficients corresponding to the original system Σ).

Remark 7
Provided that the value h is accurately chosen, it is hence assumed that the polynomial coefficients are given by p 0 , p 1 , . . . p h and p i = 0, ∀i > h. A minimal linear descriptor realization Σ poly = (C poly , E poly , A poly , B poly ) of order h + 1 corresponding to the polynomial part H poly (s) is constructed as where J h+1 is a Jordan block of dimension (h + 1) × (h + 1) with 0 eigenvalues, and I h+1 is the identity matrix of the same dimension.

Experimental results
In this section, we provide numerical results for various descriptor systems by implementing the algorithms from previous sections. We consider two simple low-dimensional examples as well as three larger scale benchmarks examples.
Next, choose the same interpolation points but ordered differently, as μ T = − 1, − 2, − 3, − 4 , and λ = 1, 2, 3, 4 . Then, put together the Loewner matrices and restart the procedure for the new matrices E 1 and A 1 , where The fact that the singularity corresponding to γ 1 = 2.3912 · 10 −1 is located on the first position for descending ordering (and, respectively on the last for the ascending ordering) influences the accurate recovery of the coefficients. Note that the results in Table 4 are correct only for k = 2, while in the previous experiment, they were correct for both k ∈ {2, 3} (see Table 2).
This shows, in particular, the sensitivity of the Loewner framework for different data sets in addition to the numerical sensitivity of the splitting scheme. That is why, for the large-scale examples, we will use so-called trust intervals for the k value. This is done in order to estimate the polynomial coefficients (for different values of k).
Finally, we assume that the sample values of the transfer function H(s) = 2s + 3 + 4 s−5 , evaluated at μ T ∪ λ, are affected by additive noise. We assume the noise level to be 10 −2 and the noisy values to follow a standard normal distribution. We repeat the experiment described above and compute the estimated polynomial coefficients for k = 2 (the value that should provide correct estimates): p 0 = 2.0623, p 1 = 1.8028, p 2 = − 0.0428, p 3 = − 0.0042, and p 4 = 0.0002. Hence, we conclude that the polynomial part extraction is indeed sensitive to noise.  The transfer function can be explicitly computed as follows:

An order n = 7 example
Hence, identify the coefficients of the polynomial part as p 0 = 3 2 , p 1 = 3 4 , and p k = 0 for k > 1. Note that the pencil λE − A has three eigenvalues at ∞ and four finite eigenvalues given by the two complex conjugate pairs: − 0.9712 ± 0.8138i and − 0.2787 ± 0.4834i.
Using the Jordan chain method in Section 3.1.1 for k = 3, the transformation matrices U ∞ and V ∞ can be explicitly obtained as We can obtain the coefficients of polynomial part as follows: Additionally, the ADTF1 method from Section 3.1.2 is applied to the original matrices in (50). The results can be found in Table 5. For k ∈ {0, 1, 2, 3}, we perfectly recover the true polynomial coefficients. For k ∈ {4, 5, 6}, the recovered values are not accurate.
Finally, in what follows, we apply this framework for data-driven splitting of the polynomial part. We show that the polynomial coefficients can be recovered from measurements of the transfer function. For example, choose 11 left and 11 right sample points on the real axis, collected in the vectors μ, λ T ∈ R 11 , μ T = − 1 2 11, 10, . . . , 1 , λ = 1 2 0, 1, . . . , 10 .
By means of the classical Loewner framework in [6], construct an 11th order system from measurements, characterized by a singular pencil (A, E). Then apply the proposed splitting technique to separate this system into two subsystems. We denote with k the dimension of the subsystem corresponding to the infinite eigenvalues.  We vary this value from 2 to 10. The results are collected in Table 6. Note that for k ∈ {2, 3, . . . , 7}, we perfectly recover the coefficients.

Descriptor system benchmarks
In this section, we apply the proposed method in Algorithm 1 to the following benchmarks problems: modified nodal analysis model 1 (MNA-1) [18], semi-discretized Oseen equations (Oseen) [26], constrained damped mass-spring system (CDMS system) [35]. The dimension and the index of these systems are collected in Table 7.
In what follows, we apply the proposed procedure in Section 3.2.1 for singular pencils.

MNA1 model
Modified nodal analysis (MNA) is a procedure used in the field of electrical engineering that determines the circuit's node voltages and branch currents. It was developed as a means to overcome the difficulty of representing voltage-defined components in classical nodal analysis (e.g., voltage-controlled sources).
As described in [18], voltage sources are connected to the network ports to obtain the admittance matrix of a multiport. The MNA equations are written as Eẋ n = Ax n + Bu p , i p = Cx n . The i p and u p vectors denote the port currents and voltages, respectively, while v and i are the MNA variables corresponding to the node voltages, inductor and voltage source currents, respectively. The conductance matrix A and the susceptance matrix E are given where N, L, and H are the matrices containing the resistors, capacitors, and inductors values, respectively. Also, G consists of 1, − 1, and 0 entries only. In this framework, if the circuit contains passive linear elements only, it follows that the E matrix is symmetric and positive semidefinite, while B = C T .
In this section, we analyze the MNA-1 system of dimension n = 578 from [18]. All system matrices of the original system are sparse and the number of inputs and outputs is m = p = 9. Next, consider instead a SISO-modified system (Ȃ,B,C, 0,Ȇ) for which the system matrices are written in terms of the original ones, asȂ = 10 −4 A,Ȇ = 10 9 E,B = Be 2 ,C = e T 2 C. Basically, we do a scaling of the original pencil λE − A and consider only the second input and output. We can directly compute the coefficients of the polynomial part of the transfer function as p true 0 = 5.5047 · 10 6 , p true 1 = 2.3010 · 10 3 and p true j = 0, j > 1. Hence, the index of the system is ν = 2.
First, depict the frequency response of the system in Fig. 1, i.e., the magnitude of the transfer function H(s) = C(sE − A) −1 B evaluated on the imaginary axis, in the interval [10 −2 , 10 1 ]i.
By inspecting the frequency response, we observe that the dominant oscillations mainly occur in the frequency range [10 −2 , 10 0 ]i. The next step is to take measurements in a high frequency range and extract the polynomial coefficients from this data.
Consequently, take N = 40 measurements in the frequency range [10 2 , 10 6 ]i and construct a real-valued 40th order Loewner model (by including the complex conjugate values in the data set, as performed in [6], page 360). By varying the value of the index k from 1 to 40, compute the values of p j , j ∈ {0, 1, 2} for each k. Afterwards, choose the trust interval [10,25] in which the values of the estimated polynomial coefficients have similar values for all k's, as can be observed in Fig. 2.
The relative deviation between the true and estimated coefficients is depicted in Fig. 3 for p 0 and p 1 only. The reason for this is because p true 2 = 0. We select the estimated polynomial coefficient values corresponding to the center of the interval, i.e., to k = 17.
Next, proceed with the remaining steps of Algorithm 1. In the next modeling step, we are using a different set of measurements. In order to match the oscillatory behavior of the system in the lower frequency range, select 600 logarithmically spaced points in the interval [10 −2 , 10 2 ]i. Again, construct a real Loewner model of Fig. 1 The frequency response of the descriptor system Adv Comput Math (2020) 46: 3 Fig. 2 The estimated polynomial coefficients of the transfer function for k ∈ [10,25] dimension 600. The decay of the normalized singular values corresponding to the augmented Loewner matrices is presented in Fig. 4.
Finally, we subtract the estimated polynomial part from this Loewner model, compress the model to order r = 28 (corresponding to σ 28 ≈ 10 −2 singular value), and then reattach the polynomial part given by the estimated polynomial coefficients-a minimal realization is considered. The dimension of the resulting system after reattaching the polynomial part is hencer = 31 (p 2 is also taken into consideration). This is true for the other experiments in this section (in Sections 5.3.2 and 5.3.3).
The frequency responses of the original large-scale model and of the two reducedorder Loewner models (with and without performing the post-processing step) are presented in Fig. 5. Notice that there is a clear mismatch between the original response and the one produced with the classical Loewner approach (which is visible starting at the frequency value 10 3 ).
Additionally, the relative approximation error is depicted in Fig. 5. This was computed using the formula H(iω) −Ĥ(iω) / H(iω) , for ω ∈ [10 −2 , 10 6 ]. For the two other numerical examples treated in this section, the relative approximation error in the frequency domain is computed using the same formula (but for different frequency ranges). Note that the advantage of using the new method consists in drastically lowering the deviation in the high frequency range.

Constrained damped mass-spring system
For the second benchmark example, we consider the holonomically constrained damped mass-spring system from [35]. Here, the ith mass of weight m i is connected to the (i + 1)st mass by a spring and a damper with constants k i and d i , respectively, and also to the ground by a spring and a damper with constants ξ i and δ i , respectively. Additionally, the first mass is connected to the last one by a rigid bar and it is influenced by the control input signal u(t). The vibration of this system is described by a descriptor system as in equation (34) from [35].
We consider the same values for the ξ i , δ i , k i , and d i parameters as in [35] and choose g = 500 for the dimension of the variables. Hence, by grouping together the position and the velocity vectors into the new x variable, we end up with a descriptor system of order n = 1001, described by Eẋ = Ax + Bu, y = Cx, where E, A ∈ R 1001×1001 and B, C T ∈ R 1001 . As opposed to the original system in [35], where it was considered that B = e g+1 ∈ R g and C = [e 1 e 2 e g−1 ] T ∈ R 3×g for g 1 (hence a system with one input and three outputs), we choose B = C T = [1 1 . . . 1] T ∈ R 1001 and analyze a SISO descriptor system instead.
For the original 1001th order system constructed as described above, we can directly identify the following coefficients of the polynomial part Hence, it follows that the original system is of index ν = 3. Take measurements corresponding to N = 100 logarithmically spaced sampling points in the frequency range [10 −2 , 10 2 ]i. Then, construct a real-valued 100th order Loewner model. The frequency response (restricted to the above sampling interval) of the system is presented in Fig. 6 together with the singular value decay of the augmented Loewner matrices . Again, we vary the value of the index k from 1 to 100 and compute the values of polynomial coefficients p j , j ∈ {0, 1, 2} for each k. The results are gathered in Fig. 7. By inspecting this figure, we observe that one can choose, for example, the trust interval [25,55]. Fig. 7 The estimated polynomial coefficients of the transfer function for k ∈ [1,100] Next, Fig. 8 shows the deviation between the true polynomial coefficients and the estimated ones for values of k restricted to the chosen trust interval. Next, we will apply Algorithm 1 for two such values, corresponding to k 1 = 27 and k 2 = 40. Denote with p ( ) j the estimated polynomial coefficients corresponding to index k for ∈ {1, 2} and j ∈ {0, 1, 2}. The absolute value of the deviation is computed as In what follows, we construct reduced-order models of order r = 8 (which corresponds to the singular value σ 8 ≈ 10 −5 ) For each of the values k 1 and k 2 , we compute the modified Loewner reduced model of order r = 8 via the procedure in Algorithm 1 as well as the reduced model using the classical Loewner model (of the same order). We also calculate the relative approximation error between the frequency response of the original large-scale system and that of the reduced models. The results are depicted in Fig. 9.
As expected, the approximation error is highly dependent on how well the polynomial coefficients are approximated. As it can be observed in the left side of Fig. 9, the approximation error is significantly larger than that in the right side of the figure (when applying Loewner ∞ ).

Semi-discretized Oseen equations model
In this section, we apply the new proposed method to a descriptor system obtained from semi-discretization of the Oseen equations, as performed and thoroughly explained in [26]. The Oseen equations are related to the linearized Navier-Stokes equations. The study of these equations is of importance in many flow systems where one is interested in the behavior of the linearized flow around a steady state. As presented in [26], after performing a spatial discretization, the dynamics is described in generalized state space by the equations where v ∈ R n v , p ∈ R n p are the states, g ∈ R n g are the inputs, and y ∈ R n y are the outputs, and E 11 is a symmetric positive definite matrix. The next step is to rewrite the system in (55) as a descriptor system in (1) of dimension n v + n p , by employing the following notation and . Hence, rewrite (55) as In the following, perform a discretization scheme corresponding to the finite element grid in the upper part of Fig. 12 (as in [26]). This corresponds to the following dimension values of the variables in (55): n v = 1352, n p = 205, n g = 6 and n y = 2. Hence, the original descriptor model constructed as in [26] has dimension n = 1557, 6 inputs, 2 outputs, and is of Hessenberg index ν = 2. The coefficients of the In the numerical simulations, consider only the first input and the second output. Hence, we analyze an index 2 descriptor system with original polynomial coefficients p true 0 = 16.6241 and p true 1 = 7.9211. Proceed to choosing data as transfer function measurements corresponding to N = 40 logarithmically spaced sampling points in the frequency range [10 1 , 10 6 ]i. Afterwards, construct a real-valued 40th order Loewner model. The frequency response (restricted to the above sampling interval) of the system is presented in Fig. 10 together with the singular value decay of the augmented Loewner matrices.
Next, vary the value of the index k from 1 to 40 and compute the values of polynomial coefficients p j , j ∈ {0, 1, 2} for each k. The results are restricted to the trust interval [6,26] and are depicted in Fig. 11.

Fig. 12
Frequency response comparison and approximation error reduced-order models of order r = 12 (which corresponds to the singular value σ 12 ≈ 10 −12 ). We compare the frequency response of the original system with that of the two reduced models of order 12 (computed via the classical Loewner approach and the new modified approach) in Fig. 12. Again, notice that there is a clear mismatch between the original response and the one of the classical Loewner (which is visible starting at around 10 3 i). The response computed with the new method follows the original curve. Moreover, in the right side of Fig. 12, we present the relative approximation error between the original and reduced-order models. Clearly, the new method performs better in the high frequency range.
Finally, we perform a time domain simulation and compute the observed output y(t) by means of an explicit forward Euler-type scheme. The control input is given by u(t) = 0.2 sin(10πt) + 0.3 sin(4πt)e −t , while the time horizon is chosen to be [0, 4]s for this experiment. We consider zero initial conditions. In the upper part of Fig. 13, we depict the observed output corresponding to the original large-scale system as well the output of the reduced-order system compute via the new proposed method. Notice that the two curves are practically indistinguishable. Additionally, we compute the absolute error between the two observed plots. The result is presented in the lower part of Fig. 13. We observe that the deviation between the two curves is indeed small. The purpose of this experiment was to test the accuracy of the proposed method outside the frequency domain for applications that are primarily considered in time domain (e.g., differential equations such as the Oseen equations).

Conclusion
The philosophy behind the classical Loewner framework is that by means of directly compressing the surrogate data model (put together solely from the available measurements), we obtain a reduced-order model that matches the interpolation conditions. In this work, we have addressed an open issue commonly encountered in this framework: the fact that the behavior at infinity of the original system is typically not matched by the reduced system.
We propose an algorithm that is based on existing methods of splitting the raw data Loewner model into two subsystems: one corresponding to the strictly proper part and one corresponding to the polynomial part. The splitting can be done for both regular and singular pairs of Loewner matrices.
Afterwards, we estimated the polynomial coefficients of the transfer function by setting up a trust interval (in which the coefficients have more or less constant values). We have shown that the choice of the coefficients plays a crucial role in the quality of approximation.
Three numerical test cases were chosen amongst well-established benchmark example from the literature. They show the practical applicability of the new proposed method to medium-and large-scale problems.
Finally, there are some open issues and possibly interesting research directions. These include investigating other splitting techniques (for singular pencils) as well as incorporating the case of noisy data in the process (and studying/quantifying the sensitivity to noise of the estimated polynomial coefficients). Finally, the proposed procedure could be extended to DAE systems with nonlinearities (e.g., bilinear control systems) .

Funding Information Open access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.