Algorithms for the rational approximation of matrix-valued functions

A selection of algorithms for the rational approximation of matrix-valued functions are discussed, including variants of the interpolatory AAA method, the RKFIT method based on approximate least squares fitting, vector fitting, and a method based on low-rank approximation of a block Loewner matrix. A new method, called the block-AAA algorithm, based on a generalized barycentric formula with matrix-valued weights is proposed. All algorithms are compared in terms of obtained approximation accuracy and runtime on a set of problems from model order reduction and nonlinear eigenvalue problems, including examples with noisy data. It is found that interpolation-based methods are typically cheaper to run, but they may suffer in the presence of noise for which approximation-based methods perform better.


Introduction
Rational approximation is a powerful tool in applied science and engineering.To give just two examples, it is very commonly used for model order reduction [2][3][4] and the solution of nonlinear eigenvalue problems [16,Section 6].Recently, several new algorithms for the rational approximation and interpolation of scalar-valued functions have been proposed, including (in reverse chronological order) the AAA algorithm [21], the RKFIT algorithm [7,8], vector fitting [13,15], and methods based on the Loewner matrices [20].The aim of this paper is to explore extensions of these methods for the purpose of approximating a matrix-valued (or block) function F : Λ → C m×n on a discrete set Λ in the complex plane.
The paper contains two key contributions.Firstly, we propose an extension of a barycentric formula to non-scalar weights and develop it into a new algorithm for the computation of rational approximants of matrix-valued functions.This algorithm, referred to as the block-AAA algorithm, generalizes the AAA algorithm in [21].Secondly, we perform extensive numerical comparisons of several rational approximation algorithms proposed by different communities (model reduction, numerical linear algebra, approximation theory).To the best of our knowledge, this is the first paper that offers a comprehensive study of these algorithms.As part of our experiments, we identify potential problems of interpolation-based algorithms when the data samples are polluted by noise.
The paper is structured as follows.Section 2 contains a brief review of existing algorithms for scalar rational approximation.In Section 3 we introduce new representations of matrix-valued rational approximants and the new block-AAA algorithm.In Section 4 we briefly review other popular algorithms for rational approximation, including methods based on partial fractions, rational Krylov spaces, and state space representations.Section 5 contains a collection of numerical experiments and discussions, guided by six examples from different application areas.Section 6 further discusses the performance of the different methods for each numerical example.We conclude in Section 7 with some ideas for future work.

Scalar rational approximation
In this section we summarize several algorithms for the rational approximation of a scalar function f : Λ → C sampled at a nonempty set of points Λ := {λ 1 , . . ., λ } ⊂ C.This serves as an introduction of our notation, but also as a template for the block variants in Sections 3 and 4.

Adaptive Antoulas-Anderson (AAA) algorithm
The AAA algorithm proposed in [21] is a practically robust and easy-to-use method for scalar rational interpolation.The degree d interpolant r d obtained after d iterations of the AAA algorithm is of the form with nonzero barycentric weights w k ∈ C, pairwise distinct support points z k ∈ C, and function values f k = f (z k ).A key ingredient of the AAA algorithm is its greedy choice of the support points, one per iteration j = 0, 1, . . ., d, intertwined with the solution of a least squares problem to determine the barycentric weights w (j) 0 , . . ., w (j) j at each iteration j.The core AAA algorithm can be summarized as follows: 1. Set j = 0, Λ (0) := Λ, and r −1 :≡ −1 i=1 f (λ i ).
MATLAB implementations of the AAA algorithm can be found in [21] and the Chebfun package [11] available at https://www.chebfun.org/In [21, Section 10] a non-interpolatory variant of the AAA algorithm is also mentioned, which is obtained by allowing the weights in the numerator and denominator of (1) to be different.

