Practical challenges in data-driven interpolation: dealing with noise, enforcing stability, and computing realizations

In this contribution, we propose a detailed study of interpolation-based data-driven methods that are of relevance in the model reduction and also in the systems and control communities. The data are given by samples of the transfer function of the underlying (unknown) model, i.e., we analyze frequency-response data. We also propose novel approaches that combine some of the main attributes of the established methods, for addressing particular issues. This includes placing poles and hence, enforcing stability of reduced-order models, robustness to noisy or perturbed data, and switching from different rational function representations. We mention here the classical state-space format and also various barycentric representations of the fitted rational interpolants. We show that the newly-developed approaches yield, in some cases, superior numerical results, when comparing to the established methods. The numerical results include a thorough analysis of various aspects related to approximation errors, choice of interpolation points, or placing dominant poles, which are tested on some benchmark models and data-sets.


Introduction
Approximation of large-scale dynamical systems is pivotal for serving the scopes of efficient simulation and designing control laws in real-time.The technique for reducing the complexity of a system is known as model order reduction (MOR) [1,5,13,14].There exist a number of methodologies for reducing large-scale models, and each method is tailored to some specific applications (mostly, but not restricted to mechanical and electrical engineering) and to achieving certain goals (stability, passivity or structure preservation), on top of the complexity reduction part.Data-driven MOR approaches are of particular importance when access to high-fidelity models is not explicitly granted.This means that a statespace formulation with access to internal variables is not available, yet input/output data are.Such methods circumvent the need to access an exact description of the original model and are applicable whenever classical projection-based MOR is not.Here, we mention system and control methodologies that are based on interpolation or least-square fit of data (i.e., frequency response measurements), such as vector fitting [31], the Loewner framework [41], or the AAA algorithm [42].Methods that use time-domain data are also of interest, including the ones that require input-output data together with those which use snapshot data (access to the state evolution), such as the classical ones in [36,58,59], followed by [8], [47] or [54].
We focus on interpolation-based or so-called moment matching (MM) methods which have emerged, were developed, and improved continuously in the last decades.The backbone of such methods represent rational Krylov-type approaches together with the Sylvester matrix equation interpretation [1,11,20].Apart from being computationally efficient and easy to implement, MM approaches have another advantage: they do not (necessarily) require access to a full-state realization of the original dynamical system.Hence, they can be viewed as data-driven methods.Here data are given by the moments of the system, i.e., samples of the underlying transfer function of the system (and of its derivatives) evaluated in a particular frequency range; for more details, we refer the readers to [5,8,41] and to Chapter 3 in [13].The notion of a moment with respect to systems and control theory is related to the unique solution of a Sylvester matrix equation [20].
The purpose of this note is twofold; first, we intend to review and to connect three important system theoretical model reduction methods based on interpolation that were introduced in the last 15 years: • The Loewner framework (LF) by Mayo and Antoulas from 2007 in [41]; • The Astolfi framework (AF) from 2010 in [8]; • The Adaptive Antoulas Anderson (AAA) algorithm by Nakatsukasa, Seté and Trefethen from 2018 in [42].
Together, these three approaches were cited multiple times in various research publications, being arguably quite popular methods.However, until now, not too many connections between them were provided, neither in the automatic control, nor in the model reduction, or numerical analysis communities.Together with the vector fitting algorithm (VF) in [31] (which is not based on interpolation, and is hence a purely optimization approach based on least-squares fitting), these methods represent arguably the most prolific rational approximation schemes developed in the system and control community.However, VF is not the object of this study since it is not based on interpolation.
The other scope of this note is to propose a new method that is based on the three methods enumerated above, and that addresses some of the shortcomings and challenges associated with them.Basically, the idea is to combine the attributes of each method, by following the steps below.
• We make use of the order-revealing property of the LF (encoded by the rank of augmented Loewner matrices); additionally, the selection of interpolation points is done via a Loewner-CUR technique proposed in [38].
• We utilize the elegant state-space parameterization of the LTI system proposed by the AF (after imposing k interpolation conditions); this is the backbone of the methods (we also show the connection between state-space forms and barycentric forms).
• We use either the fitting step from AAA (that chooses free parameters to fit the un-interpolated data in a least square sense) or we impose pole placing (dominant poles are selected from those of the Loewner model); in both cases, a linear system of equations needs to be solved.
In what follows, we consider a multiple-input multiple-output (MIMO) linear time-invariant (LTI) system Σ L of dimension n described by the following system of differential equations: with x(t) ∈ R n as the state variable, u(t) ∈ R m as the control inputs, and y(t) ∈ R p as the observed outputs.Here, we have that A ∈ C n×n , B ∈ C n×m and C ∈ C p×n .The transfer (matrix) function of the LTI system is given by H(s) ∈ C p×n , with s ∈ C, as It is to be noted that, for m = p = 1, the system becomes single-input single-output (SISO).We will sometime switch between MIMO and SISO formats while presenting the methods covered in this note, since the latter allows a more easy exposition for some of the results shown here. Let , where σ(A) denotes the spectrum of matrix A ∈ C n×n , i.e., the set of its eigenvalues.The j-moment of system Σ L at s i is given by η j (s i ) = (−1) j j! d j ds j H(s) s=si , for any integer j 1.The 0-moment is obtained by sampling the transfer function H(s) in (2) at s i , i.e., η 0 = H(s i ).In this contribution, we restrict the analysis to matching 0-moments, i.e., samples of the transfer function H(s), and not of its derivatives.However, all methodologies shown here can be expected to cope with this as well.Moreover, in practice, inferring 0-moments from time-domain data is usually a more straight-forward task; this is performed by exciting the system with harmonic inputs, and by applying spectral transformations to the outputs.Additionally, the inference of derivative values (of the transfer function) is typically susceptible to perturbations and it is more challenging to attain, from a numerical point of view.
The paper is structured in the following way; after the introduction session sets up the stage, we propose a survey of three established interpolation-based methods in Section 2.Then, the proposed methodologies are developed in Section 3, with emphasis on the one-step approach that combines optimal selection of interpolation points (chosen using CUR-DEIM) with LS fit on the rest of the data, and also the pole placement method in barycentric form that enforces dominant poles from the Loewner data-driven model.Then, Section 4 illustrates the numerical aspects of applying the methods discussed/proposed in the previous two sections to a variety of test cases (various models and data sets).Finally, Section 5 presents the conclusions and the outlook into future research.

