A tangential interpolation framework for the AAA algorithm

Interpolation‐based data‐driven methods, such as the Loewner framework or the Adaptive Antoulas‐Anderson (AAA) algorithm, are established and effective approaches to find a realization of a dynamical system from frequency response data (measurements of the system's transfer function). If a system‐theoretic representation of the original model is not available or unfeasible to evaluate efficiently, such reduced realizations enable effective analysis and simulation. This is especially relevant for models of interconnected dynamical systems, which typically have a high number of inputs and outputs to model their coupling conditions correctly.

domain is given by () ∈ ℂ × , with  ∈ ℂ, as where   is an identity matrix of size .In many applications,  is very large, hindering an efficient evaluation of Equation 1or Equation 2. Model order reduction methods seek approximations to the original dynamical system, whose output ŷ() matches the original output up to a specified tolerance  in an appropriate norm: ‖() − ŷ()‖ ≤  ⋅ ‖()‖.
A similar relation can be formulated for the transfer function Equation 2. If the dynamical system of interest is not given in state-space form Equation 1 or the state-space form cannot be evaluated conveniently, data-driven reduction or identification methods can be employed to obtain a reduced-order (surrogate) model which approximates data obtained from the original system.Some of these methods are based on interpolation [1,2] or optimization [3,4].
The Adaptive Antoulas-Anderson (AAA) algorithm [5] has emerged as a versatile tool for data-driven interpolation, combining ideas from interpolation and least-squares optimization.Originally a method for the interpolation of scalar functions, it has been extended to matrix functions in [6] by incorporating a block-wise approach.In the context of model order reduction and system inference, that is, when the matrix transfer function Equation 2 is to be interpolated from data samples, the block approach is limiting the effectiveness of the resulting surrogate model.Especially if systems with many inputs and outputs (,  ≥ 50) are considered, the resulting realizations quickly grow too large to be an efficient representation of the original system.In this contribution, we therefore modify AAA with ideas from tangential interpolation to be able to compute compact surrogates of the provided data.Tangential interpolation is widely used in projection-based model order reduction [7].
The remainder of the paper is structured as follows.We shortly review data-driven interpolation in the Loewner framework in Section 2, as it is the foundation of the proposed method.The theoretical foundation for the method is laid out in Section 3 and the resulting algorithm is presented in Section 4. We show the performance of the new method in numerical experiments in Section 5, and conclude the paper in Section 6.

DATA-DRIVEN INTERPOLATION IN THE LOEWNER FRAMEWORK
The Loewner framework [1] is an effective and established methodology that uses only input and output data, that is, transfer function measurements or evaluations, to find a surrogate model, which interpolates the transfer function of the original system Equation 2, and provides a realization of this interpolant as a state-space system.The general procedure is summarized, following [8].
Given  measurements   ∈ ℂ × ,  = 1, … , , of the transfer function at some locations   ∈ ℂ, the data are partitioned into two disjoint sets: In practical applications (when processing a fairly large number of measurements), the pencil (  , ) is typically singular or simply too large, so a post-processing step is required for the Loewner model Equation 4. Performing a singular value decomposition (SVD) of the pencil extracts the dominant features and removes inherent redundancies in the data.The left () and right () projection matrices are taken from the truncated SVDs where , ∈ ℝ × ,  ∈ ℂ × ,  ∈ ℂ × , Ỹ ∈ ℂ 2× , X ∈ ℂ ×2 .The truncation index  can, for example, be chosen as the numerical rank (based on a tolerance value  > 0).The projected Loewner system is then given by Ê =     , Â =     , B =     , Ĉ =   . (6)

ONE-SIDED TANGENTIAL INTERPOLATION AND ITS BARYCENTRIC FORMS
Assume, we want to interpolate the transfer function at  interpolation points, that is,  1 , … ,   .In what follows, the tangential directional vectors   ∈ ℂ  for 1 ≤  ≤  are included explicitly-this will be used to lower the dimension of the computed realizations.As opposed to full-matrix element-wise interpolation, we employ tangential interpolation.Note that the approach is equivalent to the full-matrix case, if the tangential directions are chosen as identity matrices.
We introduce the generalized reachability matrix with tangential directions, as corresponding to enforcing explicit interpolation, that is, matching of the following measurements: Here, the tangential vectors   ∈ ℂ ×1 are the so-called "tangential directions."The matrix  ∈ ℂ × solves a linear Sylvester equation  +  = , where the matrix pair (, ) explicitly determines the interpolation conditions enforced, that is, they are written as follows: Additionally, we consider so-called weight vectors   ∈ ℂ  for all 1 ≤  ≤ .For the moment, their entries will be considered as free parameters.As it will become evident in the subsequent sections, fixing them will be done in accordance with accommodating additional tangential interpolation conditions.
One-sided parametric realization: we show that an LTI parameterized model of dimension  =  can be constructed, having  degrees of freedom, to enforce one-sided interpolation in a realization-free manner.This means that no access to matrices , , and  is required.The free parameters are given by the entries of matrix B stored in the weight vectors ŵ 's.
Lemma 3.1.The transfer function of the reduced-order model (ROM) described by matrices above: given by: satisfies the following tangential interpolation conditions, with respect to the data in Equation 8, corresponding to the original system: Proof.The proof is quite straightforward, and it is based on a simple result.Let 1 ≤  ≤ .Then, by substituting  with   , and by multiplying the resolvent  − Â with the th unit vector   from the right, we can write that: Finally, by multiplying the above equality with Ĉ to the left, we obtain that:

□
Next, we show how to express the transfer function of the parameterized model that satisfies the  right interpolation conditions above, in a barycentric format.This is a necessary step for extending the classical AAA algorithm, to the tangential directions format, that is, in developing our tangential AAA algorithm.Lemma 3.2.With   =   − , the transfer function Equation 10of the one-sided interpolatory ROM in Equation 9, can be written in barycentric form as: Proof.To formally prove this result, we will make use of the Sherman-Morrison-Woodbury matrix inversion formula (e.g., [9]) The particular structure of the ROM allows the resolvent   − Â to be rewritten as follows: Using Equation 15, the transfer function Ĥ() is rewritten as: where X =  −1  B. Hence, the transfer function Ĥ() is explicitly written in the required format, that is: It follows that Next, we have that Putting together the last three identities proves Equation 14. □

THE TANGENTIAL AAA ALGORITHM: tAAA
The steps depicted in what follows go along the lines of the classical AAA algorithm and its block variant.Therefore, the proposed algorithm inherits the greedy structure of AAA, but we modify the least-squares problem, such that the new barycentric form and the vector-valued weights   are used.In this direction, assume we want to (approximately) enforce ℎ extra left tangential interpolation conditions given by: For any value of , so that 1 ≤  ≤ ℎ, we will use the barycentric formula of Ĥ() in Equation 14, in order to formulate equality constraints that will lead to solving a least-squares problem.Hence, we can write: ( Next, introduce the following (quasi) Loewner matrix  ∈ ℂ ℎ× , and also the data matrix Ẑ ∈ ℂ ℎ× containing the left tangentially interpolated data, and the matrix of variables B ∈ ℂ × (previously defined above): Lemma 4.1.By imposing the extra ℎ interpolation conditions in Equation ( 22), we need to solve the following linear system of ℎ equations in  variables written in matrix format, that is,  B = − F.
Remark 4.2.After solving the system of equations above, for matrix B, we put together the realization of the rational interpolant as in Equation 9.As in the single-input single-output (SISO) case, we need to impose that ℎ ≫ , so that the Loewner matrix in Equation 23 has more rows than columns.Hence, the system of equations will be solved only approximately.

A L G O R I T H M 1
The tangential AAA algorithm, tAAA.
Require: A (discrete) set of sample points Γ ⊂ ℂ with  points, function () (or the evaluations of () on the set Γ, i.e., the sample values), left and right tangential directions {  }  =1 ∈ ℂ  and {  }  =1 ∈ ℂ  , and an error tolerance  > 0. Ensure: A rational approximant () of order ( − 1, ) displayed in a barycentric form, and an equivalent dynamical system Σ in state space form.

11: end while
The resulting method in algorithm form is summarized in1.

NUMERICAL EXPERIMENTS
In the following, we demonstrate the methods discussed above by applying them on two benchmark examples available from the MOR-Wiki. 1We first consider a MIMO system with a moderate amount of inputs and outputs to be able to compare tAAA to bAAA, the block variant of AAA from [6].We also compute surrogate models using the Loewner framework in combination with an SVD (SVD-LF), as outlined in Section 2. The full potential of the tangential approach is afterward shown by applying it to compute a surrogate model for a dynamical system with  =  = 89 inputs and outputs.
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).

International Space Station
This system models the structural response of the Russian Service Module of the International Space Station (ISS) [10].
The model has  = 270 states,  = 3 inputs, and  = 3 outputs.The dataset used for the computations contains transfer function measurements at 400 logarithmically distributed points in the range  = [10 −1 , 10 2 ] ⋅ .We employ tAAA, bAAA, and SVD-LF to compute real-valued surrogate models with order  = 60 each.The transfer functions of the original model and its surrogates are given in Figure 1a, the corresponding point-wise relative errors are plotted in Figure 1b.All methods yield models with comparable accuracy, while the model obtained from SVD-LF has the lowest error and the model obtained from bAAA is least accurate.The slightly better accuracy of SVD-LF could be expected, as the SVD step considers all available data at once, while for the AAA variants, the iterative nature adds individual data points one by one.

Machine tool model
We now consider the model of a generic 3-axis machine tool with time-varying coupling conditions.Using a technique from [11], the coupling can be reformulated to a superposition of precomputed input and output mappings and therefore the model can be formulated as an LTI dynamical system.The model consists of two components ( 1 = 91181,  1 = 89,  1 = 89;  2 = 74392,  2 = 8,  2 = 8), which move relative to each other.More details about the system and the modeling process can be found in [12].Due to the high number of inputs and outputs in the first subsystem, bAAA cannot be applied here.The dataset used as input for tAAA and SVD-LF consists of the sampled transfer function at 200 logarithmically distributed frequencies in  = [1 × 10 −6 , 1 × 10 −2 ] ⋅ .We set the tolerance for tAAA to  = 1 × 10 −6 and truncate all states of the system obtained from the Loewner framework that correspond to singular values  <  max ⋅ 10 −14 .The surrogate model is investigated in the time domain, so all unstable poles have to be truncated in a post-processing step [12,13].
The resulting models are of orders  tAAA

2
= 69 for SVD-LF.The time responses of the original system and the two reduced-order models are given in Figure 2a, the corresponding relative errors are plotted in Figure 2b.
Both surrogate models approximate the original system with acceptable accuracy.The model computed with tAAA approximates the cooling phase of the system (starting after 600 s) slightly better than SVD-LF.

CONCLUSIONS
In this paper, we proposed tAAA, a modification of the AAA algorithm that incorporates ideas from tangential interpolation.This extends the capabilities of AAA-like methods to the interpolation of matrix functions of relatively large size.
In the context of this work, we used this property to compute surrogate models of multi-input, multi-output (MIMO) dynamical systems with many inputs and outputs.We demonstrated the versatility and effectiveness of the new method by comparing it to the classic Loewner framework and to a block variant of AAA.tAAA shows similar and in some cases even superior accuracy and can be applied in cases where a block-based strategy is not applicable.

1 2
Sigma plots of the transfer functions and relative errors of reduced-order models of the ISS model.ISS, International Space Station.Time responses and corresponding relative errors of surrogate models approximating the machine tool model.