Rational Krylov Fitting (RKFIT)
The RKFIT method introduced in [7,8] is based on iteratively relocating the poles of a rational function in the so-called RKFUN representation.An RKFUN is a triple , where n(z) is the unique left nullvector of zK d −H d normalized such that its first component is 1 (note that the unique existence of this nullvector is guaranteed by the fact that zK d −H d is an unreduced upper-Hessenberg matrix).
The matrices H d and K d of the RKFUN representation satisfy a rational Arnoldi decomposition AV d+1 K d = V d+1 H d associated with a rational Krylov space d+1) , and a monic polynomial q ∈ P d such that q(A) is invertible.It can be shown that the columns of V d+1 = [v 1 , . . ., v d+1 ] are all of the form v k+1 = p k (A)q(A) −1 b with polynomials p k ∈ P d for k = 0, 1, . . ., d.In other words, the rational Arnoldi decomposition encodes a basis of scalar rational functions r k := p k /q all sharing the same denominator q.Further, one can show that the roots ξ j of q(z) = ξ j =∞ (z − ξ j ) correspond to the quotients of the subdiagonal elements of the upper-Hessenberg matrices H d and K d , i.e., ξ j = h j+1,j /k j+1,j for j = 1, . . ., d.These quotients are referred to as the poles of the rational Arnoldi decomposition [7].
Given a matrix F ∈ C × , the RKFIT algorithm attempts to identify an RK-FUN d , one RKFIT iteration consists of solving a least squares problem for finding a unit-norm vector v ∈ span(V d+1 ) such that v = arg min and to use the roots of r = p/q associated with v = r(A)b as the new poles ξ d for the next iteration.This process is iterated until a convergence criterion is satisfied.

Vector Fitting (VF)
The VF algorithm, originally proposed in [15], seeks to fit a rational function in partial fraction (pole-residue) form The iterative algorithm is initiated by choosing the degree d of the rational approximant and an initial guess for the poles {ξ d }.At iteration j = 0, 1, . . .one determines parameters c is solved with least squares error over all z ∈ Λ. Afterwards, the next set of poles {ξ d } is computed as the roots of the polynomial q (j) (z) by solving a linear eigenvalue problem.The iteration continues until a convergence criterion is satisfied.Vector fitting is a non-interpolatory method.A MATLAB implementation of this method is available at https://www.sintef.no/projectweb/vectorfitting/

Loewner framework (LF)
This LF method uses the full data set in matrix format by forming (possibly very large) Loewner matrices.The original method, introduced in [20], is based on constructing a reduced-order rational function in state-space representation.In recent years, the Loewner framework has been extended to classes of mildly nonlinear systems, including bilinear systems [5].
Assuming that , the number of sampling points, is an even positive integer, the first step in the Loewner framework is to partition Λ = Λ L ∪ Λ R into two disjoint sets of the same cardinality.The set Λ L contains the left points {x 1 , . . ., x /2 }, while Λ R contains the right points {y 1 , . . ., y /2 }.Similarly, the set of scalar samples (evaluations of the function f on Λ) is partitioned into two sets.
One then defines matrices L ∈ C /2× /2 (the Loewner matrix) and L s ∈ C /2× /2 (the shifted Loewner matrix) with entries respectively.Additionally, one defines vectors The next step is to compute a rank-d truncated singular value decomposition of the Loewner matrix L ≈ XSZ * with X ∈ C /2×d , S ∈ C d×d and Z ∈ C /2×d .Finally, by means of projecting with the matrices X and Z, the fitted rational function is This is a subdiagonal rational function of type (d − 1, d) whose d poles are given by the generalized eigenvalues of the matrix pair (X * L s Z, X * LZ).
3 Matrix-valued barycentric forms and block-AAA The simplest matrix-valued barycentric form is obtained from (1) by replacing the function values f k with matrices Provided that all weights w k are nonzero, the function R d interpolates the function F at all the support points z k .Each (i, j) entry of R d is a rational function of the form p ij (z)/q(z), where p ij , q ∈ P d are polynomials of degree d.Note that all these entries share the same scalar denominator q.
A slight modification of (bary-A) yields a new matrix-valued barycentric formula with weight matrices W k ∈ C m×m .If all these W k are nonsingular, R d interpolates F at all the support points z k .Given a set of support points z 0 , . . ., z k different from any point in Λ, a linearized version of the approximation problem The weight matrices can now be obtained from a trailing left-singular block vector [W 0 , . . ., W d ] of unit norm such that Inputs: Discrete set Λ ⊂ C with points, function F , error tolerance ε > 0 Output: Rational approximant R j−1 in the form (bary-B) 1. Set j = 0, Λ (0) := Λ, and R −1 :≡ −1 i=1 F (λ i ).