A survey of established methods
In this section we discuss three established data-driven methods for rational approximation (AF, LF, and AAA, as mentioned in the previous section).The data are samples of the transfer function corresponding to the underlying dynamical system, measured on a particular frequency grid.In what follows, we mention some state-of-the-art methodologies used to measure such data, i.e., frequency response data.Typically, such measurements can be produced in practice from experiments conducted in scientific laboratories using carefully calibrated machines, called spectrum analyzers (SAs).In this category we mention swept-tuned spectrum analyzers, scalar network analyzers (SNAs), and vector network analyzers (VNAs).
The SNA is an instrument that measures microwave signals by converting them to a DC voltage using a diode detector.In a VNA, information regarding both the magnitude and the phase of a microwave signal is extracted.While there are different ways to perform such measurements, the method employed by commercial products (such as the Anritsu series described in [18]) of VNAs is to down-convert the signal to a lower intermediate frequency in a process called harmonic sampling.This signal can then be measured directly by a tuned receiver.Compared to the SNA, the VNA is a more powerful analyzer tool.The major difference is that the VNA can also measure the phase, and not only the amplitude.With this property enforced, then so-called scattering parameters (or S-parameters) can be processed.These can be used for identifying forward and reverse transmission and reflection characteristics.More details can be found in [18].
The harmonic balance method (or HBM) [43], is an established methodology in the field of electromagnetics.The HBM is used in many (if not most) commercial radio-frequency (RF) simulation tools.This is due to the fact that it has certain advantages over other common methods, namely modified nodal analysis (MNA), which makes it more appropriate to use for stiff problems and circuits containing transmission lines, nonlinearities and dispersive effects.More details can be found in the survey paper [49].
2.1.The one-sided moment-matching approach in [8] The framework introduced by Astolfi in [8] (referred to as AF throughout the paper) deals with the problem of model reduction by moment matching.Although classically interpreted as a problem of interpolation of points in the complex plane, it has instead been recast as a problem of interpolation of steady-state responses.In the following we briefly review its application to linear systems.It is to be noted that the AF was steadily extended and applied to different scenarios (including nonlinear dynamical systems, pole-zero placement, and least-squares fit) [33-35, 45, 53, 55].
The moments of a linear system can be characterized in terms of the solution of Sylvester equations.By using this observation, it has been shown that the moments are in one-to-one relation with the steady-state output response of the interconnection between a signal generator and the original linear system.
In what follows, for simplicity of exposition, it is assumed that Σ L is a minimal system (both fully controllable and fully observable).For exact definitions on minimality, controllability, or observability of LTI systems, we refer the reader to [1].
Let k n and S ∈ C k×k a non-derogatory matrix (for which the characteristic and minimal polynomials coincide) with σ(S) ∩ σ(A) = ∅ and R ∈ C 1×k so that (S, R) is observable.Consider the signal generator system Σ sg described by the equations Then, the explicit solution of (3) can be written as ω(t) = e St ω(0).Hence, the control input is written as u(t) = Re St ω(0).In addition, the eigenvalues of S are called interpolation points.For a linear system Σ L , and interpolation points s i ∈ C \ σ(A), for i = 1, . . ., k, consider a nonderogatory matrix S ∈ R k×k .It follows that there exists a one-to-one relation between the moments of the system Σ L and 1. the matrix CΠ, where Π is the (unique) solution of the Sylvester equation AΠ + BR = ΠS, for any row vector R ∈ R 1×k so that (R, S) is observable; 2. the steady-state response of the output y of the interconnection of system Σ L and the system Σ sg , for any R and ω(0) such that the triplet (R, S, ω(0)) is minimal.
More precisely, let ∆ ∈ R k be a column vector containing k free parameters (denoted here by δ 1 , δ 2 , . . ., δ k , with δ i = 0, 1 ≤ i ≤ k).Then, as stated in [8], the family of linear time-invariant systems that interpolates the moments of system Σ L at the eigenvalues of matrix S, is given by Σ∆ : where the matrices S and R are as before and Π is the unique solution of the Sylvester equation AΠ + BR = ΠS.Additionally, the condition σ(S) ∩ σ(S − ∆R) = ∅ needs to be enforced.It is to be noted that the free parameters explicitly enter the vector B = ∆, but also the matrix Â, as Â = S − ∆R.Finally, Ĉ = CΠ has fixed entries.
The Sylvester matrix equation for the reduced-order system is written as Â Π + BR = ΠS.This can be explained by the fact that the reduced-order system matches the prescribed moments of the original system, hence the same format of the two equations.Without loss of generality, one can consider that the matrix Π is the identity matrix, i.e.Π = I k (this can be achieved by applying similarity transformations).By replacing this value into the reduced-dimension Sylvester matrix equation above, the formula Â = S − ∆R directly follows.
Afterwards, the free parameters collected in the vector ∆ can be chosen in order to enforce or impose additional conditions as mentioned in [8], such as: matching with imposing additional k interpolation conditions, matching with prescribed eigenvalues, matching with prescribed relative degree, matching with prescribed zeros, matching with a passivity constraint, matching with L 2 -gain, or matching with a compartmental constraint.
An important aspect of the AF is the characterization of all, i.e., infinitely many families of reducedorder models that satisfy k prescribed interpolation conditions.This is done by explicitly computing such parameterized models, for which the free parameters are the variables entering in the vector ∆.The main parameterization developed here will be used as a "backbone" of the methods developed in Section 3.
As stated in the original paper, the main advantage of the AF (characterization of moments in terms of steady-state responses) is that it allows the definition of moments for systems which do not admit a clear/immediate representation in terms of transfer function(s).Hence, the author provides as examples the case of linear time-varying systems, and the case of nonlinear systems.Moreover, it is stated in [8] that one disadvantage of the framework is that it requires the existence of steadystate responses.Consequently, the original system has to be (exponentially) stable.However, in most practical applications, this is a realistic requirement.[41] In this section we present a short summary of the Loewner framework (LF), as introduced in [41].It is to mentioned that LF has its roots in the earlier work of [4], and that LF can be considered to be a double-sided moment-matching approach (as opposed to AF, which is one-sided).