Find matrix [W
j ] of unit Frobenius norm such that is solved with least squares error over all z ∈ Λ (j+1) .6. Set j := j + 1 and go to step 2.
Figure 1: Pseudocode of the block-AAA algorithm has smallest possible norm.Note that the matrix L is a block Loewner matrix. 1  It is now easy to derive an AAA-like algorithm based on (bary-B).It is shown in Figure 3 and referred to as the block-AAA algorithm.
A MATLAB implementation of the block-AAA algorithm is available at https://github.com/nla-group/block_aaa As it will become clearer from the discussions in Section 4, the block-AAA algorithm is different from the set-valued AAA [19] and the fast-AAA [18] algorithms: the entries of a block-AAA approximant R d do not share a common scalar denominator of degree d.Indeed, a block-AAA approximant R d of order d can have a larger number of up to dm poles, i.e., its McMillan degree can be as high as dm.For this reason we refer to the order of a block-AAA approximant instead of its degree.Finally, we mention for completeness two other barycentric formulas which might be of interest elsewhere.One is given as In the Loewner framework discussed in Section 2.4, the Loewner matrix entries are defined in terms of left and right directions (µi, i, vi) and (λj, rj, wj), but here we have a special case where the left and right directions i and rj are chosen as the unit vectors.
with the matrices C k ∈ C m×n and D k ∈ C m×m chosen such that R (λ i ) ≈ F (λ i ).We assume that the z k and λ i are pairwise distinct and consider the linearized problem This least squares problem is solved by a trailing left-singular block vector given by in general the approximation is non-interpolatory.Such approximants might be useful for the solution of nonlinear eigenvalue problems; see also the appendix.
It is also possible to have matrix-valued "support points" This leads to tangential interpolation determined by the eigenvalues and eigenvectors of the matrices Z k .We do not further explore this form here.

Other block methods
In this section we list methods for the rational approximation of matrix-valued functions using the barycentric formula (bary-A), as well as methods based on other representations of their approximants.These methods will be compared numerically in Sections 5 and 6.

Set-valued AAA (fast-AAA):
The set-valued AAA algorithm presented in [19], and similarly the fast-AAA algorithm in [18], applies the standard AAA algorithm to each component of F using common weights and support points for all of them, thereby effectively producing a barycentric interpolant in the form (bary-A).
Surrogate AAA: This method presented in [12] applies the standard AAA algorithm to a scalar surrogate function f (z) := a T F (z)b, with vectors a, b chosen at random.The resulting block interpolant is of the form (bary-A), obtained by replacing the scalar values

RKFIT (block):
The RKFIT algorithm has been generalized in [8] to the problem of finding a family r d are represented in the RKFUN format and share the common scalar denominator q of degree d.In our setting, we simply associate with each (i, j) entry of F a corresponding function f [k] (z) and run RKFIT on that family.Matrix Fitting: The extension of vector fitting to matrix-valued data has been implemented in the MATLAB Matrix Fitting Toolbox available at https://www.sintef.no/projectweb/vectorfitting/downloads/matrix-fitting-toolbox/ As described in the manual [14] associated with the toolbox, the primary intention of this software is the rational multi-port modelling of data in the frequency domain.
In particular, the toolbox deals with fitting so-called admittance values (known also as Y parameters) and scattering values (known also as S parameters).
The fitted rational function is provided in the pole-residue form or as an equivalent state space model and B ∈ C nd×n .For both representations, stable poles can be enforced.Note that the matrix elements of R d share a common set of d poles.
The above mentioned matrix fitting implementation [14] works with the elements of the upper triangular part of the matrix samples F i := F (λ i ) (i = 1, . . ., ), i.e., the entries [F i ] a,b for a ≤ b, which are stacked into a single vector.Hence, this implementation can be used only for symmetric matrix-valued functions.The elements of the vector are fitted using an implementation of the "relaxed version" of vector fitting [13], provided in the function vectfit3.m of the toolbox.The main driver function VFdriver.maccepts as input arguments sample points Λ = Λ = {λ 1 , . . ., λ } on the imaginary axis in complex conjugate pairs (see Section 3.1, page 7 of the toolbox manual).In our context of fitting on arbitrary discrete sets Λ in the complex plane, this requirement might represent an inconvenience.Also, the available driver seems to be intended to be run for even orders d only, although this is not an inherent restriction of vector fitting.Block Loewner framework: The extension of the Loewner framework to matrixvalued functions is covered in [20].Therein, the case of directionally-sampled matrix data is treated in detail.Introducing left directional vectors i ∈ C m and right directional vectors r j ∈ C n , one defines a Loewner matrix L ∈ C /2× /2 and a shifted Loewner matrix as follows: The definition of the matrices V ∈ C /2×n , W ∈ C m× /2 is also different from those in the scalar case in (2), namely the directional vectors appear in the form The rational matrix-valued function R d is then computed exactly as in the scalar case (3) in the form where the truncated singular vector matrices X, Z are the same as before.

Numerical comparisons
In all our tests we assume that Λ = {λ 1 , . . ., λ } and all the samples F (λ i ) are given.We will compare the six algorithms introduced in the previous sections, using MATLAB implementations provided by their authors whenever available.In order to compare the algorithms with respect to their approximation performance, we use the root mean squared error As before, d generally refers to the order of the rational approximant R d , but it is the same as the degree for all algorithms except block-AAA.The RMSE values are reported for various examples and different different orders d in Tables 1, 2, 3, 4, and 5.We also report timings for all algorithms under consideration.Although MAT-LAB timings might not always be reliable, the differences in the runtimes are usually large enough to give some indication of algorithmic complexity.In order to reduce random fluctuations in the timings, we run each algorithm for 20 times and compute the average execution time for a single run.All experiments were performed on a desktop computer with 8 GB RAM and an Intel(R) Core(TM) i7-7500U CPU running at 2.70 GHz.

Two toy examples
We first consider a rational function given by The entries of F can be written with a common scalar denominator of degree d = 6, e.g., q(z) = (z + 1)(z 2 + z − 5)(z 3 + 3z 2 − 1).We use = 100 logarithmically spaced samples of F on the imaginary axis in the interval [1,100]i.As in addition the samples of F are symmetric matrices, this example is suitable for the Matrix Fitting Toolbox [14].
The RMSE values obtained with all discussed algorithms for varying orders d = 0, 1, . . ., 20 are shown in Figure 2 (left).Note that all methods eventually recover the function accurately.The Loewner approach requires a slightly higher order than the other methods, i.e., d 1 = 8.The block-AAA algorithm correctly identifies F when the order is d 2 = 5.
Next, we modify the (1,2) entry of the function F by replacing its denominator by z 2 + z + 5.The entries of F can then be written with a common denominator of degree d = 8, e.g., q(z) = (z + 1)(z 2 + z − 5)(z 2 + z + 5)(z 3 + 3z 2 − 1).The RMSE values are shown in Figure 2 (right).Now the VF approach fails as the modified function F is no longer symmetric.The block-AAA algorithm again identifies F correctly with an order of d 1 = 5, while all other methods require a degree of d 2 = 8 as expected.

An example from the MF Toolbox
We now choose an example from the Matrix Fitting Toolbox [14].We use the file ex2 Y.m provided therein, containing 300 samples of 3 × 3 admittance matrices for a three-conductor overhead line.The sampling points are all on the imaginary axis within the interval 62.8319 • [1, 10 4 ]i. Figure 3 (left) depicts the magnitude of the nine matrix entries over the sampled frequency range.1.

A first model order reduction example
We now consider the CD player example from the SLICOT benchmark collection [1].The mechanism to be modelled is a swing arm on which a lens is mounted by means of two horizontal leaf springs.The challenge is to find a low-cost controller that can make the servo-system faster and less sensitive to external shocks.The LTI system that models the dynamics has 60 vibration modes, hence the dimension