The Loewner framework in
For a tutorial paper on LF for LTI systems, we refer the reader to [7], and for a recent extension that uses time-domain data, we refer the reader to [48].The Loewner framework has been recently extended to certain classes of nonlinear systems, such as bilinear systems in [6], and quadratic-bilinear (QB) systems in [24], but also to linear parameter-varying systems in [28].Additionally, issues such as stability preservation or enforcement, or passivity preservation, were tackled in the LF in [23,29], for the former, and in [2,12], for the latter.
The interpolation problem is formulated as shown below (for convenience of exposition, we show here only the SISO formulation).We are given data nodes and data points in the set D, partitioned into two disjoint subsets D L and D R , with and we seek to find a rational function Ĥ(s), such that the following interpolation conditions hold: The Loewner matrix L ∈ C q×k and the shifted Loewner matrix L s ∈ C q×k play an important role in the LF, and are given by while the data vectors V ∈ C q , W T ∈ C k are given by Moreover, the following Sylvester matrix equations ([1, Ch. 6]) are satisfied by the Loewner and shifted Loewner matrices (here, where M = diag(µ 1 , . . ., µ q ) and Λ = diag(λ 1 , . . ., λ k ) are diagonal matrices.The following relation holds true The unprocessed Loewner surrogate model, provided that k = q, is composed of the matrices and if the pencil (L, L s ) is regular, then the function Ĥ(s) satisfying the interpolation conditions in (6) can be explicitly computed in terms of the matrices in (11), as Ĥ(s In practical applications (when processing a fairly large number of measurements), the pencil (L s , L) is often singular.Hence, a post-processing step is required for the Loewner model in (11).In such cases, one needs to perform a singular value decomposition (SVD) of augmented Loewner matrices, to extract the dominant features and remove inherent redundancies in the data.By doing so, projection matrices X, Y ∈ C k×r are obtained, as left, and respectively, right truncated singular vector matrices: where The truncation index r can be chosen as the numerical rank (based on a tolerance value τ > 0) or as the exact rank of the Loewner pencil (in exact arithmetic), depending on the application and data size.More details can be found in [7].
The system matrices corresponding to a projected Loewner model of dimension r can be computed as follows: We note that MIMO extensions of the LF were already proposed in the original contribution [41].There, a tangential interpolation framework is considered.Instead of imposing interpolation of full p × m blocks, the authors prefer to interpolate the original transfer matrix function samples along certain vectors (or tangential directions).We also note that a first attempt of re-interpreting the LF in [41] as a one-sided method was made in [25].In the latter, the main difference to the classical work in [4] was that a compression of the left (un-interpolated) data set was enforced.However, in [25], it was still unclear how to split the data, i.e., what the right data set should be (where interpolation is enforced).Finally, it is to be noted that the choice of interpolation points is crucial in the LF.An exhaustive study of different choices was proposed in [37], while a greedy strategy was proposed in [17], for scenarios in which limited experimental data are available.

The AAA algorithm in [42]
The AAA algorithm introduced in [42] is an adaptive and iterative extension of the interpolation method based on Loewner matrices, originally proposed in [4].The main steps are as follows 1. Express the fitted rational approximants in a barycentric representation, which represents a numerically stable way of expressing rational functions [15].
Algorithm 1 The AAA algorithm.
Require: A (discrete) set of sample points Γ ⊂ C with N points, function f (or the evaluations of f on the set Γ, i.e., the sample values), and an error tolerance > 0. Ensure: A rational approximant r n (s) of order (n, n) displayed in a barycentric form.
10: j ← j + 1. 11: end while 2. Select the next interpolation (support) points via a greedy scheme; basically, interpolation is enforced at the point where the (absolute or relative) error at the previous step is maximal.
3. Compute the other variables (the so-called barycentric weights) in order to enforce least squares approximation on the non-interpolated data.
In recent years, the AAA algorithm has proven to be an accurate, fast, and reliable rational approximation tool with a fairly large range of applications.Here, we will mention only a few: nonlinear eigenvalue problems [39], MOR of parameterized linear dynamical systems [16], MOR of linear systems with quadratic outputs [26], rational approximation of periodic functions [10], representation of conformal maps [22], rational approximation of matrix-valued functions [27], or signal processing with trigonometric rational functions [60].The procedure is sketched in Algorithm 1.
It is to be mentioned that a modified version of AAA that enforces real-valued and strictlyproper rational appoximants was recently proposed in [30].There, the format of the function in (13) was modified by inserting a 1 into the denominator, as follows rj (s Consequently, the equation in (14) becomes L (j) ω (j−1) = −f (j−1) , where the vector f (j−1) ∈ C j is given by f . It is to be noted that rj (s) in ( 15) is theoretically a rational approximant of order (j − 1, j), if we do not take into account pole/zero cancellations or any other zero cancellations of coefficients in the numerator or denominator.

Skeleton of the main methods
Similar to the methods reviewed in Section 2 we want to find an LTI system with a transfer function of the structure (1) that interpolates data provided as measurements H (s i ) , i = 1, . . ., k of the transfer function of the original system.We can directly put together an LTI parametrized model of dimension r = km, having km 2 degrees of freedom with transfer function with the underlying data concatenated to a matrix of weights Ŵi B = ŴH and Â ∈ C r×r formed from a diagonal matrix populated with the interpolation points λ i disturbed by B, such that Making use of the Woodbury matrix identity and denoting Λ s = sI km − Λ, the transfer function ( 16) can be rewritten as A complete derivation of ( 20) is given in Appendix A.1.
In the single-input single-output case (m = p = 1, hence r = k), the barycentric weights reduce to scalars and the matrices for a ROM of structure (16) are given by By inserting the formulae in ( 21) into (20), and using the notation h i := H (λ i ), leads to Hence, the transfer function of the model in ( 21) is given in barycentric representation by This can be performed analogously for a multi-input multi-output case (m = p, r = km).The first part of (20) becomes the second part Consequently, the transfer function in (20) has also a barycentric form given by The transfer function is defined by the choice of the interpolation points and of the weights.The interpolation points can be chosen as dominant parts of the available data or based on their location in the frequency spectrum.The weights can be computed such that the data which are not interpolated, are approximated in an optimal way.Alternatively, the weights can be chosen to enforce poles at specific locations.In the following, we show different strategies for both choices.

Automatic choice of interpolation points
The approximation quality of a surrogate model of the form ( 16) is greatly influenced by the choice of the interpolation points λ.This choice is not always obvious, so automatic strategies are frequently employed.The Loewner framework uses the SVD to identify dominant subsets of the available data to enforce interpolation on.Alternatively, the AAA algorithm uses a greedy scheme to minimize the error between surrogate and original data.Another approach, originally introduced by [37], makes use of the CUR decomposition to extract interpolation points from a relevant subset of the available data.
The CUR decomposition approximates a matrix A by a product of three low-rank matrices Ǎ = Č Ǔ Ř, where Č and Ř represent subsets of the columns respectively rows of A [40,56].In our case the three matrices are only a byproduct, we are more interested in the interpolation points λ and µ that are associated to the columns and rows extracted as Č and Ř.In combination with the skeleton for a realization described in Section 3.1, Algorithm 2 computes a surrogate model approximating a set of given transfer function data.We use the algorithm from [56] to compute the CUR decomposition and thus identify dominant parts of the original data set and their corresponding left and right interpolation points.Contrary to [37], we decompose the original Loewner matrix L rather than the augmented Loewner matrices L L s and L H L H s H . Using all interpolation points obtained from the CUR decomposition would introduce redundant data into the surrogate.Therefore we choose only a subset of the interpolation points: either only the left points, only the right points, or every other entry from a concatenated and sorted vector of left and right points.Together with the data associated to the chosen interpolation points they are used to populate a rectangular Loewner matrix.We now need to compute weights for barycentric interpolation as described in the following section.After having obtained the weights, a surrogate model ( 16) can be computed from ( 17)- (19).6: Compute the weights , where L † is the pseudo-inverse of L and 7: Compute Â, B, Ĉ with ( 17)-( 19).

Computing the barycentric weights 3.3.1. Least-squares approach
The matrix-valued weights Ŵi can be computed similarly to AAA [27] by solving the minimization problem This solution can, for example, be obtained from an optimization in least-squares sense.The weights for the SISO case are computed analogously.Here, the matrix-values weights and transfer function values reduce to scalars.

Pole placement
The next step would be to take advantage of the degrees of freedom in the vector B from (21), so that the ROM thus constructed has particular (stable) poles [21,35,46].These will be denoted with ζ 1 , ζ 2 , . . ., ζ k .The following derivations assume a SISO model.To enforce that this happens, we need to make sure that the matrix ζ j I k − Â loses rank for all 1 ≤ j ≤ k.In what follows, we show how to enforce this property in an elegant, straightforward way.Remember that the transfer function of the parameterized AF model is given by: Now, let's say we would like this transfer function to have k poles at the selected values ζ j 's.Clearly, the condition is D(ζ j ) = 0 and hence we need to enforce: where C ζ,λ is a Cauchy matrix defined by: (C ζ,λ ) i,j = 1 ζi−λj .Details on how to obtain the above expression by following the procedure in [3] are given in Appendix A.2.We note that placing poles is a difficult numerical problem which requires the inversion of a Cauchy matrix, which is highly ill-conditioned, by nature.

Require: Transfer function samples {H (s
and {s i } N i=1 using the Loewner framework (cf.Section 2.2).
Instead of doing this, we could solve C ζ,λ B = −1 r , without inverting the Cauchy matrix explicitly, i.e., by solving a linear systems of equations.Algorithm 3 summarizes this procedure in a data-driven context.The required underlying model is obtained from a set of transfer function evaluations by applying the Loewner framework.The method is illustrated for SISO systems, but can readily be extended to the MIMO case.

Automatic choice of poles and interpolation points
A reasonable choice of poles and interpolation points for Algorithm 3 is not always readily available, but the approximation of the surrogate is heavily influenced by this choice.In the following, we show an extension to Algorithm 3 which computes a surrogate model ( 16) without requiring sets of poles and interpolation points as input parameters.Algorithm 4 sketches the skeleton of such automatic algorithm.Similar to Algorithm 3 it employs the Loewner framework to obtain a realization of a surrogate interpolating the provided data.Subsequently, a generalized eigendecomposition of the Loewner realization of the original data is computed to find suitable locations for poles.From this is is possible to compute the dominance of all eigenvalues; for details, see, e.g.[51].The algorithm now chooses the k most dominant eigenvalues as poles to enforce in the surrogate.It should be noted that only eigenvalues with negative real parts should be considered, if the stability of the surrogate is important.The required interpolation points can now be chosen similar to Algorithm 2 by computing a CUR decomposition and using the interpolation points associated to the rows or columns of the decomposition as interpolation points for the new surrogate.
The approximation of dominant poles of the underlying model from data is less robust if the transfer function samples are disturbed by noise.This leads to a reduced approximation quality.For a better performance if applied to noisy data, Algorithm 4 can be modified as follows: To obtain poles which should be enforced, first choose manually the most prominent features in the transfer function, e.g.peaks, which should be approximated by the surrogate model.Now choose the eigenvalues which imaginary parts are closest to the frequencies, where the chosen features of the transfer function are located.The CUR decomposition also fails at extracting the most dominant rows and columns of the Loewner matrix if noisy data is assessed.Therefore another heuristic is employed to choose the interpolation points: Use the value s i which corresponds to the lowest amplitude of the transfer function between the locations of two enforced poles.This leads to reasonable approximations, especially for lightly damped systems.Other approaches include choosing simply the middle between the location of two poles or specifying an offset between pole and interpolation point location.

Numerical results
In the following, we demonstrate the methods discussed in Section 3 by applying them on three benchmark examples available from the MOR-Wiki 1 : Algorithm 4 Loewner framework with automatic pole placement (LFaPP).

Require: Transfer function samples {H (s
and {s i } N i=1 using the Loewner framework (cf.Section 2.2). 2: Compute the generalized eigenvalue decompositions AX = EXα and Y H A = αY H E for the matrices of right and left eigenvectors X, Y and the matrix of eigenvalues α = diag (α 1 , . . ., α n ).ISS This system models the structural response of the Russian Service Module of the International Space Station (ISS) [52].The model has n = 270 states, m = 3 inputs, and p = 3 outputs.The dataset used for the computations contains transfer function measurements at 400 logarithmically distributed points in the range 10 −1 , 10 2 • ı.The model is also part of the SLICOT benchmark collection [44].
Flexible aircraft This system models lift and drag along the flexible wing of an aircraft.The system matrices are not available, we only have access to a dataset of 420 transfer functions samples at linearly distributed frequencies between 0.1 and 42.0 Hz.The original dataset has one input (the gust disturbance) and 92 outputs.For the following experiments, we choose the 91st output which corresponds to the first flexible mode [50].The dataset is available from [57].

Sound transmission
This system models the sound transmission through a system of two brass plates with an air enclosure between them.The transfer function measures the sound pressure in an adjacent acoustic cavity.The geometry is based on [32]; the data-transfer function evaluations at 1000 linearly-distributed frequency values between 1 and 1000 Hz-is available from [9].
We note that no tangential interpolation (as described in [41]) is applied for the MIMO model.Instead, the Loewner matrices are constructed in a block-wise manner.The case of tangential interpolation, within the proposed approaches in this note, will be investigated in future works.
We enforce realness of all surrogate models (all matrices contain only real entries) by applying the transformation described in [7].For this, all data must be available in complex conjugate pairs.The required transformation matrix is given by with = k 2 and the real-valued quantities are obtained from For some of the experiments we add artificial noise to the measurements, in order to obtain perturbed data.The modified measurements are given by where Z i ∈ C is the ith sample drawn from a set of random numbers Z ∼ CN µ, σ 2 following a complex normal distribution with mean µ and standard deviation σ 2 .Here, the real and imaginary parts of Z are independent normally distributed variables [19].We assess the approximation error of the surrogate models with an approximated L ∞ norm, because many surrogates have unstable poles and hence, the H ∞ can not be computed.For a given reduced order r, the L ∞ error in the considered frequency range ω ∈ [ω min , ω max ] is approximated by Note that strategies to post-process surrogates to obtain stable models have been studied in [29].
The numerical experiments have been conducted on a laptop equipped with an AMD Ryzen™ 7 PRO 5850U and 12 GB RAM running Linux Mint 21 as operating system.All algorithms have been implemented and run with MATLAB R2021b Update 2 (9.11.0.1837725).

Code and data availability
The data that support the findings of this study are openly available in Zenodo at doi:10.5281/zenodo.7490158under the BSD-2-Clause license, authored by Quirin Aumann and Ion Victor Gosea.

Case of exact measurement data
In the following, we compare the performance of the new approach LS-Loewner to the following established strategies: • Loewner-SVD: Truncate Loewner matrices populated with the complete dataset to order r using an SVD [7].
• Loewner-CUR: Construct a purely interpolatory model of order r using all data points chosen by the CUR decomposition, similar to [37].
• Modified AAA: Apply the strictly-proper variant of AAA [29] to the complete dataset to compute a reduced-order model of size r.
We first consider the original MIMO ISS example and a SISO variant where we select the first input and output, respectively, from the MIMO system.To evaluate the overall performance of the different methods related to the size of a surrogate model, we compute the approximated L ∞ errors for models with orders 6 ≤ r ≤ 60.The approximation error versus the dimension of the respective surrogate model is depicted in Figure 1 for all four methods.
Since tangential interpolation was not employed here, the order of the MIMO surrogates rises by m for each additional interpolation point, i.e., r = km.This explains the lower accuracy of the MIMO surrogate.For the maximum reduced order r = 60, k = 20 interpolation points are considered.The errors of the SISO surrogates for r = 20, i.e., k = 20, is in a similar range as in the MIMO case.The SISO surrogates reach similar levels of approximation for all employed methods.In the MIMO case, Loewner-SVD performs best.This can be explained by the following observation: the other methods always consider the complete transfer function measurement H(λ i ) ∈ C p×m per interpolation point, while Loewner-SVD extracts only the r most dominant singular vectors for projection, regardless of to which interpolation point they belong to.In turn, the other methods also consider probably less important parts of the data as long as one input/output combination of the respective sample is relevant for approximation.It can also be noted that LS-Loewner and Loewner-CUR perform very similar.This was expected, as both methods rely on the same interpolation points.
All four methods are now employed to compute a surrogate model of size r = 108 to approximate the transfer function of the flexible aircraft model.The size of the surrogate model is determined by truncating all singular values τ < 1•10 −6 of an underlying Loewner matrix.
The transfer functions of all resulting models and their respective relative errors are given in Figure 2. Again, all methods succeed in computing a sufficiently accurate surrogate.However, the approximation quality of Loewner-CUR is noticeably worse than that of the other three methods.Given that both Loewner-CUR and LS-Loewner use the same interpolation points, the weights computed from the least squares problem show a better performance compared to the partitioning approach used in Loewner-CUR.

Perturbed measurement data
Analyzing measurement data perturbed by noise is a challenging task for interpolatory methods, such as the Loewner framework and the AAA algorithm (as pointed out in, e.g., [27]).In this experiment we investigate the effect of noise to the performance of the four methods described above and show, how enforcing poles and/or interpolation points can increase the approximation quality.In the first experiment we consider transfer function data from the ISS model perturbed by noise with mean µ = 0 and standard deviation σ 2 = 0.15.We employ LFaPP and enforce poles at ı[.77, 2, 4, 5.6, 9.33, 37.9] near peaks of the transfer function.The resulting real-valued surrogate model has order r = 12.The transfer functions of the surrogate model with enforced poles and reduced models computed from the same noisy data with LS-Loewner, Loewner-SVD, Loewner-CUR, and Modified AAA are given in Figure 3.
Enforcing the poles near peaks in the transfer function of the underlying data allows the surrogate to capture the behavior of the original data in a wider frequency range than applying LS-Loewner, Loewner-SVD, and Loewner-CUR.The choice of the locations, in which vicinity the poles should be chosen is, however, not automatized.Figure 3 also shows the relative errors of all surrogate models referenced to the original data without noise.While the enforced poles all have a negative real part, the models computed from the variants of the LF and AAA exhibit unstable eigenvalues.Thus, pole placement can be seen also as a means to enforce stability of the surrogate models.Alternatively, a post-processing step can be added to enforce stable models (for both LF and AAA methods), as performed in [29].We now evaluate the performance of the algorithms by applying them to heavily distorted transfer function measurements of the sound transmission problem.Noise with a standard deviation of σ 2 = 0.25 is considered and three algorithms are employed to compute surrogates: Loewner-SVD, LFPP (Algorithm 3), and LFaPP (Algorithm 4).We also test the modifications to LFaPP described in Section 3.4.These results are denoted by "LFaPP mod.".For LFPP we enforce poles at the eigenvalues of the underlying Loewner model which imaginary parts are near 2πı [72,189,392,401,706,856].These locations correspond to characteristic peaks in the transfer function.Further, we choose the in- terpolation points at 2πı [138, 339, 369, 569, 712, 954], which lie at the dips between the enforced poles.
Loewner-SVD and LFaPP do not require input parameters in addition to the measured data.Figure 4 shows the transfer function of the resulting surrogate models in comparison to the original and noisy underlying data.It can be observed, that the automatic approaches Loewner-SVD and LFaPP (mod.)cannot approximate the transfer function well after the first two peaks, i.e., for frequencies higher than 200 Hz, while LFPP approximates the original data over the complete frequency range with decent accuracy.The importance of reasonable interpolation points can be seen in the difference of LFPP and LFaPP mod., which have the same poles.It should be noted that the surrogate model computed by Loewner-SVD has two unstable poles while the other three surrogate models are stable.It is, however, not always clear a priori how to choose the poles and interpolation points for LFPP in order to achieve the best approximation quality possible.In this example, the noise level is too high for one of the automatic approaches to yield reasonable dominant interpolation points or poles.

Conclusion and outlook
In this contribution, we have proposed an extensive study of interpolation-based data-driven approaches for approximating the response of linear dynamical systems.All methods require input and output data, i.e., transfer function measurements, while direct access to the system operators or the states is not required.We showed different approaches how to achieve compact surrogate models approximating the input/output behavior of the original system and how to ensure various properties of the surrogate models, such as stability.Strategies how to work with noisy measurement data have also been addressed.
A natural extension of the framework described here is to apply the ideas of tangential interpolation as a means of modeling a MIMO system from data.Here, the tangential directions need to be incorporated in the parameterized one-sided realization.Further topics include enforcing different structures of the original model in the surrogate model, e.g., second-order or delay structures.It would also be interesting to study the possibility of placing certain stable poles while achieving interpolation in a least-squares sense.Application cases for the proposed methodology could include damping optimization.Here, a family of parameterized interpolants could be used to find optimal positions for viscous dampers in a structural system.

A. Appendix
A.1.The Woodbury matrix identity We can expand the right part of (19), such that: The Woodbury matrix identity is as follows: where M, Û, T and V are conformable matrices: M is n × n, T is k × k, Û is n × k, and V is k × n.This can be derived using blockwise matrix inversion.By denoting with Λ s = sI km − Λ, then the first Then, the next step is to choose projection matrices W, V ∈ C n×k as As explained in [3], the choice of W H above is explained by imposing the required poles for the reduced model, while V is chosen to match the interpolation conditions at the λ i 's.Moreover, using these notations, it follows that CV = 0. Next, put together the following matrices Ẽ = W H EV, Ã = W H AV.Then, it follows that (s Ẽ − Ã) loses rank when s ∈ {ζ 1 , . . ., ζ r }.To show this, we simply write Let H ζ (s) = C ζ (sE − A) −1 B be a rational function in s and we note that Ê and Â are a special type of diagonally scaled Cauchy matrices, with the following exact definition: From the definition in (39) Hence, the reduced-order linear descriptor system Σ pp : ( Ẽ, Ã, B, C) matches k interpolation conditions and has the required poles.Next, we show that this model can be written equivalently in the AF format.We first note that Ĉ = C.For next step, provided that the matrix Ẽ is non-singular, we remove it by incorporating it into the other matrices, as: Ȃ = Ẽ−1 Ã, B = Ẽ−1 B, Ȇ = I k , C = C.We note that the two realizations of the interpolatory ROM, i.e., ( Â, B, Ĉ) in ( 21 Hence, the above choice of vector B in (21) imposes the required poles.

Figure 1 :
Figure 1: The L ∞ errors of reduced-order models of order r computed from the ISS data.Left: SISO with first input and output, respectively.Right: MIMO with three inputs and three outputs m = p = 3 (the number of interpolation points is k = r m ).

Figure 2 :
Figure 2: Transfer function (top) and relative pointwise errors (bottom) for reduced-order models of size r = 108 for the aircraft model.The error is plotted only at frequencies which do not coincide to interpolation points of the respective method.

Figure 3 :
Figure 3: Transfer function of a surrogate with enforced poles compared to the noisy and original transfer function values.The transfer function of a model r = 12 computed from Loewner-SVD is given for reference.