A second model order reduction example
We consider another model reduction example, namely the ISS model from the SLICOT benchmark collection [1].There, an LTI system is used as structural model for the component 1R (Russian service module) in the International Space Station (ISS).The state space dimension of the linear system is 270 with m = n = 3 inputs and outputs.The matrix transfer function of this system is sampled at 400 logarithmically spaced points in the interval [10 −1 , 10 2 ]i.In Figure 5 (left) we depict the absolute value of the matrix entries.3.

A nonlinear eigenvalue problem
We consider the buckling example in [17], a 3×3 nonlinear eigenvalue problem that arises from a buckling plate model.Since we are interested in the approximation of non-rational functions by means of rational functions, we select only the nonconstant part of F (z). Hence, we consider the following 2 × 2 symmetric matrixvalued function We choose 500 logarithmically spaced sampling points in the interval [10 −2 , 10]i. Figure 6 shows the RMSE values achieved by the tested algorithms for varying orders d = 0, 1, . . ., 30.As before, we report the numerical RMSE values and the corresponding timings for two selected orders d = 10 and d = 20 in Table 4.

A scalar example with noise
In this experiment we investigate the effects of noisy perturbations on the approximation quality.As a test case, we use a scalar function f (z) = (z − 1)/(z 2 + z + 2).
Only two methods will be considered for this experiment, namely RKFIT (a noninterpolatory method) and AAA (an interpolatory method).We sample the function  f (z) at 500 logarithmically spaced points in the interval [10 −1 , 10 1 ]i, and then add normally distributed noise with a standard deviation of τ = 10 −2 to these samples.We compute rational approximants of degree d = 5 using the (scalar) RKFIT and AAA algorithms.For RKFIT, we perform three iterations starting with the default initialization of all poles at infinity.As it can be observed in Figure 7 (left), the red curve corresponding to the RKFIT approximants follows the noisy measurements very well on average.At the same time, the blue curve corresponding to the AAA approximant shows considerable deviations from the measurements.By inspecting the deviation between the two rational approximants and the original function f (z) in Figure 7 (right), we find that RKFIT approximation error for is comparable to the noise level.The RKFIT approximant has effectively estimated the additive noise rather accurately, which is also confirmed visually by the true noise curve overlaid as a dotted line on top of the RKFIT error curve.On the other hand, the approximation error attained by the AAA method is between 1-2 orders of magnitude larger than the added noise level.
Finally, we vary the degree d = 0, 1, . . ., 5 and compute the corresponding RMSE values of the RKFIT and AAA approximants.The results are depicted in Figure 8.Note that the RKFIT method achieves lower RMSE values than the AAA method and exhibits a more "regular" convergence.The RKFIT error stagnates at the noise level of approximately τ = 10 −2 when the degree d = 2 is reached.We believe that this preferable approximation performance is due to the non-interpolatory nature of RKFIT, and that AAA suffers from the fact that noisy values are being interpolated  (resulting in the observed oscillations of the approximants).

A second example with noise
For this case, we analyze a noisy version of the ISS model considered in Section 5.4.We modify the original LTI system (A, B, C) by choosing C = B T in order to obtain a symmetric transfer function and be able to apply the Matrix Fitting Toolbox.
The sample values are corrupted by additive normally distributed noise with standard deviation τ = 10 −2 .In Figure 9 (left) we depict the magnitude of each entry of the perturbed 3 × 3 matrix transfer function.
The RMSE values for varying orders d = 0, 1, . . ., 80 are plotted in Figure 9 (right).Now, for d = 10 and d = 20 we record the numerical RMSE values for all methods in Table 5.In accordance with our observations on the previous scalar example (Section 5.6) we find that the non-interpolatory methods like RKFIT and VF exhibit the most steady convergence behaviour, while the interpolation-based methods (AAA, Loewner) produce approximants whose accuracy varies wildly as the order d changes.The only methods that reliably attain an RMSE close to the noise level of τ = 10 −2 are RKFIT and VF, with RKFIT (10 iterations) consistently attaining the lowest RMSE values for all considered orders.
Recall that we considered the original (nonsymmetric) ISS model without noise in Section 5.4.As the timings of all methods are similar to those reported in Table 2, we do not report them separately for the noisy case considered here.

Discussion
this section we present a more detailed discussion of the numerical results reported in Section 5.The discussion will address each method separately, taking into consideration both the approximation quality and also the amount of time needed to run the method (runtime).

Set-valued AAA
A main benefit of the set-valued AAA algorithm is its fast runtime.With one exception, the set-valued AAA method was the second fastest method in our tests, surpassed only by the surrogate AAA method.
With respect to the approximation quality, this method produced models with similar accuracy to RKFIT for the example in Section 5.2 (see Figure 3, right).Similar behaviour is observed for the following three examples, producing good and very good results.For the example in Section 5.3 it even produced results comparable in accuracy to the block-AAA method when the order d was very high; see Figure 4 (right).For the examples in Sections 5.6 and 5.7 with added noise, however, the method produced poor results.In particular for the noisy ISS example, the RMSE values produced by the set-valued AAA method are the highest (see Figure 9, right).
As discussed above, we believe that the reason for this poor approximation is the interpolation of noisy data values.

Surrogate AAA
The surrogate AAA method performs best in terms of timings, and it might be an attractive approach when the problem dimension (m, n) is very large and approximation accuracy is not the main concern.The algorithm produces, in general, quite poor approximants compared to the other methods.For the example in Section 5.2, it did not reach an RMSE value below 10 −9 , while the set-valued AAA, RKFIT, and block-AAA methods were able to converge below 10 −13 ; see Figure 3 (right).A similar behaviour was observed for both MOR examples in Section 5.3 and 5.4.For the later, the method was the second worst (after VF).For the buckling plate example in Section 5.5 it produced the poorest approximation results (see Figure 6).Finally, for the ISS example with noise in Section 5.7, we observed that the method produced quite poor results, comparable to those produced by the Loewner framework.

Vector Fitting
For this method we used the implementation provided in [14].One limitation of this toolbox is that it can only be applied for symmetric matrix-valued functions, i.e., the samples must be symmetric matrices.Consequently, this implementation of VF could not be applied for two examples, namely the ones in Sections 5.3 and 5.4.
For the example in Section 5.2 we observed that the VF method produced poor approximants compared to some of the other methods.On the other hand, for the buckling plate example in Section 5.5 the VF method surpassed three methods and reached an RMSE of about 10 −12 for the degree d = 30 (see Figure 6).Note also that VF together with RKFIT are the only methods that reliably attain an RMSE value (τ = 0.03) close to that of the noise level (τ = 10 −2 ), as could be observed in Figure 9 (right).
Another observation is that it is generally not worth performing 10 VF iterations instead of just 5. We did not observe any significant approximation enhancement and the total runtime of VF will approximately double for 10 VF iterations.Throughout the experiments performed in Section 5, we found that VF was generally slower than the set-valued and surrogate AAA methods, but faster than the RKFIT and block-AAA methods.

RKFIT
For the RKFIT method we used the implementation in the RKToolbox.For each experiment, we ran the RKFIT algorithm for a fixed prescribed number of iterations.When stopping after five iterations, RKFIT was the second slowest method for the experiments in Sections 5.2, 5.3, and 5.5 (faster only than the block-AAA method) and the slowest for the experiments in Sections 5.4 and 5.7.
On the other hand, the RKFIT method proved to yield the second lowest RMSE values for the first five experiments, beaten only by block-AAA.This is not surprising given that RKFIT produces approximants R d where all matrix entries share a common scalar denominator of degree d, say, while a block-AAA approximant of the same order d can have up to dm singularities.Hence, among the scalar denominator methods, RKFIT was the most accurate.Moreover, RKFIT proved to be the most robust when dealing with noise.For both variants (with 5 and 10 iterations), the RMSE values obtained by RKFIT were the lowest among all methods (see Figure 9, right), comparable only to those of VF.
The runtime of the RKFIT method linearly depends on the number of iterations.It was observed in some cases that the RMSE values would not significantly decrease when increasing the number of RKFIT iterations from five to ten.For example, it can be seen in Figure 3 (right) and also in Figure 6 (right), that the two RMSE curves for RKFIT (5 iterations) and RKFIT (10 iterations) are practically indistinguishable.In practice, a dynamic stopping criterion based on the stagnation of the RMSE should therefore be used.

The Loewner framework
The Loewner framework, as introduced in Section 2.4, is the only direct method out of the six included in this study, i.e., it does not rely on an iteration.This could be an advantage with respect to the speed of execution.In general, this method performs very well in terms of running times for small, and medium to large data sets.The reason for this is that it relies on computing a full singular value decomposition of a matrix with dimension half of the data set.
It was observed that, with some exceptions, the Loewner framework was the third fastest method, surpassed only by the set-valued and surrogate AAA methods.This is not surprising since all the three methods rely on multiplying the matrix-valued samples with left and right tangential vectors, thereby reducing the dimension of the SVD problem to be solved at each iteration.
In terms of approximation accuracy, the Loewner framework produced rather poor results as compared to other methods.For example, in Section 5.2 it was observed that the RMSE values achieved by this method were only lower than those of VF (see Figure 3, right).For the experiments in Sections 5.4 and 5.5, it produced results better only than those of surrogate AAA.Note that, for the CD player example in Section 5.3, the approximation quality was comparable to that of the set-valued AAA method and the RKFIT (with 5 iterations).Finally, the Loewner method failed for the example with added noise (see Figure 9).

Block-AAA
The block-AAA method introduced in Section 4 is the only method in this comparison which produces rational approximants with nonscalar denominators.As such, the approximation quality obtained for a certain order d might not be fully comparable to that of another method.In terms of the execution times, the method proved to be the second slowest for the examples in Sections 5.2, 5.3, and 5.5 (in these cases, only the RKFIT method with 10 iterations was slower).Additionally, for the example presented in Section 5.4, the block-AAA method was indeed the slowest.The reason for this inefficiency is the repeated solution of the SVD problem in Step 5 of the algorithm (see Figure 3), involving a matrix of size m(j + 1) × n at each iteration j = 0, 1, . . ., d, as well as the greedy search in Step 2 which requires the evaluation of R j−1 (z) at all remaining sampling points in Λ (j) to determine the next interpolation point z j .We have tried to speed up Step 5 using the updating trick of the SVD decomposition described in [19], but noticed some numerical instabilities for sufficiently large degrees.We have therefore not used any SVD updating strategy for the block-AAA results reported here.
In terms of approximation accuracy, the block-AAA method usually produced the best results for a given order d, with the noisy ISS MOR example being the only exception.A representative illustration is Figure 3 (right), where block-AAA requires an order of d = 25 to reach an RMSE value of 10 −14 , while RKFIT and the set-valued AAA method require d = 50 (twice the order).The other methods do not even reach this accuracy.A similar behaviour is observed in Figure 4 (right) for the CD player example and in Figure 5 for the ISS example.For the noisy data experiment in Section 5.7, block-AAA was outperformed by RKFIT and, again, we believe this is due to RKFIT's non-interpolatory character.

Conclusion
We proposed an extension of the AAA algorithm to the case of matrix-valued interpolation.The block-AAA method uses a new generalized barycentric representation of the interpolant with matrix-valued weights.Our numerical experiments indicate that this modification allows for the computation of accurate approximants of low order.Having a low order might be particularly interesting in cases where the approximant needs to be linearized for further processing, e.g., in the context of nonlinear eigenvalue computations.In the appendix we show how the generalized barycentric interpolants can be linearized into a corresponding eigenvalue problem of size proportional to the order.The use of block-AAA for this application will be explored in future work.In terms of algorithmic complexity and runtime, the block-AAA method is currently inefficient and can be practically applied only for problems of small dimensions.Further work is required to improve its performance.One approach to deal with this shortcoming could be to replace the full singular value decomposition performed at each step of block-AAA with a CUR decomposition (selecting only a small number of relevant columns and rows in the block Loewner matrix), and to somehow restrict the search for the next interpolation point to a smaller set.
Our comparisons indicated that there was no method that always performed best with respect to accuracy and runtime.Interpolation-based methods are often cheaper to run, but they may suffer in the presence of noise for which approximationbased methods perform better.It seems fair to say that, despite the exciting recent developments in the area, there is still a lot of work to be done to design robust and fast methods for the rational approximation of matrix-valued functions.
while we obtain (bary-B) from (bary-C) by setting C k = W k F (z k ) and D k = W k .Hence it is sufficient to focus on the most general representation (bary-C).
If λ ∈ C is a nonlinear eigenvalue of R d (z), then clearly R d (z) −1 has a pole at z = λ, and vice versa.It is therefore sufficient to be able to find the nonlinear eigenvalues of R d defined in (bary-C) as those points λ ∈ C for which the "numerator" is a singular matrix.It is well known that for given support points z 0 , z 1 , . . ., z d , we can choose nonzero weights w 0 , w 1 , . . ., w d such that N (z) is an interpolating matrix polynomial: ; see, e.g., [9, eq.(3. 2)] and [22, eq. (3.4)].With this choice of the weights, we can immediately apply [22,Thm. 4.8], which states that the matrix pencil L(z) = L 0 − zL 1 , where with θ j = w j−1 /w j for j = 1, . . ., d is a strong linearization of N (z).In other words, the nonlinear eigenvalues of N (z) can be computed as the generalized eigenvalues of the matrix pair (L 0 , L 1 ).As long as R d does not also have a pole at a nonlinear eigenvalue λ of N (z), λ will be a nonlinear eigenvalue of R d (z).Finally, it is perhaps interesting to note that we have freedom to choose the nonzero weights w j and this choice will likely influence the numerical stability of the linearization (L 0 , L 1 ).This might be explored in some future work.

Figure 2 :
Figure 2: Accuracy comparison for the two toy examples in Section 5.1.Left: approximating a symmetric 2 × 2 matrix-valued function.Right: modified nonsymmetric case.

Figure 3 :
Figure 3: Example with admittance matrices from the MF Toolbox, described in Section 5.2.Left: entries of the 3 × 3 matrix F (z) over the frequency range.Right: accuracy performance.

Figure 3 (
Figure 3 (right) shows the RMSEs achieved by each of the tested algorithms for varying orders d = 0, 1, . . ., 50.For two particular orders, namely d = 10 and d = 20, we report the numerical RMSE values and the corresponding timings in Table1.

Figure 4 (
left) shows the magnitude of the 4 matrix entries over the frequency range.
Figure 4 (right) shows the RMSE values achieved by each of the algorithms for varying orders d = 0, 1, . . ., 80.The numerical RMSE values and the corresponding timings for orders d = 10 and d = 20 are listed in Table 2 .

Figure 4 :
Figure 4: player example from Section 5.3.Left: entries of the 2 × 2 matrix F (z) evaluated at 200 points on the imaginary axis.Right: accuracy performance.

Figure 5 (
right) shows the RMSE achieved by each of the algorithms for varying orders d = 0, 1, . . ., 80, and selected numerical RMSE values and corresponding timings for orders d = 10 and d = 20 are given in Table

Figure 5 :
Figure 5: ISS example from Section 5.4.Left: entries of the 3 × 3 matrix F (z) evaluated at 400 points on the imaginary axis.Right: accuracy performance.

Figure 6 :
Figure 6: The buckling problem in Section 5.5.Left: entries of the 2 × 2 (symmetric) matrix F (z) evaluated at 500 points on the imaginary axis.Right: accuracy performance.

Figure 7 :
Figure 7: Scalar example with noisy measurements from Section 5.6.Left: the original scalar function f (z) and the computed RKFIT and AAA approximants of order d = 5.Right: the deviation between the two approximants and the original function f (z).

Figure 8 :
Figure8: RMSE convergence of the approximants computed via RKFIT (noninterpolatory) and AAA (interpolatory) for increasing degree d.The original function is rational of degree 2, and noise with a standard deviation of τ = 10 −2 has been added to perturb it.

Table 1 :
Selected RMSE values and timings for all tested algorithms -MF Toolbox example The system's transfer function is a rational function of size 2×2 which we sample at 200 logarithmically spaced points in the interval [10 1 , 10 5 ]i.

Table 2 :
Selected RMSE values and timings for all tested algorithms -CD player example

Table 3 :
Selected RMSE values and timings for all tested algorithms -ISS example

Table 4 :
Selected RMSE values and timings for all tested algorithmsbuckling problem

Table 5 :
RMSE values for all tested methods -noisy ISS example