∎

11institutetext: D. Palitta 22institutetext: Max Planck Institute for Dynamics
of Complex Technical Systems
Sandtorstraße 1, 39106 Magdeburg, Germany
22email<EMAIL_ADDRESS>33institutetext: S. Lefteriu
44institutetext: IMT Lille Douai, Institut Mines-Télécom, Univ. Lille, Centre
for Education, Research and Innovation (CERI) Digital Systems,
F-59000 Lille, France
44email<EMAIL_ADDRESS>
# An efficient, memory-saving approach for the Loewner framework

Davide Palitta  Sanda Lefteriu

(Received: date / Accepted: date)

###### Abstract

The Loewner framework is one of the most successful data-driven model order
reduction techniques. If $N$ is the cardinality of a given data set, the so-
called Loewner and shifted Loewner matrices $\mathbb{L}\in\mathbb{C}^{N\times
N}$ and $\mathbb{S}\in\mathbb{C}^{N\times N}$ can be defined by solely relying
on information encoded in the considered data set and they play a crucial role
in the computation of the sought rational model approximation. In particular,
the singular value decomposition of a linear combination of $\mathbb{S}$ and
$\mathbb{L}$ provides the tools needed to construct accurate models which
fulfill important approximation properties with respect to the original data
set. However, for highly-sampled data sets, the dense nature of $\mathbb{L}$
and $\mathbb{S}$ leads to numerical difficulties, namely the failure to
allocate these matrices in certain memory-limited environments or excessive
computational costs. Even though they do not possess any sparsity pattern, the
Loewner and shifted Loewner matrices are extremely structured and, in this
paper, we show how to fully exploit their Cauchy-like structure to reduce the
cost of computing accurate rational models while avoiding the explicit
allocation of $\mathbb{L}$ and $\mathbb{S}$. In particular, the use of the
_hierarchically semiseparable_ format allows us to remarkably lower both the
computational cost and the memory requirements of the Loewner framework
obtaining a novel scheme whose costs scale with $N\log N$.

###### Keywords:

Loewner framework data-driven model order reduction Cauchy-like matrices HSS
matrices

##  1 Introduction

The Loewner framework, originally proposed in artAJMACA  for solving the
generalized realization problem coupled with tangential interpolation, was
successfully employed for data-driven model order reduction from frequency
domain data SLACATCAD09 . Measurements of the frequency response are available
in several communities: electrical engineering (impedance, admittance or
scattering parameters SLACATCAD09 ), mechanical and civil engineering
(structural and vibro-acoustic frequency response functions PolyMAX  or
frequency response measurements of thermal systems Thermalsys ), to name a
few. The first step in the Loewner framework consists in setting up the data
matrices and building the Loewner and shifted Loewner matrices entry-wise
based on the chosen partition into right and left data, followed by computing
the singular value decomposition (SVD) of a linear combination of these
matrices and forming the model by projection, using the dominant singular
triplets. The main advantages of the Loewner framework over existing
approaches are, on the one hand, its system identification capabilities, in
the sense that the order of the system can be deduced from the singular value
drop, and, on the other hand, its potential in dealing with systems with a
large number of inputs and outputs efficiently, thanks to incorporating the
concept of tangential interpolation. The main drawbacks, however, are the
large storage requirements paired with the significant CPU cost inherent to
the full SVD computation for data sets with a large number of measurements
(values in the range $10^{5}$ are common in industrial applications). To
bypass these inconveniences, greedy-type approaches were proposed in
SLACATCAD09 , thus reducing memory requirements, from $\mathcal{O}(N^{2})$ for
storing the dense Loewner and shifted Loewner matrices to
$\mathcal{O}(N+n^{2})$, and the computational cost, from $\mathcal{O}(N^{3})$
for computing the SVD to $\mathcal{O}(Nn^{3})$ and $\mathcal{O}(Nn^{4})$,
where $N$ is the size of the data set and $n$ is the order of the model.

Taking advantage of numerical linear algebra tools to reduce storage and
computational requirements for the Loewner framework is another avenue worth
exploring due to the inherent structure embedded in the albeit dense Loewner
and shifted Loewner matrices. The factored ADI-Galerkin method for computing
these matrices as solutions to certain Sylvester equations with a factored
right-hand side was investigated in FSVDL . Such a scheme computes low-rank
approximations to the dense Loewner matrix to speed-up the SVD computation.
However, in FSVDL  no results about the accuracy of the computed reduced
models are reported. Moreover, the memory constraints coming from the
allocation of $\mathbb{L}$ and $\mathbb{S}$ are still present. Alternatively,
one can focus on accelerating solely the step of the SVD calculation by
employing Krylov methods (see, e.g., Sto12 ; Hoc01 ; BagRei05 ; Lar98  to name
a few), by using the randomized SVD nakatsukasa2020fast  to compute the
dominant singular triplets instead of the full SVD or other types of inexact
SVD-type decompositions (adaptive cross approximation ACA , particularly
suited for hierarchical matrices, or a CUR decomposition CUR  as in LoewCUR ;
9073015 ).

The novel approach proposed in this paper tackles the issue of the memory
requirements, at the same time as reducing the CPU cost of the Loewner
framework while maintaining the accuracy of the standard approach for large
values of the number of measurements. As the Loewner and shifted Loewner
matrices satisfy Sylvester equations with diagonal coefficient matrices, they
are, in fact, Cauchy-like matrices, obtained as the Hadamard product between a
Cauchy matrix $\mathcal{C}$ and low-rank right-hand sides. Extensive research
has been devoted to fully exploiting the rich structure of Cauchy matrices.
Several algorithms for computing the matrix-vector product $\mathcal{C}{\bf
x}$ can be found in the literature and many avoid assembling the full matrix
$\mathcal{C}$ (see, e.g., Pan2014 ; GREENGARD1987325 ; Gohberg1994 ;
Carrier1988 ). Hierarchically semiseparable matrices (HSS) have deemed
efficient for approximating Cauchy matrices with a low off-diagonal rank
PAN2015107 ; Pan2014 . HSS and other rank-structured matrices are widely used
in developing fast algorithms for algebraic operations (matrix-vector
multiplications, matrix factorizations, matrix inversion, etc., see, e.g.,
Vandebril2005 ; Pan2014 ; PAN2015107 ; Chandrasekaran2006 ; Xia2010  and
references therein) used as building blocks for the solution of certain
problems like linear systems of equations XiaCGL10 , eigenvalue problems
VogXCB16 , linear and quadratic matrix equations MassPR18 ; Kressner2019 , and
many more. For our application, the approximation of the Cauchy matrix in HSS
format considerably decreases the computational cost of matrix-vector products
involving a linear combination of the Loewner and shifted Loewner matrices
needed for the partial SVD computation, while avoiding to form them. All
results involving HSS-matrices presented in this paper have been obtained by
means of the hm-toolbox hmtoolbox .

The employment of an HSS-representation of $\mathcal{C}$ may introduce some
inexactness in our scheme and this has to be taken into account in the
iterative SVD computation. The use of inexact matrix-vector products within
iterative procedures has been the subject of numerous research papers: Krylov
techniques for solving linear systems and matrix equations Simoncini2003 ;
Bouras2005 ; Eshof2004 ; Kuerschner2019 ; Kuerschner2018 , eigenvalue problems
Freitag2007 ; Simoncini2002 , or an inexact variant of the Lanczos
bidiagonalization for the computation of leading singular triplets of a
generic matrix function GAAF2017 . In our case, we do not need an accurate
approximation of the singular triplets, but rather have meaningful spaces
spanned by the computed left and right singular vectors so that the obtained
reduced model inherits the desired approximation properties (see, e.g.,
IonitaPhD ; TutorialLoewner ).

The remainder of the paper is structured as follows. Section 2 provides a
review of the Loewner framework, whereas section 3 presents results showcasing
the special structure of the Loewner and shifted Loewner matrices as Cauchy-
like matrices and their approximation as hierarchically semiseparable matrices
allowing for efficient, inexact matrix-vector products in the partial SVD
computation. Section 4 presents the results of our numerical experiments and
section 5 concludes the paper.

##  2 Review of the Loewner framework

The Loewner framework has been proposed to address the rational
interpolation/approximation problem. In the control community, this is
referred to as system identification from frequency domain measurements and is
stated below.

######  Problem 1 (Rational approximation)

Given pairs of points representing the frequency $f_{j}$, and the
corresponding transfer function measurement at that frequency ${\bf
H}_{j}\in\mathbb{C}^{p\times q}$ for a system with $q$ inputs and $p$ outputs:

$(f_{j};{\bf H}_{j}),\leavevmode\nobreak\ \leavevmode\nobreak\
j=1,\ldots,N,\vspace{-.5mm}$ (1)

with $p$ and $q$ assumed to be much smaller than $N$, the problem amounts to
finding the rational transfer function ${\bf H}(\mathrm{s})$ which
approximates the data:

${\bf H}_{j}\approx{\bf
H}(\mathrm{s}=\mathrm{i}\omega_{j}),\leavevmode\nobreak\ \leavevmode\nobreak\
\forall j=1,\ldots,N.$ (2)

Thus, the transfer function evaluated for the Laplace variable
$\mathrm{s}=\mathrm{i}\omega_{j}$, where $\mathrm{i}^{2}=-1$ and
$\omega_{j}=2\pi f_{j}$, should be close (in some norm) to the corresponding
measurement ${\bf H}_{j}$. Several equivalent representations are possible for
the rational transfer function, namely pole-residue, pole-zero, state-space or
descriptor-form.

Most systems of interest are real, with their transfer function satisfying the
complex conjugate condition ${\bf H}(\bar{\mathrm{s}})=\overline{{\bf
H}(\mathrm{s})}$. Hence, we assume that the given data set satisfies this
condition and is of the following form:

$\left(\mathrm{i}\omega_{j},-\mathrm{i}\omega_{j};{\bf H}_{j},\overline{{\bf
H}}_{j}\right),\leavevmode\nobreak\ \leavevmode\nobreak\
j=1,\ldots,N.\vspace{-.5mm}$ (3)

We proceed by presenting the Loewner framework as a solution scheme addressing
the rational approximation Problem 1. The first step in the Loewner framework
artAJMACA ; SLACATCAD09  is partitioning the data in two disjoint sets. This
partition influences the conditioning of the problem (IonitaPhD, , Ch. 2.1)
and finding the optimal partition for each data set is beyond the scope of
this paper. The most natural partitions are summarized in the following
(assuming an even number of measurements $N$ and frequencies sorted in
ascending order):

  * •

Half&Half: the first half of the data in one set and the other half in the
second set:

$\left\\{\mathrm{i}\omega_{1},-\mathrm{i}\omega_{1},\ldots,\mathrm{i}\omega_{\frac{N}{2}},-\mathrm{i}\omega_{\frac{N}{2}}\right\\}\cup\left\\{\mathrm{i}\omega_{\frac{N}{2}+1},-\mathrm{i}\omega_{\frac{N}{2}+1},\ldots,\mathrm{i}\omega_{N},-\mathrm{i}\omega_{N}\right\\}$
(4)

and, correspondingly,

$\left\\{{\bf H}_{1},\overline{{\bf H}}_{1},\ldots,{\bf
H}_{\frac{N}{2}},\overline{{\bf H}}_{\frac{N}{2}}\right\\}\cup\left\\{{\bf
H}_{\frac{N}{2}+1},\overline{{\bf H}}_{\frac{N}{2}+1},\ldots,{\bf
H}_{N},\overline{{\bf H}}_{N}\right\\},$ (5)

  * •

Odd&Even: data with odd indices in the first set and data with even indices in
the second set:

$\left\\{\mathrm{i}\omega_{1},-\mathrm{i}\omega_{1},\ldots,\mathrm{i}\omega_{N-1},-\mathrm{i}\omega_{N-1}\right\\}\cup\left\\{\mathrm{i}\omega_{2},-\mathrm{i}\omega_{2},\ldots,\mathrm{i}\omega_{N},-\mathrm{i}\omega_{N}\right\\}$
(6)

and, correspondingly,

$\left\\{{\bf H}_{1},\overline{{\bf H}}_{1},\ldots,{\bf
H}_{N-1},\overline{{\bf H}}_{N-1}\right\\}\cup\left\\{{\bf
H}_{2},\overline{{\bf H}}_{2},\ldots,{\bf H}_{N},\overline{{\bf
H}}_{N}\right\\}.$ (7)

The first set on the right in (4) and (6) comprises the _right points_ ,
denoted by $\lambda_{k}$ , $k=1,\ldots,N$, while the second set comprises the
_left points_ $\mu_{h}$ , $h=1,\ldots,N$.

The following step in the Loewner framework is choosing tangential directions
as vectors which transform matrix data ${\bf H}_{j}$ into vector data: _right
tangential directions_ are column vectors ${\bf r}_{k}\in\mathbb{C}^{q}$ such
that ${\bf H}_{k}{\bf r}_{k}={\bf w}_{k}$, whereas _left tangential
directions_ are row vectors $\ell_{h}\in\mathbb{C}^{1\times p}$ such that
$\ell_{h}{\bf H}_{h}={\bf v}_{h}$. The column vectors ${\bf
w}_{k}\in\mathbb{C}^{p}$ are referred to as _right vector data_ , while the
row vectors ${\bf v}_{h}\in\mathbb{C}^{1\times q}$ are referred to as _left
vector data_. For simplicity, tangential directions can be chosen as
alternating columns/rows of the identity matrix SLACATCAD09 , resulting in
vector data being column and row vectors of the original matrix data ${\bf
H}_{j}$ in (1).

###### Remark 1

For scalar data obtained from single-input single-output (SISO) systems
($p=q=1$), tangential directions ${\bf r}_{k}$, $\ell_{h}$ are simply equal to
1.

###### Remark 2

If the loss of information due to utilizing a single tangential direction per
measurement, instead of the whole matrix ${\bf H}_{j}$, does not allow to
obtain an accurate approximation, one can employ the original matrix ${\bf
H}_{j}$. This is equivalent to considering several tangential directions for
the same point. To obtain block right matrix data for ${\bf
H}_{j}\in\mathbb{C}^{p\times q}$, the corresponding frequency should be
repeated $q$ times as a right point and all columns of the identity matrix of
size $q\times q$ should be considered as right directions. Similarly, to
obtain block left matrix data for ${\bf H}_{j}\in\mathbb{C}^{p\times q}$, the
corresponding frequency should be repeated $p$ times as a left point and all
rows of the identity matrix of size $p\times p$ should be considered as left
directions.

With this notation in place, the Loewner matrix is defined entry-wise as

$\mathbb{L}_{hk}=\frac{{\bf v}_{h}{\bf r}_{k}-\ell_{h}{\bf
w}_{k}}{\mu_{h}-\lambda_{k}},\leavevmode\nobreak\ \leavevmode\nobreak\
h,k=1,\ldots,N,$ (8)

and the shifted Loewner matrix is defined as

$\mathbb{S}_{{hk}}=\frac{\mu_{h}{\bf v}_{h}{\bf r}_{k}-\lambda_{k}\ell_{h}{\bf
w}_{k}}{\mu_{h}-\lambda_{k}},\leavevmode\nobreak\ \leavevmode\nobreak\
h,k=1,\ldots,N.$ (9)

Note that the numerators are scalar quantities as they are obtained by taking
inner products.

The quantities defined previously are collected into the following matrices

$\displaystyle\bm{\Lambda}$
$\displaystyle=\mbox{diag}\left(\left[\begin{array}[]{ccc}\lambda_{1},&\ldots&,\lambda_{N}\end{array}\right]\right)\in\mathbb{C}^{N\times
N},$ $\displaystyle{\bf R}$ $\displaystyle=\left[\begin{array}[]{ccc}{\bf
r}_{1},&\ldots&,{\bf r}_{N}\end{array}\right]\in\mathbb{C}^{q\times N},$
$\displaystyle{\bf W}$ $\displaystyle=\left[\begin{array}[]{ccc}{\bf
w}_{1},&\ldots&,{\bf w}_{N}\end{array}\right]\in\mathbb{C}^{p\times N},$ (13)
$\displaystyle\bm{M}$
$\displaystyle=\mbox{diag}\left(\left[\begin{array}[]{ccc}\mu_{1},&\ldots&,\mu_{N}\end{array}\right]\right)\in\mathbb{C}^{N\times
N},$ $\displaystyle{\bf L}$ $\displaystyle=\left[\begin{array}[]{c}\ell_{1}\\\
\vdots\\\ \ell_{N}\end{array}\right]\in\mathbb{C}^{N\times p},$
$\displaystyle{\bf V}$ $\displaystyle=\left[\begin{array}[]{c}{\bf v}_{1}\\\
\vdots\\\ {\bf v}_{N}\end{array}\right]\in\mathbb{C}^{N\times q},$ (21)
$\displaystyle\mathbb{L}$ $\displaystyle=\left[\begin{array}[]{ccc}\frac{{\bf
v}_{1}{\bf r}_{1}-\ell_{1}{\bf w}_{1}}{\mu_{1}-\lambda_{1}}&\ldots&\frac{{\bf
v}_{1}{\bf r}_{N}-\ell_{1}{\bf w}_{N}}{\mu_{1}-\lambda_{N}}\\\
\vdots&\ddots&\vdots\\\ \frac{{\bf v}_{N}{\bf r}_{1}-\ell_{N}{\bf
w}_{1}}{\mu_{N}-\lambda_{1}}&\ldots&\frac{{\bf v}_{N}{\bf r}_{N}-\ell_{N}{\bf
w}_{N}}{\mu_{N}-\lambda_{N}}\end{array}\right]\in\mathbb{C}^{N\times N},$
$\displaystyle\mathbb{S}$ $\displaystyle=\left[\begin{array}[]{ccc}\frac{{\bf
v}_{1}{\bf r}_{1}-\ell_{1}{\bf w}_{1}}{\mu_{1}-\lambda_{1}}&\ldots&\frac{{\bf
v}_{1}{\bf r}_{N}-\ell_{1}{\bf w}_{N}}{\mu_{1}-\lambda_{N}}\\\
\vdots&\ddots&\vdots\\\ \frac{{\bf v}_{N}{\bf r}_{1}-\ell_{N}{\bf
w}_{1}}{\mu_{N}-\lambda_{1}}&\ldots&\frac{{\bf v}_{N}{\bf r}_{N}-\ell_{N}{\bf
w}_{N}}{\mu_{N}-\lambda_{N}}\end{array}\right]\in\mathbb{C}^{N\times N}.$ (28)

By construction, the Loewner and shifted Loewner matrices satisfy the
following Sylvester equations:

$\displaystyle\bm{M}\mathbb{L}-\mathbb{L}\bm{\Lambda}$ $\displaystyle={\bf
V}{\bf R}-{\bf L}{\bf W},$
$\displaystyle\bm{M}\mathbb{S}-\mathbb{S}\bm{\Lambda}$
$\displaystyle=\bm{M}{\bf V}{\bf R}-{\bf L}{\bf W}\bm{\Lambda},$ (29)

as well as the following relations:

$\displaystyle\mathbb{S}-\mathbb{L}\bm{\Lambda}$ $\displaystyle={\bf V}{\bf
R},$ $\displaystyle\mathbb{S}-\bm{M}\mathbb{L}$ $\displaystyle={\bf L}{\bf
W},$ (30)

which will prove useful in our proposed matrix-free matrix-vector product
approach.

Assuming that the data is generated from a real system (3), to avoid complex
arithmetic, a change of basis can be performed. By defining

$\bm{\Pi}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{rr}1&-\mathrm{i}\\\
1&\mathrm{i}\end{array}\right]\mbox{ and }{\bf
P}=\mbox{blkdiag}\left(\left[\begin{array}[]{ccc}\bm{\Pi},&\ldots&,\bm{\Pi}\end{array}\right]\right)\in\mathbb{C}^{N\times
N},$ (31)

we obtain matrices with real entries:

$\bm{\Lambda}_{r}:={\bf P}^{*}\bm{\Lambda}{\bf P},\leavevmode\nobreak\
\bm{M}_{r}:={\bf P}^{*}\bm{M}{\bf P},\leavevmode\nobreak\ {\bf L}_{r}:={\bf
P}^{*}{\bf L},\leavevmode\nobreak\ {\bf V}_{r}:={\bf P}^{*}{\bf
V},\leavevmode\nobreak\ {\bf R}_{r}:={\bf R}{\bf P},\leavevmode\nobreak\ {\bf
W}_{r}:={\bf W}{\bf P},\leavevmode\nobreak\ \mathbb{L}_{r}:={\bf
P}^{*}\mathbb{L}{\bf P},\leavevmode\nobreak\ \mathbb{S}_{r}:={\bf
P}^{*}\mathbb{S}{\bf P},$

where ${\bf P}^{*}$ stands for the complex conjugate transpose of the matrix
${\bf P}$ and ${\bf P}^{-1}={\bf P}^{*}$. These quantities satisfy the same
equations as in (29) and (30). Unfortunately, $\bm{\Lambda}_{r}$ and
$\bm{M}_{r}$ are no longer diagonal and this represents a major drawback in
taking advantage of the Sylvester equations (29) for a fast computation of
$\mathbb{L}_{r}$ and $\mathbb{S}_{r}$. However, $\bm{\Lambda}_{r}^{2}$ and
$\bm{M}_{r}^{2}$ are diagonal and given by

$\displaystyle\bm{\Lambda}_{r}^{2}$
$\displaystyle=\mbox{blkdiag}\left[\begin{array}[]{rr}-\omega_{k}^{2}&0\\\
0&-\omega_{k}^{2}\end{array}\right],\;k=1,\ldots,N,$
$\displaystyle\bm{M}_{r}^{2}$
$\displaystyle=\mbox{blkdiag}\left[\begin{array}[]{rr}-\omega_{h}^{2}&0\\\
0&-\omega_{h}^{2}\end{array}\right],\;h=1,\ldots,N.$ (36)

By multiplying the first equation in (29) by $\bm{M}_{r}$ on the left and,
afterwards, multiplying it by $\bm{\Lambda}_{r}$ on the right and adding the
results together, a new Sylvester equation with diagonal coefficient matrices
is obtained:

$\displaystyle\bm{M}_{r}^{2}\mathbb{L}_{r}-\mathbb{L}_{r}\bm{\Lambda}_{r}^{2}$
$\displaystyle=\bm{M}_{r}\left({\bf V}_{r}{\bf R}_{r}-{\bf L}_{r}{\bf
W}_{r}\right)+\left({\bf V}_{r}{\bf R}_{r}-{\bf L}_{r}{\bf
W}_{r}\right)\bm{\Lambda}_{r}.$ (37)

By performing the same operations on the second equation in (29), a similar
Sylvester equation is obtained for the shifted Loewner matrix:

$\displaystyle\bm{M}_{r}^{2}\mathbb{S}_{r}-\mathbb{S}_{r}\bm{\Lambda}_{r}^{2}$
$\displaystyle=\bm{M}_{r}\left(\bm{M}_{r}{\bf V}_{r}{\bf R}_{r}-{\bf
L}_{r}{\bf W}_{r}\bm{\Lambda}_{r}\right)+\left(\bm{M}_{r}{\bf V}_{r}{\bf
R}_{r}-{\bf L}_{r}{\bf W}_{r}\bm{\Lambda}_{r}\right)\bm{\Lambda}_{r}.$ (38)

In the following, we say we employ the Odd&Even (real) partition whenever the
approach above is adopted.

After introducing notation, we are ready to state the solution provided by the
Loewner framework to the approximation Problem 1. A (non minimal) model for
the transfer function in descriptor-form ${\bf H}(\mathrm{s})={\bf
C}\left(\mathrm{s}{\bf E}-{\bf A}\right)^{-1}{\bf B}+{\bf D}$ is given by

${\bf H}(s)={\bf W}\left(\mathbb{S}-s\mathbb{L}\right)^{-1}{\bf V},$ (39)

or, alternatively, by ${\bf H}(s)={\bf
W}_{r}\left(\mathbb{S}_{r}-s\mathbb{L}_{r}\right)^{-1}{\bf V}_{r}$ if real
arithmetic was enforced. As we have recast the original problem as a
tangential interpolation problem, this transfer function satisfies the right
and left interpolation conditions artAJMACA  ${\bf H}(\lambda_{k}){\bf
r}_{k}={\bf w}_{k}$ and $\ell_{h}{\bf H}(\mu_{h})={\bf v}_{h}$,
$h,k=1,\ldots,N$ exactly. To obtain a minimal model, we perform a singular
value decomposition

$[{\bf Y},\bm{\Sigma},{\bf
X}]=\mbox{svd}(\mathbb{S}-x\mathbb{L}),\leavevmode\nobreak\
\leavevmode\nobreak\ x\in\\{f_{i}\\},$ (40)

where $\bm{\Sigma}$ is diagonal and ${\bf Y}$, ${\bf X}$ contain the left and
right singular vectors, respectively. Choosing the order $n$ of the truncated
SVD ($n$ is application-dependent), we define (in Matlab notation) ${\bf
X}_{n}={\bf X}(:,1:n)$ and ${\bf Y}_{n}={\bf Y}(:,1:n)^{*}$. Finally, the
model of size $n$ in descriptor form is

${\bf E}=-{\bf Y}_{n}\mathbb{L}{\bf
X}_{n}=-\mathbb{L}_{n},\leavevmode\nobreak\ {\bf A}=-{\bf Y}_{n}\mathbb{S}{\bf
X}_{n}=-\mathbb{S}_{n},\leavevmode\nobreak\ {\bf B}={\bf Y}_{n}{\bf V}={\bf
V}_{n},\leavevmode\nobreak\ {\bf C}={\bf W}{\bf X}_{n}={\bf
W}_{n},\leavevmode\nobreak\ {\bf D}={\bf 0}.$ (41)

When employing real arithmetic, the SVD trunctation step is analogous, in
terms of the $\mathbb{L}_{r}$, $\mathbb{S}_{r}$, ${\bf V}_{r}$ and ${\bf
W}_{r}$ matrices. In the following section, we exploit the Cauchy-like
structure of the Loewner and shifted Loewner matrices to design efficient
approaches, both in terms of memory storage and CPU time, to compute the SVD
in (40) by making use of hierarchical matrices.

##  3 Exploiting the structure of $\mathbb{L}$ and $\mathbb{S}$

For data sets with a sizable number $N$ of measurements $\mathbf{H}_{j}$, the
construction of the large, dense Loewner and shifted Loewner matrices is
demanding, both in terms of computational efforts as well as storage
requirements. The computation of each entry of $\mathbb{L}$ and $\mathbb{S}$
using (8) and (9) yields a total cost of $\mathcal{O}(N^{2}\left(p+q\right))$
floating point operations (FLOPs) for assembling the entire $\mathbb{L}$ and
$\mathbb{S}$ matrices. The number of nonzero entries in $\mathbb{L}$ and
$\mathbb{S}$ is $\mathcal{O}(N^{2})$, much larger than the memory requirements
for storing the data in $\bm{\Lambda}$, $\bm{M}$, ${\bf R}$, ${\bf W}$, ${\bf
L}$, and ${\bf V}$111The number of nonzero entries in the data matrices ${\bf
V}$ and ${\bf R}$ amounts to $\mathcal{O}(qN)$ and to $\mathcal{O}(pN)$ for
${\bf W}$ and ${\bf L}$.. Besides these excessive storage requirements, there
are also considerations to be made regarding the CPU time required for the SVD
computation of the matrix $\mathbb{S}-x\mathbb{L}$, $x\in\\{f_{i}\\}$ in (40).
Especially for large dimensional problems, for which we expect a fast decay,
it is preferred to compute only the first $n$ singular triplets, thus avoiding
wasting resources in computing the full SVD. To this end, many iterative
methods have been developed for computing partial SVDs; see, e.g., Sto12 ;
Hoc01 ; BagRei05 ; Lar98  to name a few. The bottleneck in these approaches is
the matrix-vector product with the coefficient matrix, namely
$\mathbb{S}-x\mathbb{L}$ in our case. This operation costs
$\mathcal{O}(N^{2})$ FLOPs due to the dense pattern of
$\mathbb{S}-x\mathbb{L}$.

This section tackles the cost reduction of performing a matrix-vector product
with $\mathbb{S}-x\mathbb{L}$ while avoiding the explicit allocation of
$\mathbb{L}$ and $\mathbb{S}$. The proposed strategy is supported by a
thorough analysis of the computational cost, showing that, for very large data
sets for which carrying out the full SVD is intractable, our strategy leads to
remarkable reductions in both the computational efforts and the storage demand
for building minimal realizations in the Loewner framework.

###  3.1 Hadamard product and Cauchy matrices

We present novel results which exploit the particular structure of the Loewner
and shifted Loewner matrices. These developments involve the Sylvester
equations (29) with diagonal coefficient matrices $\bm{\Lambda}$ and ${\bf
M}$.

###### Theorem 1

The Loewner and shifted Loewner matrices $\mathbb{L}$ and $\mathbb{S}$
satisfying the Sylvester equations in (29) are such that

$\mathbb{L}=\sum_{j=1}^{q}\mbox{diag}(\widetilde{\bf
v}_{j})\mathcal{C}\mbox{diag}(\widetilde{\bf
r}_{j}^{*})-\sum_{j=1}^{p}\mbox{diag}(\widetilde{\ell}_{j})\mathcal{C}\mbox{diag}(\widetilde{\bf
w}_{j}^{*}),$ (42)

and

$\mathbb{S}=\sum_{j=1}^{q}\mbox{diag}(\bm{M}\widetilde{\bf
v}_{j})\mathcal{C}\mbox{diag}(\widetilde{\bf
r}_{j}^{*})-\sum_{j=1}^{p}\mbox{diag}(\widetilde{\ell}_{j})\mathcal{C}\mbox{diag}(\bm{\Lambda}^{*}\widetilde{\bf
w}_{j}^{*}),$ (43)

where $\mathcal{C}$ denotes the following Cauchy matrix

$\displaystyle\mathcal{C}=\left[\begin{array}[]{ccc}\frac{1}{\mu_{1}-\lambda_{1}}&\ldots&\frac{1}{\mu_{1}-\lambda_{N}}\\\
\vdots&\ddots&\vdots\\\
\frac{1}{\mu_{N}-\lambda_{1}}&\ldots&\frac{1}{\mu_{N}-\lambda_{N}}\end{array}\right],$

while the vectors $\widetilde{\bf v}_{j}\in\mathbb{C}^{N}$ and
$\widetilde{\ell}_{j}\in\mathbb{C}^{N}$ denote the $j$-th columns of ${\bf V}$
and ${\bf L}$, respectively, so that

$\bm{V}=[\widetilde{\bf v}_{1},\ldots,\widetilde{\bf
v}_{q}],\quad\bm{L}=[\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{q}].$

Similarly, the vectors $\widetilde{\bf r}_{j}\in\mathbb{C}^{1\times N}$ and
$\widetilde{\bf w}_{j}\in\mathbb{C}^{1\times N}$ are the $j$-th rows of ${\bf
R}$ and ${\bf W}$, respectively, namely

$\bm{R}=\left[\begin{array}[]{c}\widetilde{\bf r}_{1}\\\ \vdots\\\
\widetilde{\bf r}_{q}\\\
\end{array}\right],\quad\bm{W}=\left[\begin{array}[]{c}\widetilde{\bf
w}_{1}\\\ \vdots\\\ \widetilde{\bf w}_{q}\\\ \end{array}\right].$

###### Proof

The Loewner and shifted Loewner matrices $\mathbb{L}$ and $\mathbb{S}$ are
Cauchy-like matrices as they are obtained by taking the Hadamard product
$\circ$ between the Cauchy matrix $\mathcal{C}$ and the right-hand sides of
the Sylvester equations in (29). In particular,

$\displaystyle\mathbb{L}$ $\displaystyle=\mathcal{C}\circ\left({\bf V}{\bf
R}-{\bf L}{\bf W}\right),$ $\displaystyle\mathbb{S}$
$\displaystyle=\mathcal{C}\circ\left(\bm{M}{\bf V}{\bf R}-{\bf L}{\bf
W}\bm{\Lambda}\right).$ (44)

An important property of the Hadamard product reads as follows. For any
vectors ${\bf x},{\bf y}\in\mathbb{C}^{N}$, it holds

$\mathcal{C}\circ({\bf x}{\bf y}^{*})=\text{diag}({\bf
x})\mathcal{C}\text{diag}({\bf y}^{*}).$

This, along with the low-rank structure of ${\bf V}{\bf R}-{\bf L}{\bf W}$ and
$\bm{M}{\bf V}{\bf R}-{\bf L}{\bf W}\bm{\Lambda}$, yields the results in (42)
and (43).

###### Corollary 1

Given a vector ${\bf y}\in\mathbb{C}^{N}$ and $x\in\\{f_{i}\\}$, we have

$(\mathbb{S}-x\mathbb{L}){\bf y}=\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ(\mathcal{C}(\widetilde{\bf r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ(\mathcal{C}(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf y},$

where $\widetilde{\bf y}=(\bm{\Lambda}-x{\bf I}){\bf y}$, with ${\bf I}$, the
identity matrix.

###### Proof

Thanks to (30), we can write

$\displaystyle(\mathbb{S}-x\mathbb{L}){\bf y}=$
$\displaystyle(\mathbb{L}\bm{\Lambda}+{\bf V}{\bf R}-x\mathbb{L}){\bf
y}=\mathbb{L}(\bm{\Lambda}-x{\bf I}){\bf y}+{\bf V}{\bf R}{\bf y}.$

The result follows by substituting the expression of $\mathbb{L}$ given in
Theorem 1 in the equation above.

Similar results to those in Theorem 1 and Corollary 1 can be obtained for
$\mathbb{L}_{r}$ and $\mathbb{S}_{r}$ solving the Sylvester equations in (37)
and (38), respectively. The developments follow the same line of proof as
above with straightforward adjustments.

Corollary 1 shows that the majority of the computational cost of performing
the matrix-vector multiplication $(\mathbb{S}-x\mathbb{L}){\bf y}$ amounts to
computing $p+q$ matrix-vector products with the Cauchy matrix $\mathcal{C}$.

Extensive research has been devoted to fully exploiting the rich structure of
Cauchy matrices. Several algorithms for computing the matrix-vector product
$\mathcal{C}{\bf y}$ can be found in the literature and many avoid assembling
the full matrix $\mathcal{C}$ (see, e.g., Pan2014 ; GREENGARD1987325 ;
Gohberg1994 ; Carrier1988 ). In the next section we recall the strategy
presented by Pan in PAN2015107 ; Pan2014  to represent $\mathcal{C}$ in terms
of a hierarchically semiseparable (HSS) matrix. Even though the novel scheme
proposed in this paper does not depend on the strategy employed for performing
the matrix-vector product $\mathcal{C}{\bf y}$ – as long as it is efficient –
we believe that the HSS framework may be advantageous as, in principle, many
matrix-vector products with $\mathcal{C}$ are needed for computing a (partial)
SVD of the matrix $\mathbb{S}-x\mathbb{L}$.

We conclude this section with the following remarks.

###### Remark 3

The number $n$ of singular triplets needed to be computed to achieve the
minimal realization $({\bf E},{\bf A},{\bf B},{\bf C},{\bf D})$ in (41) is
difficult to estimate a-priori222In (Beckermann2019, , Section 4.3), some
results on the numerical rank of $\mathbb{L}$ are presented provided
$[\min\lambda_{k},\max\lambda_{k}]\cap[\min\mu_{h},\max\mu_{h}]=\emptyset$..
However, the expression of $\mathbb{L}$ and $\mathbb{S}$ in terms of the
Hadamard product can be useful to this end. Indeed, another important property
of the Hadamard product is that, for any matrices ${\bf A}$ and ${\bf B}$,
$\text{rank}({\bf A}\circ{\bf B})\leq\text{rank}({\bf A})\text{rank}({\bf
B})$. Therefore,

$\displaystyle\text{rank}(\mathbb{L})$
$\displaystyle\leq\text{rank}(\mathcal{C})\cdot\text{rank}({\bf V}{\bf R}-{\bf
L}{\bf W})\leq(p+q)\cdot\text{rank}(\mathcal{C}),$

and similarly for $\mathbb{S}$. Thus, we have

$\displaystyle\text{rank}(\mathbb{S}-x\mathbb{L})$ $\displaystyle\leq
2(p+q)\cdot\text{rank}(\mathcal{C}),\quad\forall\,x\in\\{f_{i}\\}.$ (45)

In general, the Cauchy matrix $\mathcal{C}$ is full rank so this inequality is
trivially satisfied. However, depending on the partitioning of the points into
$\lambda_{k}$ and $\mu_{h}$ (as in (4) and (6)), it can be _numerically_ low-
rank (see, e.g., (Pan2014, , Theorem 5), Beckermann2019 ; Chandrasekaran2007
). If $\pi_{\mathcal{C}}$ denotes the numerical rank of $\mathcal{C}$, then
$2(p+q)\pi_{\mathcal{C}}$ is a rough estimate for the numerical rank of
$\mathbb{S}-x\mathbb{L}$333For $\mathbb{L}_{r}$ and $\mathbb{S}_{r}$
satisfying (37) and (38), the value $4(p+q)\pi_{\mathcal{C}}$ can be used as
an estimate for the number of singular triplets to compute.. Oftentimes, the
underlying dynamical system is of much lower complexity, thus allowing for the
computation of a minimal realization of reduced order $n$. One can also use
insight of the system itself or count the number of peaks in the frequency
response to estimate $n$ (for systems with poles having dominant imaginary
parts).

###### Remark 4

The expression of $\mathbb{L}$ and $\mathbb{S}$ in terms of the Hadamard
product provides us with an upper bound of the spectral norm of the Loewner
and shifted Loewner matrix. Indeed, the spectral norm is submultiplicative
with respect to the Hadamard product (Horn1991, , Theorem 5.5.1), hence

$\begin{array}[]{rll}\|\mathbb{L}\|&=&\|\mathcal{C}\circ({\bf V}{\bf R}-{\bf
L}{\bf W})\|\leq\|\mathcal{C}\|\cdot\|{\bf V}{\bf R}-{\bf L}{\bf W}\|\\\
&\leq&\|\mathcal{C}\|_{F}\cdot\|{\bf V}{\bf R}-{\bf L}{\bf W}\|=\|{\bf V}{\bf
R}-{\bf L}{\bf
W}\|\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{1}{\mu_{i}-\lambda_{j}}\right|^{2}},\\\
\end{array}$

where $\|\mathcal{C}\|_{F}$ denotes the Frobenius norm of $\mathcal{C}$. Note
that $\|{\bf V}{\bf R}-{\bf L}{\bf W}\|$ can be computed cheaply, e.g., by a
power method exploiting the low rank of ${\bf V}{\bf R}-{\bf L}{\bf W}$.

Similarly,

$\|\mathbb{S}\|\leq\|{\bf M}{\bf V}{\bf R}-{\bf L}{\bf
W}\bm{\Lambda}\|\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{1}{\mu_{i}-\lambda_{j}}\right|^{2}}.$

###### Remark 5

Low-rank approximations to $\mathbb{L}$ and $\mathbb{S}$ may be computed by
adaptive cross approximation ACA , particularly suited for hierarchical
matrices, the CUR decomposition CUR  as in LoewCUR ; 9073015 , or related
schemes. These approaches select a certain number of columns and rows of the
original matrices in a greedy fashion based on various heuristics, and a
_core_ matrix is utilised to compute a low-rank approximation. If a given
threshold on the desired accuracy of the computed approximation is provided as
an input, these algorithms often construct matrices whose rank is much larger
than the one of the target matrices $\mathbb{L}$ and $\mathbb{S}$. On the
other hand, by fixing the rank $k$ of the approximation,
$k\approx\mathtt{rank}(\mathbb{L}),\mathtt{rank}(\mathbb{S})$ \- assuming we
know an estimate of $\mathtt{rank}(\mathbb{L})$, $\mathtt{rank}(\mathbb{S})$
\- the accuracy we achieve may be very low affecting the reliability of the
computed reduced models.

###  3.2 Hierarchically semiseparable (HSS) representation of a Cauchy matrix

The literature on HSS matrices is rather vast and technical (see, e.g.,
Vandebril2005 ; Pan2014 ; PAN2015107 ; Chandrasekaran2006 ; Xia2010  and
references therein). Here we recall only the main properties of this class of
matrices and their role in the efficient representation of Cauchy matrices.
Such a technique is also closely related to the Fast Multipole Method (FMM).
We refer the interested reader to, e.g., Chandrasekaran2007 ;
Chandrasekaran2006  for more details on the interconnection between HSS
matrices and FMM.

######  Definition 1 ((PAN2015107, , Definition 27))

Let ${\bf A}$ be an $N\times N$ matrix with $\alpha$ being the maximum rank of
all its subdiagonal blocks, namely the blocks of all sizes lying strictly
below the block diagonal, and $\beta$ the maximum rank of all its
superdiagonal blocks, namely the blocks of all sizes lying strictly above the
block diagonal, respectively. Then, ${\bf A}$ is $(\alpha,\beta)$-HSS if its
diagonal blocks consist of $\mathcal{O}((\alpha+\beta)N)$ entries.

The $(\alpha,\beta)$-HSS representation of a matrix ${\bf A}$ is very
advantageous whenever $\alpha$ and $\beta$ are small. For instance, it allows
us to express ${\bf A}$ in terms of $\mathcal{O}((\alpha+\beta)N)$ parameters
avoiding storing its $N^{2}$ entries. Moreover, a whole, efficient HSS
arithmetic has been developed in the last decades (see, e.g.,
Chandrasekaran2006 ; Xia2010 ). For instance, the computational cost of the
matrix-vector product ${\bf A}{\bf y}$ amounts to
$\mathcal{O}((\alpha+\beta)N)$ FLOPs. If ${\bf A}$ is nonsingular, its inverse
is also a $(\alpha,\beta)$-HSS matrix that can be computed in
$\mathcal{O}((\alpha+\beta)^{3}N)$ FLOPs (see, e.g., (PAN2015107, , Section
6)).

To fully exploit the HSS framework for our purposes, we wish to represent the
Cauchy matrix $\mathcal{C}$ in terms of a HSS matrix with a low off-diagonal
rank. In light of Corollary 1, this would considerably decrease the
computational cost of the matrix-vector products involving
$\mathbb{S}-x\mathbb{L}$ while avoiding forming the dense matrices
$\mathbb{S}$ and $\mathbb{L}$.

The construction of an HSS approximation $\widetilde{\mathcal{C}}$ to
$\mathcal{C}$ is rather involved and the magnitude of the
$(\alpha,\beta)$-rank of the computed $\widetilde{\mathcal{C}}$ strictly
depends on the partitioning of the frequencies along with the accuracy that
has been selected for the actual computation of
$\widetilde{\mathcal{C}}$444Roughly speaking, such a threshold is related to
the computation of the low-rank approximations to the off-diagonal blocks of
$\mathcal{C}$ (see, e.g., (Xi2014, , Corollary 4.3), (Kressner2019, , Theorem
4.7)). (see, e.g., (PAN2015107, , Section 8) for further details on the
computation of an HSS-representation of a Cauchy matrix). In this paper we
employ the readily available hm-toolbox hmtoolbox .

###### Example 1

We investigate the impact of the most commonly-used frequency partitions
(Half&Half, Odd&Even, Odd&Even (Real)) on the HSS-rank of the computed
$\widetilde{\mathcal{C}}$ for a mechanical structure. We emphasize that the
most effective partition is problem-dependent and is still an open problem,
beyond the scope of this paper. We consider the _Flexible Aircraft_ data set
POUSSOTVASSAL2018559  from the MORwiki morWiki . This dataset contains 421
frequency values $\omega_{j}$ expressed in rad/s and the corresponding
measurements of the transfer function ${\bf H}_{j}$. We disregard the last
data point and consider the remaining frequencies ranging from $f_{1}=0.1$Hz
to $f_{420}=42$Hz. As this is a mechanical structure, frequencies considered
are in the low spectrum, as opposed to electrical systems, for which
frequencies typically span the GHz range.

We recall the three different partitions of the frequencies
$\\{\omega_{j}\\}_{j=1}^{j=420}$:

  * •

Half&Half:
$\bm{\Lambda}=\text{diag}([\mathrm{i}\omega_{1},-\mathrm{i}\omega_{1},\ldots,\mathrm{i}\omega_{210},-\mathrm{i}\omega_{210}])$,
${\bf
M}=\text{diag}([\mathrm{i}\omega_{211},-\mathrm{i}\omega_{211},\ldots,\mathrm{i}\omega_{420},-\mathrm{i}\omega_{420}])$.

  * •

Odd&Even:
$\bm{\Lambda}=\text{diag}([\mathrm{i}\omega_{1},-\mathrm{i}\omega_{1},\ldots,\mathrm{i}\omega_{419},-\mathrm{i}\omega_{419}])$,
${\bf
M}=\text{diag}([\mathrm{i}\omega_{2},-\mathrm{i}\omega_{2},\ldots,\mathrm{i}\omega_{420},-\mathrm{i}\omega_{420}])$.

  * •

Odd&Even (Real):
$\bm{\Lambda}_{r}=\text{diag}([-\omega_{1}^{2},-\omega_{1}^{2},\ldots,-\omega_{419}^{2},-\omega_{419}^{2}])$,
${\bf
M}_{r}=\text{diag}([-\omega_{2}^{2},-\omega_{2}^{2},\ldots,-\omega_{420}^{2},-\omega_{420}^{2}])$.

| Half&Half | Odd&Even | Odd&Even (Real)
---|---|---|---
$\mathtt{hssrank}(\widetilde{\mathcal{C}})$ | 32 | 30 | 13
$\mathtt{rank}(\mathcal{C})$ | 36 | 420 | 210
$\|\widetilde{\mathcal{C}}-\mathcal{C}\|/\|\mathcal{C}\|$ | 2.68e-12 | 2.61e-11 | 6.62e-13
Table 1: Example 1. HSS-rank and relative error of the HSS representation
$\widetilde{\mathcal{C}}$ of $\mathcal{C}$ for different frequency partitions
along with the rank of $\mathcal{C}$.

For each partition, we compute the corresponding Cauchy matrix in HSS format
$\widetilde{\mathcal{C}}$ without assembling the full $\mathcal{C}$
beforehand, by means of the function hss of the hm-toolbox:

$\widetilde{\mathcal{C}}=\mathtt{hss(^{\prime}cauchy^{\prime},dM,-dL,N,N)}$

where $\mathtt{dM}$ and $\mathtt{dL}$ are $N$ dimensional vectors containing
the frequencies $\mu_{h}$ and $\lambda_{k}$, respectively. We then calculate
its rank by $\mathtt{hssrank}(\widetilde{\mathcal{C}})$555Following Definition
1, this function returns $\max\\{\alpha,\beta\\}$.. In Table 1 we report the
HSS-rank of the matrix $\widetilde{\mathcal{C}}$ for the partitions mentioned
above. Thanks to the small dimension of the dataset, we are able to compute
the full Cauchy matrix $\mathcal{C}$ and document its (standard) rank along
with the relative error
$\|\widetilde{\mathcal{C}}-\mathcal{C}\|/\|\mathcal{C}\|$. As expected, having
two disjoint sets of frequencies like in the Half&Half partition leads to a
Cauchy matrix $\mathcal{C}$ whose (standard) rank is low. This does not happen
in the other two scenarios we examine so that taking advantage of the HSS
format is necessary to achieve memory-saving representations of $\mathcal{C}$.
The results in Table 1 show that a good accuracy in terms of the relative
error can be achieved for all three frequency partitions. Nevertheless, the
HSS rank of $\widetilde{\mathcal{C}}$ is significantly lower for the Odd&Even
(Real) partition, most likely due to the squaring of the frequencies performed
in Odd&Even (Real), which leads to a fast decay in the magnitude of the off-
diagonal entries of $\mathcal{C}$. Hence, for a fixed threshold, the off-
diagonal blocks of the Cauchy matrix associated to the Odd&Even (Real)
partition can be approximated by matrices having a smaller rank than those
associated to the other two scenarios we examined.

In Figure 1 we display the absolute value – on a logarithmic scale – of the
entries of the Cauchy matrix $\mathcal{C}$ stemming from the different
partitions. The same scale has been used in all the three figures, enforcing
the observation that the Odd&Even (Real) partition exhibits the fastest decay
in the magnitude of the off-diagonal entries of $\mathcal{C}$ .

Half&Half

Odd&Even

Odd&Even (real)

Figure 1: Example 1. Absolute value – on a logarithmic scale – of the entries
of the Cauchy matrix $\mathcal{C}$ stemming from the three different frequency
partitions we have examined.

###  3.3 Efficient, inexact matrix-vector products

Whenever the matrix $\mathcal{C}$ admits an accurate approximation in terms of
a low-rank HSS matrix $\widetilde{\mathcal{C}}$, the computational cost of
performing the matrix-vector product $(\mathbb{S}-x\mathbb{L}){\bf y}$ can be
significantly reduced.

###### Proposition 1

Let $\widetilde{\mathcal{C}}$ be an $(\alpha,\beta)$-HSS matrix that
approximates the Cauchy matrix $\mathcal{C}$ accurately. If $\mathbb{L}$ and
$\mathbb{S}$ satisfy the Sylvester equations in (29), then

$(\mathbb{S}-x\mathbb{L}){\bf y}=\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf
y}+\mathcal{\bm{E}}\widetilde{\bf y},$ (46)

where
$\|\mathcal{\bm{E}}\|\leq(p+q)\|\mathcal{C}-\widetilde{\mathcal{C}}\|\max^{2}\\{\|{\bf
V}\|_{\infty},\|{\bf R}\|_{\infty},\|{\bf L}\|_{\infty},\|{\bf
W}\|_{\infty}\\}$. Moreover, the computational cost of performing

$\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf y},$ (47)

amounts to $\mathcal{O}((p+q)(\alpha+\beta+1)N)$ FLOPs.

###### Proof

From the result in Corollary 1, we can write

$\displaystyle(\mathbb{S}-x\mathbb{L}){\bf y}=$
$\displaystyle\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf y}$
$\displaystyle+\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ((\mathcal{C}-\widetilde{\mathcal{C}})(\widetilde{\bf
r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ((\mathcal{C}-\widetilde{\mathcal{C}})(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\widetilde{\bf
v}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
r}_{j}^{*}\circ\widetilde{\bf
y}))-\sum_{j=1}^{p}\widetilde{\ell}_{j}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
w}_{j}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf
y}+\mathcal{\bm{E}}\widetilde{\bf y},$

where $\mathcal{\bm{E}}:=\displaystyle\sum_{j=1}^{q}\text{diag}(\widetilde{\bf
v}_{j})(\mathcal{C}-\widetilde{\mathcal{C}})\text{diag}(\widetilde{\bf
r}_{j}^{*})-\sum_{j=1}^{p}\text{diag}(\widetilde{\ell}_{j})(\mathcal{C}-\widetilde{\mathcal{C}})\text{diag}(\widetilde{\bf
w}_{j}^{*})$. Therefore,

$\displaystyle\|\mathcal{\bm{E}}\|\leq$
$\displaystyle\|\mathcal{C}-\widetilde{\mathcal{C}}\|\left(\sum_{j=1}^{q}\|\text{diag}(\widetilde{\bf
v}_{j})\|\|\text{diag}(\widetilde{\bf
r}_{j}^{*})\|+\sum_{j=1}^{p}\|\text{diag}(\widetilde{\ell}_{j})\|\|\text{diag}(\widetilde{\bf
w}_{j}^{*})\|\right)$ $\displaystyle=$
$\displaystyle\|\mathcal{C}-\widetilde{\mathcal{C}}\|\left(\sum_{j=1}^{q}\|\widetilde{\bf
v}_{j}\|_{\infty}\|\widetilde{\bf
r}_{j}\|_{\infty}+\sum_{j=1}^{p}\|\widetilde{\ell}_{j}\|_{\infty}\|\widetilde{\bf
w}_{j}\|_{\infty}\right)$ $\displaystyle\leq$
$\displaystyle(p+q)\|\mathcal{C}-\widetilde{\mathcal{C}}\|\max{{}^{2}}\\{\|{\bf
V}\|_{\infty},\|{\bf R}\|_{\infty},\|{\bf L}\|_{\infty},\|{\bf
W}\|_{\infty}\\}.$

This proves the first part of Proposition 1. To conclude, by making use of the
property that the matrix-vector product with a $(\alpha,\beta)$-HSS matrix
costs $\mathcal{O}((\alpha+\beta)N)$ FLOPs and that ${\bf V}{\bf R}$ has rank
$q$, a direct computation shows that the number of operations needed to
perform (47) amounts to $\mathcal{O}((p+q)(\alpha+\beta+1)N)$ FLOPs, which
proves the second claim in Proposition 1.

As before, analogous results can be obtained for $\mathbb{L}_{r}$ and
$\mathbb{S}_{r}$ satisfying (37) and (38), respectively.

Proposition 1 shows that, whenever $\|\mathcal{C}-\widetilde{\mathcal{C}}\|$
is small, the matrix-vector product $(\mathbb{S}-x\mathbb{L}){\bf y}$ can be
well-approximated by the expression in (47) while dramatically reducing the
computational complexity from $\mathcal{O}(N^{2})$ FLOPs to
$\mathcal{O}((p+q)(\alpha+\beta+1)N)$ FLOPs. However, when this approximation
is used within our favorite iterative procedure for computing a partial SVD of
$\mathbb{S}-x\mathbb{L}$, the inexactness introduced by neglecting the term
$\mathcal{\bm{E}}\widetilde{\bf y}$ should be taken into account.

The use of inexact matrix-vector products within certain iterative procedures
has been the subject of numerous research papers: Krylov techniques for
solving linear systems and matrix equations Simoncini2003 ; Bouras2005 ;
Eshof2004 ; Kuerschner2019 ; Kuerschner2018 , eigenvalue problems Freitag2007
; Simoncini2002 , or an inexact variant of the Lanczos bidiagonalization for
the computation of some leading singular triplets of a generic matrix function
$f({\bf A})$ can be found in GAAF2017 . With the goal to decrease the
computational cost of the overall procedure, these studies show that the
accuracy of the matrix-vector product can be _relaxed_ (becoming more and more
inaccurate) as iterations proceed. In our framework, the inexactness
introduced by approximating $(\mathbb{S}-x\mathbb{L}){\bf y}$ with (47) is
fixed throughout the entire iterative procedure and mainly depends on
$\|\mathcal{C}-\widetilde{\mathcal{C}}\|$, which is often small, as shown in
Example 1. Therefore, the approximation

$(\mathbb{S}-x\mathbb{L}){\bf y}\approx\sum_{i=1}^{q}\widetilde{\bf
v}_{i}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
r}_{i}^{*}\circ\widetilde{\bf
y}))-\sum_{i=1}^{p}\widetilde{\ell}_{i}\circ(\widetilde{\mathcal{C}}(\widetilde{\bf
w}_{i}^{*}\circ\widetilde{\bf y}))+{\bf V}{\bf R}{\bf y},$

does not greatly affect the accuracy of the computed singular triplets (see
section 4). Moreover, in our case, we do not need an accurate approximation of
the singular triplets of $\mathbb{S}-x\mathbb{L}$. The main goal is to have
meaningful spaces spanned by the computed left and right singular vectors so
that the obtained reduced model $({\bf E},{\bf A},{\bf B},{\bf C})$ inherits
the desired approximation properties. Moreover, as shown in (IonitaPhD, ,
Corollary 1.4), (TutorialLoewner, , Proposition 8.25), in the case of noise-
free measurements of a low-order rational function, even general projectors,
not necessarily obtained from the SVD, can be employed for identifying the
underlying function.

###### Remark 6

In Remark 3 we suggested to use the value $2(p+q)\pi_{C}$, where $\pi_{C}$ is
the numerical rank of $\mathcal{C}$, to decide on the number $n$ of singular
triplets of $\mathbb{S}-x\mathbb{L}$ needed for the reduced model. For
interlaced partitions, as it is the case with Odd&Even and Odd&Even (real)
(see Table 1), the numerical (standard) rank of the Cauchy matrix is large, in
general. Hence, the value $n=2(p+q)\max\\{\alpha,\beta\\}$ may instead be
employed for the computation of a meaningful reduced model whenever
$\mathcal{C}$ can be well-approximated by a $(\alpha,\beta)$-HSS matrix
$\widetilde{\mathcal{C}}$666As before, the value
$4(p+q)\max\\{\alpha,\beta\\}$ should be preferred whenever $\mathbb{L}_{r}$
and $\mathbb{S}_{r}$ solve (37) and (38), respectively.. Moreover, the HSS-
rank of $\widetilde{\mathcal{C}}$ is obtained as a byproduct of the
construction of $\widetilde{\mathcal{C}}$.

###### Remark 7

If $\mathcal{C}$ admits an accurate approximation in terms of an
$(\alpha,\beta)$-HSS matrix $\widetilde{\mathcal{C}}$, the expression in
Theorem 1 shows that $\mathbb{L}$ can also be well-approximated by a HSS
matrix $\widetilde{\mathbb{L}}$ whose rank is at most
$((p+q)\alpha,(p+q)\beta)$. Even though the computational cost of
$\widetilde{\mathbb{L}}{\bf y}$ would still be
$\mathcal{O}((p+q)(\alpha+\beta)N)$ FLOPs, using the HSS approximation
$\widetilde{\mathbb{L}}$ of $\mathbb{L}$ may be very advantageous whenever
linear systems with $\mathbb{L}$ need to be solved (see, e.g., the procedure
presented in Embree2019  for the pseudospectra computation of
$\mathbb{S}-x\mathbb{L}$). Indeed, as mentioned in section 3.2, the
computation of the inverse $\widetilde{\mathbb{L}}^{-1}$ of
$\widetilde{\mathbb{L}}$ costs $\mathcal{O}((p+q)^{3}(\alpha+\beta)^{3}N)$
FLOPs. Once $\widetilde{\mathbb{L}}^{-1}$ is computed, we need only
$\mathcal{O}((p+q)(\alpha+\beta)N)$ FLOPs to perform
$\widetilde{\mathbb{L}}^{-1}{\bf y}$.

##  4 Numerical results

In this section we present numerical experiments illustrating the potential of
the proposed approach.

In Example 2, we compare our approach to standard procedures employed in the
Loewner framework. Recall that the main steps in the standard approach involve
forming the full Loewner and shifted Loewner matrices $\mathbb{L}$ and
$\mathbb{S}$ and computing the SVD of $\mathbb{S}-x\mathbb{L}$. This SVD can
be either computed in full, followed by keeping only the $n$ dominant singular
vectors, or only these $n$ singular vectors can be obtained by means of an
iterative procedure, where the matrix-vector product with
$\mathbb{S}-x\mathbb{L}$ is needed777We employ the Matlab functions svd and
svds, respectively.. In the following, we report the overall running time,
considering the _construction_ step (Construction), i.e., the computation of
$\mathbb{L}$ and $\mathbb{S}$ in the standard approach and of
$\widetilde{\mathcal{C}}$ in our approach, as well as the _reduction_ step
(Reduction), involving the SVD computation followed by projection to obtain
the reduced matrices in (41). In terms of memory requirements, for our
approach, this involves the allocation of $\widetilde{\mathcal{C}}$ in the HSS
format, while for the standard approach, we report the storage required for
$\mathbb{L}$ and $\mathbb{S}$.

In Table 2 we recall the computational cost of the construction and reduction
steps of both the standard approach, based on either a full or a partial SVD,
and the novel one presented in this paper along with their memory
requirements.

| Construction | Reduction | Storage
---|---|---|---
Full svd | $\mathcal{O}(N^{2}(p+q))$ | $\mathcal{O}(N^{3})$ | $\mathcal{O}(N^{2})$
svds w/ $\mathbb{S}-x\mathbb{L}$ | $\mathcal{O}(N^{2}(p+q))$ | $\mathcal{O}(nN^{2})$ | $\mathcal{O}(N^{2})$
svds w/ $\widetilde{\mathcal{C}}$ | $\mathcal{O}(\max\\{\alpha,\beta\\}N\log N)$ | $\mathcal{O}(n(p+q)(\alpha+\beta+1)N)$ | $\mathcal{O}((\alpha+\beta+p+q)N)$
Table 2: Computational cost of the construction (Construction) and reduction
(Reduction) steps of the different approaches we test along with their storage
demand (Storage). The computational cost of the construction of
$\widetilde{\mathcal{C}}$ can be found, e.g., in (hmtoolbox, , Table 1).

Lastly, the accuracy of the reduced models is reported in terms of the
normalized $\mathcal{H}_{2}$-error:

$\mathcal{H}_{2}-\text{error}=\sqrt{\frac{\sum_{j=1}^{N}\|{\bf H}_{j}-{\bf
H}(i\omega_{j})\|_{F}^{2}}{\sum_{j=1}^{N}\|{\bf H}_{j}\|_{F}^{2}}},$

where $\|\cdot\|_{F}$ denotes the Frobenius norm. Similar results in terms of
accuracy are attained for the $\mathcal{H}_{\infty}$-error, however, we
decided not to document them here, for the sake of brevity.

In Example 3, we compare our novel strategy to the one presented in FSVDL ,
which makes use of the low-rank ADI-Galerkin method for computing the Loewner
matrix as the solution to (29). Such a scheme computes low-rank approximations
to the dense Loewner matrix to speed-up the SVD computation, however, the
memory constraints originating from the allocation of $\mathbb{L}$ and
$\mathbb{S}$ are still present.

Results were obtained by running Matlab R2020b MATLAB2020b  on a MacBook Pro
with an Intel Core i9 processor running at 2.3GHz using 16GB of RAM. All
computations involving HSS matrices employed the hm-toolbox hmtoolbox  with
the default settings and the threshold for off-diagonal truncation set to
$10^{-14}$.

###### Example 2

We consider a synthetic problem for which we can control the order of the
original system ($n$), the number of inputs and outputs ($p=q$), as well as
the number of measurements ($N$). The system dynamics is generated randomly,
with poles in complex conjugate pairs. In particular:

  * •

the real part of the poles is random with mean $-10^{4}$ and standard
deviation $-2\cdot 10^{3}$; the imaginary part is also random, with mean
$10^{4}$ and standard deviation $10^{6}$.

  * •

residues associated to each pole are rank-1 matrices, obtained as outer
products between two random vectors, both having the real part with mean $0$
and standard deviation $10$, while the imaginary part has mean $0$ and
standard deviation $10^{2}$.

Measurement points $\\{\omega_{j}\\}_{j=1}^{j=N}$ are logarithmically
distributed between $10^{4}$ and $10^{7}$ rad/sec. Last, but not least, random
noise with a signal-to-noise ratio $SNR=100$ was added to the transfer
function evaluation ${\bf H}(\mathrm{i}\omega_{j})$ to obtain the measurement
matrices ${\bf H}_{j}$. We adopt the Odd&Even (real) partition of the
frequencies as it achieves satisfactory approximation results while
eliminating complex arithmetic.

We compare the proposed approach to the traditional Loewner framework, in
which the Loewner and shifted Loewner matrices $\mathbb{L}$ and $\mathbb{S}$
are formed and the full SVD of $\mathbb{S}-x\mathbb{L}$ is computed, as well
as the alternative approach in which, after building $\mathbb{L}$ and
$\mathbb{S}$, a partial SVD of $\mathbb{S}-x\mathbb{L}$ using the Matlab svds
function is computed for various instances of the data set described above for
different values of $N$, $p$, and $n$. The command svds was employed with the
left starting vector $\widetilde{\bf v}_{1}$ (same notation as in Theorem 1)
instead of a random starting vector, which is the default setting.

Figure 2 presents the memory requirements for storing the Loewner and shifted
Loewner matrices $\mathbb{L}$ and $\mathbb{S}$ (in red), as opposed to storing
the HSS approximation $\widetilde{\mathcal{C}}$ in our approach (in blue),
along with the storage needed to allocate the data in $\bm{\Lambda}_{r}$,
$\bm{M}_{r}$, ${\bf R}_{r}$, ${\bf L}_{r}$, ${\bf V}_{r}$, ${\bf W}_{r}$, for
increasing values in the number of inputs and outputs $p$ (in black). We point
out that for values of $N$ larger than $40\,000$, we were not able to allocate
the full matrices $\mathbb{L}$ and $\mathbb{S}$ on the employed laptop (this
value, however, depends on the available RAM memory of the machine). For
instances when these matrices can be allocated, Figure 2 shows that the memory
requirements for the proposed approach are always much lower than for the
standard scheme. Moreover, in contrast to what happens to the memory required
for the data matrices, the storage demanded by the allocation of
$\widetilde{\mathcal{C}}$ in HSS format is independent of $p$.

Figure 2: Example 2. Memory requirements in Megabytes to store $\mathbb{L}$,
$\mathbb{S}$, $\widetilde{\mathcal{C}}$, and the data matrices ($\Lambda_{r}$,
$M_{r}$, ${\bf V}_{r}$, ${\bf W}_{r}$, ${\bf L}_{r}$, ${\bf R}_{r}$) for
different values of $N$ and $p$.

We report the results of the comparison between the different approaches in
terms of run time in Table 3 for the number of measurements $N$ varying
between $1\,000$ to $100\,000$, the number of inputs and outputs $p$ taking
values $1$, $5$ and $10$, and the number of poles being $50$ or $100$. The “–”
is used to indicate the instances for which we were not able to compute the
reduced model (41): for $N>50\,000$, we cannot allocate the full matrices
$\mathbb{L}$ and $\mathbb{S}$, and for $N=30\,000,40\,000$ we could not
compute the full SVD of $\mathbb{S}-x\mathbb{L}$ . Such constraints are not
relevant to our proposed strategy. It is pertinent to remark the following:

  1. 1.

the CPU time of the full SVD approach does not depend on $p$ and $n$, only on
$N$, as expected from Table 2: indeed, the cost of building $\mathbb{L}$ and
$\mathbb{S}$ is quadratic in $N$ whereas the full SVD demands
$\mathcal{O}(N^{3})$ FLOPs; the full SVD approach is rarely the fastest method
(it can happen for very modest values of $N$ in the considered range);

  2. 2.

the CPU time of the full assembly of $\mathbb{S}-x\mathbb{L}$ followed by the
svds Matlab command does not depend on $p$, only on $N$ and $n$, as expected
from Table 2: the construction of $\mathbb{L}$ and $\mathbb{S}$ costs
$\mathcal{O}(N^{2})$ FLOPs, whereas the computational effort for the partial
SVD depends on $n$, leading to a more demanding procedure for large $n$; it is
usually the fastest approach for (very) modest values of $N$ in the considered
range and $p>1$;

  3. 3.

the HSS rank of the Cauchy matrix approximation $\widetilde{\mathcal{C}}$ only
depends on the frequency samples, hence on $N$ because, in our scenario, the
sampling interval is the same, but the distribution of points inside the
interval is different for each $N$; there may be instances when, for the same
samples, the HSS rank of $\widetilde{\mathcal{C}}$ may produce slightly
different results due to the randomness induced by the adaptive cross
approximation procedure used in constructing $\widetilde{\mathcal{C}}$ (for
instance, for $N=50\,000$, $p=5$, $n=50$ and $n=100$, the rank is $28$, while
for the rest of the values considered for $n$ and $p$, the rank is $27$);
moreover, the HSS rank increases with $N$;

  4. 4.

our proposed approach is as accurate as the first two approaches, highlighting
the fact that the HSS approximation $\widetilde{\mathcal{C}}$ does not lead to
significant losses in the approximation properties of the reduced model (41);
clearly, our approach cannot be more accurate than the traditional Loewner
framework, especially when the full SVD is performed;

  5. 5.

last, but not least, the CPU time of the proposed solution depends linearly on
$p$, $n$ and as $N\log N$ (Table 2), thus being the fastest method for large
values of $N$; moreover, no memory constraints are present for $N$ up to
$100\,000$.

|  |  | Full svd |  svds w/ $\mathbb{S}-x\mathbb{L}$ |  svds w/ $\widetilde{\mathcal{C}}$
---|---|---|---|---|---
$N$ | $p$ | $n$ | Time (s) |  $\mathcal{H}_{2}$-error | Time (s) |  $\mathcal{H}_{2}$-error | $\mathtt{hssrank}(\widetilde{\mathcal{C}})$ | Time (s) |  $\mathcal{H}_{2}$-error
1 000 | 1 | 50 | 0.28 | 3.62e-10 | 0.20 | 3.62e-10 | 15 | 1.07 | 3.62e-10
3 000 | 10.69 | 3.71e-10 | 3.51 | 3.71e-10 | 19 | 3.16 | 3.71e-10
5 000 | 20.79 | 3.7e-10 | 12.29 | 3.7e-10 | 21 | 6.99 | 3.7e-10
10 000 | 158.41 | 3.71e-10 | 68.40 | 3.71e-10 | 22 | 15.55 | 3.71e-10
15 000 | 554.98 | 3.73e-10 | 209.31 | 3.73e-10 | 24 | 22.86 | 3.73e-10
29 000 | 4674.35 | 3.74e-10 | 1590.37 | 3.74e-10 | 26 | 52.33 | 3.74e-10
30 000 | – | – | 1827.45 | 3.74e-10 | 26 | 50.68 | 3.74e-10
40 000 | – | – | 11214.41 | 3.74e-10 | 27 | 71.44 | 3.74e-10
50 000 | – | – | – | – | 27 | 91.88 | 3.74e-10
100 000 | – | – | – | – | 30 | 189.71 | 3.75e-10
1 000 | 1 | 100 | 0.21 | 9.63e-11 | 0.22 | 9.63e-11 | 15 | 1.59 | 9.64e-11
3 000 | 10.66 | 1.01e-10 | 5.84 | 1.01e-10 | 19 | 5.65 | 1.01e-10
5 000 | 20.62 | 1.01e-10 | 18.89 | 1.01e-10 | 21 | 12.77 | 1.01e-10
10 000 | 156.76 | 1.01e-10 | 93.16 | 1.01e-10 | 22 | 28.47 | 1.02e-10
50 000 | – | – | – | – | 27 | 155.46 | 1.03e-10
100 000 | – | – | – | – | 30 | 321.84 | 1.03e-10
1 000 | 5 | 50 | 0.27 | 3.71e-10 | 0.19 | 3.71e-10 | 15 | 1.20 | 3.72e-10
3 000 | 10.69 | 3.17e-10 | 3.56 | 3.17e-10 | 19 | 4.59 | 3.17e-10
5 000 | 20.81 | 3.01e-10 | 12.29 | 3.01e-10 | 21 | 9.51 | 3.02e-10
10 000 | 157.75 | 2.92e-10 | 68.43 | 2.92e-10 | 22 | 19.67 | 2.92e-10
50 000 | – | – | – | – | 28 | 107.74 | 2.88e-10
100 000 | – | – | – | – | 30 | 230.17 | 2.88e-10
1 000 | 5 | 100 | 0.26 | 2.52e-10 | 0.25 | 2.52e-10 | 15 | 2.21 | 2.52e-10
3 000 | 10.69 | 1.41e-10 | 5.90 | 1.41e-10 | 19 | 8.62 | 1.41e-10
5 000 | 20.78 | 1.35e-10 | 18.97 | 1.35e-10 | 21 | 17.73 | 1.35e-10
10 000 | 157.75 | 1.33e-10 | 93.58 | 1.33e-10 | 22 | 36.84 | 1.33e-10
50 000 | – | – | – | – | 28 | 197.72 | 1.31e-10
100 000 | – | – | – | – | 30 | 421.02 | 1.3e-10
1 000 | 10 | 50 | 0.26 | 6.57e-10 | 0.18 | 6.57e-10 | 15 | 1.61 | 6.57e-10
3 000 | 11.05 | 3.2e-10 | 3.67 | 3.2e-10 | 19 | 5.58 | 3.2e-10
5 000 | 20.83 | 2.73e-10 | 12.25 | 2.73e-10 | 21 | 10.94 | 2.73e-10
10 000 | 159.13 | 2.68e-10 | 69.34 | 2.68e-10 | 22 | 23.09 | 2.68e-10
50 000 | – | – | – | – | 27 | 132.64 | 2.56e-10
100 000 | – | – | – | – | 30 | 293.55 | 2.54e-10
1 000 | 10 | 100 | 0.24 | 5.27e-10 | 0.24 | 5.27e-10 | 15 | 2.88 | 5.25e-10
3 000 | 10.68 | 1.78e-10 | 5.94 | 1.78e-10 | 19 | 10.45 | 1.78e-10
5 000 | 20.84 | 1.73e-10 | 18.97 | 1.73e-10 | 21 | 20.50 | 1.73e-10
10 000 | 157.48 | 1.65e-10 | 93.63 | 1.65e-10 | 22 | 42.97 | 1.65e-10
50 000 | – | – | – | – | 27 | 248.91 | 1.58e-10
100 000 | – | – | – | – | 30 | 552.66 | 1.58e-10
Table 3: Example 2. Computational time (in seconds) and
$\mathcal{H}_{2}$-error achieved by each approach for different values of $N$
(number of samples), $p$ (number of inputs and outputs), and $n$ (order of the
underlying system and of the model) on the employed laptop.

Figure 3: Example 2. Left: Frequency response obtained with the standard
approach (in black) and the measurements (in red) for $p=1$, $n=50$, and
$N=10\,000$. Right: Frequency response obtained with the proposed approach (in
black) and the measurements (in red) for $p=1$, $n=50$, and $N=10\,000$.

In Figure 4 (left) we plot the computational time of the three approaches for
$p=1$, $n=50$, and different values of $N$. Even though these are the same
results as those reported in Table 3, Figure 4 (left) clearly shows the
$\mathcal{O}(N^{3})$ trend of the full SVD scheme versus the
$\mathcal{O}(N^{2})$ trend of the svds scheme versus the $\mathcal{O}(N)$
behaviour of the proposed approach. In Figure 4 (right) we depict, on a
logarithmic scale, the running time of the proposed procedure for $n=50$ and
different values of $N$ and $p$, clearly exhibiting a linear dependency on $p$
and an $N\log N$ dependency with respect to $N$.

Figure 4: Example 2. Left: Computational time achieved by the different
approaches for $p=1$, $n=50$, and $N$. Right: Computational time achieved by
our novel procedure for $n=50$, and different values of $N$ and $p$.

###### Example 3

In this example we compare the novel strategy presented in this paper to the
fast Loewner SVD scheme illustrated in FSVDL . We consider the same data set
as the one in Example 2, this time with $SNR=120$ and a random ${\bf
D}\in\mathbb{R}^{p\times p}\neq{\bf 0}$. Due to the fact that the models
resulting from the Loewner framework have ${\bf D}={\bf 0}$, a realization of
size $n+p$ is needed to approximate the system with ${\bf D}\neq{\bf 0}$
artAJMACA ; SLACATCAD09 .

In FSVDL , a Galerkin-ADI method is applied to the Sylvester equation (29)
satisfied by the Loewner matrix. At the $k$-th iteration, a low-rank
approximation $P_{k}L_{k}Q_{k}^{*}$,
$P_{k},Q_{k}\in\mathbb{C}^{N\times\bar{k}}$,
$L_{k}\in\mathbb{C}^{\bar{k}\times\bar{k}}$, to $\mathbb{L}$ is thus computed.
If $U_{k}S_{k}V_{k}^{*}=L_{k}$ denotes the SVD of $L_{k}$, then the matrices
$P_{k}U_{k}$ and $V_{k}^{*}Q_{k}^{*}$ can be used in place of ${\bf X}_{n}$
and ${\bf Y}_{n}$ in (41) to compute the reduced model. The method is stopped
whenever the norm of the residual matrix
$\bm{M}P_{k}L_{k}Q_{k}^{*}-P_{k}L_{k}Q_{k}^{*}\bm{\Lambda}-{\bf V}{\bf R}+{\bf
L}{\bf W}$, consisting of the left-hand side of the Sylvester equation with
$\mathbb{L}$ replaced by its low-rank approximation $P_{k}L_{k}Q_{k}^{*}$, is
smaller than a certain threshold $\varepsilon$. In the results that follow we
employ $\varepsilon=10^{-4}$, as done in FSVDL . At each iteration step, the
SVD of $L_{k}$ is truncated to keep only the $n+p$ significant values.

We consider the Half&Half partition of the frequencies as this is the best
scenario for the scheme coming from FSVDL . The Half&Half partition often
leads to a rather fast convergence of the Galerkin-ADI method in terms of
number of iterations so that a quite small approximation space is constructed.
If different partitions were used, the Galerkin-ADI method could be equipped
with a quite involved divide-and-conquer scheme; see FSVDL . On the other
hand, as illustrated in Example 1, the Half&Half partition leads to higher
values of the HSS-rank of $\widetilde{C}$ than for the Even&Odd partition with
a consequent increment in the computational efforts of our scheme. In
addition, as for FSVDL , our tests employed complex arithmetic and did not
solve the corresponding Sylvester equation (37) for real-coefficient matrices.

In Table 4 we report the results for $p=10$, $n=50$, and different values of
$N$. Notice that even though the Galerkin-ADI approach efficiently computes
the approximation spaces, the construction of the reduced model (41) still
requires the allocation of both $\mathbb{L}$ and $\mathbb{S}$. Therefore, also
for the Galerkin-ADI scheme severe memory constraints hold and for
$N>30\,000$, we are not able to allocate the $\mathbb{L}$ and $\mathbb{S}$
matrices with complex entries on the machine used for running the tests.

| Galerkin-ADI with $\varepsilon=10^{-4}$ |  svds w/ $\widetilde{\mathcal{C}}$
---|---|---
| # of | Scheme | Total |  |  | Total |
$N$ | Iter. | Time(s) | Time (s) |  $\mathcal{H}_{2}$-error | $\mathtt{hssrank}(\widetilde{\mathcal{C}})$ | Time (s) |  $\mathcal{H}_{2}$-error
5 000 | 5 | 2.93 | 6.61 | 1.55e-2 | 42 | 100.21 | 2.06e-9
10 000 | 5 | 5.4 | 23.21 | 2.76e-2 | 46 | 174.04 | 1.81e-9
15 000 | 5 | 10.80 | 138.2 | 4.06e-2 | 49 | 280.11 | 1.31e-9
20 000 | 5 | 14.67 | 342.34 | 9.45e-2 | 50 | 340.42 | 1.17e-9
25 000 | 6 | 23.38 | 668.56 | 1.18e-2 | 52 | 426.34 | 9.55e-10
30 000 | 6 | 30.67 | 1198.81 | 1.17e-1 | 52 | 540.20 | 9.06e-10
Table 4: Example 3. Number of iterations, computational time (in seconds)
solely of the Galerkin-ADI iteration scheme together with the total time
(including building the data, the full Loewner and shifted Loewner matrices
and the projection step) as well as the $\mathcal{H}_{2}$-error achieved by
the Galerkin-ADI approach. In comparison, we list the HSS-rank, the total time
(in seconds) as well as the $\mathcal{H}_{2}$-error of the novel scheme
presented in this paper for different values of $N$ (number of samples),
$p=10$, and $n=50$.

Even though the Galerkin-ADI approach is faster for $N<20\,000$, the computed
approximation spaces are quite poor. Indeed, the computed reduced models are
always $7$ orders of magnitude less accurate than the ones constructed by our
approach. The paper FSVDL  validates the Galerkin-ADI scheme on a system with
randomly generated poles for various orders $n$ and number of samples $N$ but
does not mention the accuracy of the resulting models. Moreover, in terms of
CPU time, our results are comparable to the ones in FSVDL  when considering
the computational time solely of the Galerkin-ADI iteration, disregarding the
steps involving building the full matrices and projecting these to obtain the
reduced model.

The remarkable difference in the accuracy attained by the two approaches make
any sort of computational comparison rather pointless. However, we would like
to point out that the computational time of the Galerkin-ADI approach grows
quadratically with $N$ due to the need to assemble and store the full Loewner
and shifted Loewner matrices, while an $N\log N$ dependency of the
computational cost of our novel approach can be evidenced once again from the
timings reported in Table 4.

Several ideas could be implemented to improve the accuracy of the models
obtained with the Galerkin-ADI approach. In order to have the fairest
comparisons with respect to our novel approach, each of these ideas will be
tested separately to explore all the possibilities to enhance the Galerkin-ADI
approach from FSVDL .

First, the tolerance $\varepsilon$ for solving the Lyapunov equation via
Galerkin-ADI can be chosen to a value comparable to the noise level for an
$SNR$ of $120$, namely $\varepsilon=10^{-12}$. Results are detailed in Table 5
only for the case $N=5\,000$, $p=10$, and $n=50$ as the trend is obvious from
this one example. While the accuracy of the model has slightly improved with
respect to results obtained for $\varepsilon=10^{-4}$, the number of
iterations has also considerably increased, leading to matrices $L_{k}$ of
much larger dimensions for which the SVD $L_{k}=U_{k}S_{k}V_{k}^{*}$ becomes
costly. Hence, the CPU cost of the scheme has exploded and is no longer
viable. In any case, even for a tolerance value close to the noise level, the
accuracy of the model is several orders of magnitude worse than with our
proposed technique ($10^{-3}$ versus $10^{-9}$).

$\varepsilon$ | # of Iter. | Scheme Time(s) | Total Time(s) |  $\mathcal{H}_{2}$-error
---|---|---|---|---
$10^{-4}$ | 5 | 2.93 | 6.61 | 1.55e-2
$10^{-12}$ | 51 | 1648.02 | 1655.01 | 2.20e-3
Table 5: Example 3. Number of iterations, computational time (in seconds)
solely of the Galerkin-ADI iteration scheme together with the total time
(including building the data, the full Loewner and shifted Loewner matrices
and the projection step) as well as the $\mathcal{H}_{2}$-error achieved by
the Galerkin-ADI approach for $N=5\,000$, $p=10$, and $n=50$.

Second, it is always advisable to compute the projection subspaces from a
linear combination of $\mathbb{S}$ and $\mathbb{L}$, namely
$\mathbb{S}-x\mathbb{L}$ rather than only $\mathbb{L}$, as the Loewner matrix
$\mathbb{L}$ encodes the strictly rational part and the addition of
$\mathbb{S}$ provides all the information on the system, including its
polynomial part (the ${\bf D}$-term). We apply the low-rank Galerkin-ADI
method to the Sylvester equation fulfilled by $\mathbb{S}-x\mathbb{L}$ thus
computing a matrix $P_{k}Z_{k}Q_{k}^{*}$ such that
$P_{k}Z_{k}Q_{k}^{*}\approx\mathbb{S}-x\mathbb{L}$. Results are detailed in
Table 6 for the case $\varepsilon=10^{-4}$, $N=5\,000$, $p=10$, and $n=50$.
For all instances considered, results were comparable in terms of CPU time to
those obtained when considering solely the Sylvester equation satisfied by
$\mathbb{L}$ in the Galerkin-ADI iteration (listed in the first line of Table
6 for reference), while in terms of accuracy, they are slightly worse. For
this example, the sole benefit of using a linear combination
$\mathbb{S}-x\mathbb{L}$ might be the system identification properties as, in
principle, a sharp drop in the SVD of $Z_{k}$ reveals the degree of the
underlying system.

| # of Iter. | Scheme Time(s) | Total Time(s) |  $\mathcal{H}_{2}$-error
---|---|---|---|---
$\mathbb{L}$ | 5 | 2.93 | 6.61 | 1.55e-2
$\mathbb{S}$ | 6 | 3.96 | 6.57 | 5.22e-2
$\mathbb{S}-x\mathbb{L}$, $x=f(1)$ | 4 | 2.97 | 9.51 | 9.23e-2
$\mathbb{S}-x\mathbb{L}$, $x=f(N/2)$ | 4 | 2.89 | 9.34 | 9.23e-2
$\mathbb{S}-x\mathbb{L}$, $x=f(N)$ | 4 | 2.95 | 9.55 | 9.23e-2
Table 6: Example 3. Number of iterations, computational time (in seconds)
solely of the Galerkin-ADI iteration scheme together with the total time
(including building the data, the full Loewner and shifted Loewner matrices
and the projection step) as well as the $\mathcal{H}_{2}$-error achieved by
the Galerkin-ADI approach on the Sylvester equations satisfied by
$\mathbb{L}$, $\mathbb{S}$ and $\mathbb{S}-x\mathbb{L}$, for
$\varepsilon=10^{-4}$, $N=5\,000$, $p=10$, and $n=50$.

The third avenue worth exploring is employing real arithmetic and the
corresponding Sylvester equations (37) and (38). Table 7 shows the results
obtained using real arithmetic, both for the Galerkin-ADI scheme, as well as
our proposed method. For reference, the first line in Table 7 lists the
results previously obtained in complex arithmetic. For the method in FSVDL ,
the cost of the scheme has mostly increased, due to more complicated Sylvester
equations in (37) and (38). The CPU cost of building the data matrices, the
full Loewner and shifted Loewner matrices has also increased, yielding a total
cost far superior to that obtained in complex arithmetic. In some instances,
the accuracy has improved slightly. On the other hand, the real arithmetic
causes the HSS-rank of the Cauchy matrix approximation to be much smaller with
a remarkable impact on the CPU time and almost no effects on the model
accuracy when using our novel approach.

| Galerkin-ADI with $\varepsilon=10^{-4}$ |  svds w/ $\widetilde{\mathcal{C}}$
---|---|---
| # of | Scheme | Total |  |  | Total |
$N$ | Iter. | Time(s) | Time (s) |  $\mathcal{H}_{2}$-error | $\mathtt{hssrank}(\widetilde{\mathcal{C}})$ | Time (s) |  $\mathcal{H}_{2}$-error
$\mathbb{L}$ complex | 5 | 2.93 | 6.61 | 1.55e-2 | 42 | 100.21 | 2.06e-9
$\mathbb{L}$ | 3 | 1.76 | 50.43 | 1.03e-2 | 24 | 81.56 | 2.05e-9
$\mathbb{S}$ | 6 | 3.61 | 52.42 | 2.66e-3 | 24 | 80.97 | 2.05e-9
$\mathbb{S}-x\mathbb{L}$, $x=f(1)$ | 7 | 18.57 | 65.84 | 2.53e-3 | 24 | 80.84 | 2.05e-9
$\mathbb{S}-x\mathbb{L}$, $x=f(N/2)$ | 5 | 15.13 | 63.63 | 1.87e-2 | 24 | 82.11 | 2.05e-9
$\mathbb{S}-x\mathbb{L}$, $x=f(N)$ | 6 | 16.47 | 66.72 | 9.01e-2 | 24 | 82.95 | 2.05e-9
Table 7: Example 3. Number of iterations, computational time (in seconds)
solely of the Galerkin-ADI iteration scheme together with the total time
(including building the data, the full Loewner and shifted Loewner matrices
and the projection step) as well as the $\mathcal{H}_{2}$-error achieved by
the Galerkin-ADI approach. In comparison, we list the HSS-rank, the total time
(in seconds) as well as the $\mathcal{H}_{2}$-error of the novel scheme
presented in this paper for different values of $N$ (number of samples),
$p=10$, and $n=50$ when employing real arithmetic.

We conclude this example by mentioning that the use of a _hybrid_ approach may
be fruitful. In particular, our novel approach can be employed to avoid
storing the large and dense Loewner and shifted Loewner matrices. Then, the
Galerkin-ADI scheme can be used to compute the first dominant singular vectors
of $\mathbb{S}-x\mathbb{L}$, instead of employing svds, thus also being able
to identify the order of the underlying system. However, the accuracy will not
be comparable to that of our proposed approach. We implemented this idea and
list the CPU times of the various steps in Table 8 together with the resulting
accuracy for Galerkin-ADI applied to solving the Sylvester equation (37) for
$\mathbb{L}$ in real arithmetic with $\varepsilon=10^{-4}$ for $N=5\,000$,
$p=10$, and $n=50$. Plots of the responses of our proposed approach, together
with the Galerkin-ADI scheme as proposed in FSVDL  and the hybrid approach are
shown in Figure 5. Even though the general shape of the response is well
captured, some resonances are not modeled accurately, as expected from the
much higher model errors reported earlier. This can be noticed better from the
error plots in Figure 6.

Data matrices Time(s) | Galerkin-ADI Scheme Time(s) | Projection Time(s) | Total Time(s) |  $\mathcal{H}_{2}$-error
---|---|---|---|---
0.8 | 1.76 | 4.65 | 7.21 | 2.1e-2
Table 8: Example 3. Computational time (in seconds) of the three individual
steps in the _hybrid_ approach: setting up of the data matrices, the Galerkin-
ADI iteration scheme and projection to obtain the reduced model, together with
the total time as well as the $\mathcal{H}_{2}$-error for $N=5000$, $p=10$,
and $n=50$ in real arithmetic.

Figure 5: Example 3. Frequency response of the model (in black) and the
measurements (in red) for $N=5\,000$, $p=10$, and $n=50$ using our proposed
approach, Galerkin-ADI as in FSVDL  and the hybrid approach, employing real
arithmetic.

Figure 6: Example 3. Error plots for $N=5\,000$, $p=10$, and $n=50$ using our
proposed approach, Galerkin-ADI as in FSVDL  and the hybrid approach,
employing real arithmetic.

##  5 Conclusion

By exploiting the Cauchy-like structure of the Loewner and shifted Loewner
matrices, a novel strategy for reducing the computational costs and the memory
requirements of the Loewner framework has been proposed. In particular, the
use of the HSS-format leads to tremendous savings in the storage demand and
computational efforts of the overall scheme. Indeed, except for the
construction of $\widetilde{\mathcal{C}}$ whose cost is polylogarithmic in
$N$, both the memory requirements and the computational cost of iteratively
performing the SVD now linearly depend on the cardinality of the considered
data set.

The success of our procedure strongly relies on the capability of representing
the Cauchy matrix $\mathcal{C}$ in terms of an HSS-matrix
$\widetilde{\mathcal{C}}$ with low $\left(\alpha,\beta\right)$ rank of the
off-diagonal blocks. Even though we restricted ourselves to showing how
different, but common, partitions of the frequencies affect the HSS-rank of
$\widetilde{\mathcal{C}}$, a thorough analysis of their connection may be
beneficial. Moreover, we have always computed $\widetilde{\mathcal{C}}$ at
high accuracy. We believe that the employment of more inexact, and thus with a
lower rank, HSS-representations of $\mathcal{C}$ and its effects on the
accuracy of the overall scheme may be another interesting research direction
which is worth pursuing depending on the application at hand.

The strategy presented in this paper can be applied to more sophisticated
problems as long as the Loewner and shifted Loewner matrices maintain a
Cauchy-like structure. In particular, our approach can be employed with minor
modifications in model order reduction of parametrized IonitaPhD , linear
switched GosPA18 , and bilinear systems AntGI16 .

## Acknowledgements

We are in debt with Leonardo Robol for some help with hmtoolbox  and fruitful
discussions about the topic of this paper. His assistance is greatly
appreciated. We also thank Peter Benner and Jens Saak for insightful comments
on earlier versions of the manuscript.

The first author is member of the Italian INdAM Research group GNCS.

## Declarations

The research presented in this paper is based upon work supported by the
National Science Foundation under Grant No. DMS-1439786 while both the authors
were in residence at the Institute for Computational and Experimental Research
in Mathematics (ICERM) in Providence, RI, during the _Model and Dimension
Reduction in Uncertain and Dynamic Systems_ program. Even though the second
half of the program had to be performed virtually due to the restrictions
caused by the COVID-19 pandemic, we are extremely grateful to the organizers
of the program and the whole staff of ICERM for doing whatever possible to
maintain an exciting, fruitful, and high-quality working environment.

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

The datasets and algorithms generated during and/or analysed during the
current study are available from the corresponding author on reasonable
request. Moreover, the approach presented in this paper will be included in
the hm-toolbox in the near future.

## References

  * (1) A. C. Antoulas, I. V. Gosea, and A. C. Ioniţă, Model reduction of bilinear systems in the Loewner framework, SIAM Journal on Scientific Computing, 38 (2016), pp. B889–B916.
  * (2) A. C. Antoulas, S. Lefteriu, and A. C. Ioniţă, A tutorial introduction to the Loewner framework for model reduction, ch. 8, pp. 335–376.
  * (3) J. Baglama and L. Reichel, Augmented implicitly restarted Lanczos bidiagonalization methods, SIAM J. Sci. Comput., 27 (2005), pp. 19–42.
  * (4) M. Bebendorf, Approximation of boundary element matrices, Numer. Math., 86 (2000), pp. 565–589.
  * (5) B. Beckermann and A. Townsend, Bounds on the singular values of matrices with displacement structure, SIAM Rev., 61 (2019), pp. 319–344.
  * (6) A. Bouras and V. Frayssé, Inexact matrix-vector products in Krylov methods for solving linear systems: a relaxation strategy, SIAM J. Matrix Anal. Appl., 26 (2005), pp. 660–678.
  * (7) J. Carrier, L. Greengard, and V. Rokhlin, A fast adaptive multipole algorithm for particle simulations, SIAM J. Sci. Statist. Comput., 9 (1988), pp. 669–686.
  * (8) S. Chandrasekaran, P. Dewilde, M. Gu, W. Lyons, and T. Pals, A fast solver for HSS representations via sparse matrices, SIAM J. Matrix Anal. Appl., 29 (2006/07), pp. 67–81.
  * (9) S. Chandrasekaran, M. Gu, X. Sun, J. Xia, and J. Zhu, A superfast algorithm for Toeplitz systems of linear equations, SIAM J. Matrix Anal. Appl., 29 (2007), pp. 1247–1266.
  * (10) A. S. Derakhtenjani, J. A. Candanedo, Y. Chen, V. R. Dehkordi, and A. K. Athienitis, Modeling approaches for the characterization of building thermal dynamics and model-based control: A case study, Science and Technology for the Built Environment, 21 (2015), pp. 824–836.
  * (11) P. Drineas, M. W. Mahoney, and S. Muthukrishnan, Relative-Error $CUR$ Matrix Decompositions, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 844–881.
  * (12) M. Embree and A. C. Ioniţă, Pseudospectra of Loewner matrix pencils, (2019).  ArXiv preprint: 1910.12153.
  * (13) M. A. Freitag and A. Spence, Convergence theory for inexact inverse iteration applied to the generalised nonsymmetric eigenproblem, Electron. Trans. Numer. Anal., 28 (2007/08), pp. 40–64.
  * (14) S. W. Gaaf and V. Simoncini, Approximating the leading singular triplets of a large matrix function, Applied Numerical Mathematics, 113 (2017), pp. 26 – 43.
  * (15) I. Gohberg and V. Olshevsky, Fast algorithms with preprocessing for matrix-vector multiplication problems, J. Complexity, 10 (1994), pp. 411–427.
  * (16) I. V. Gosea, M. Petreczky, and A. C. Antoulas, Data-driven model order reduction of linear switched systems in the Loewner framework, SIAM Journal on Scientific Computing, 40 (2018), pp. B572–B610.
  * (17) L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, Journal of Computational Physics, 73 (1987), pp. 325 – 348.
  * (18) A. Hochman, Fast singular-value decomposition of Loewner matrices for state-space macromodeling, in 2015 IEEE 24th Electrical Performance of Electronic Packaging and Systems (EPEPS), 2015, pp. 177–180.
  * (19) M. E. Hochstenbach, A Jacobi–Davidson type SVD method, SIAM J. Sci. Comput., 23 (2001), pp. 606–628.
  * (20) R. Horn and C. Johnson, Topics in Matrix Analysis, Cambridge Univ. Press, Cambridge, UK, 1991.
  * (21) A. C. Ioniţă, Lagrange rational interpolation and its applications to approximation of large-scale dynamical systems, PhD thesis, Rice University, Aug. 2013.
  * (22) D. Karachalios, I. Gosea, and A. Antoulas, Data-driven approximation methods applied to non-rational functions, Proc. Appl. Math. Mech., 18 (2018).
  * (23) D. Kressner, S. Massei, and L. Robol, Low-rank updates and a divide-and-conquer method for linear matrix equations, SIAM J. Sci. Comput., 41 (2019), pp. A848–A876.
  * (24) P. Kürschner and M. Freitag, Inexact methods for the low rank solution to large scale Lyapunov equations, BIT Numerical Mathematics, (2020).
  * (25) R. Larsen, Lanczos bidiagonalization with partial reorthogonalization, DAIMI Report Series, 27 (1998).
  * (26) S. Lefteriu and A. C. Antoulas, A New Approach to Modeling Multiport Systems From Frequency-Domain Data, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 29 (2010), pp. 14–27.
  * (27) S. Massei, D. Palitta, and L. Robol, Solving rank-structured Sylvester and Lyapunov equations, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 1564–1590.
  * (28) S. Massei, L. Robol, and D. Kressner, hm-toolbox: MATLAB software for HODLR and HSS matrices, SIAM J. Sci. Comput., 42 (2020), pp. C43–C68.
  * (29) MATLAB, version 9.9.0.1467703 (R2020b), The MathWorks Inc., Natick, Massachusetts, 2020.
  * (30) A. J. Mayo and A. C. Antoulas, A framework for the solution of the generalized realization problem, Linear Algebra and Its Applications, 405 (2007), pp. 634–662.
  * (31) Y. Nakatsukasa, Fast and stable randomized low-rank matrix approximation, 2020.  ArXiv preprint: 2009.11392.
  * (32) D. Palitta and P. Kürschner, On the convergence of low-rank Krylov methods, (2021).  Accepted for publication in Numerical Algorithms. ArXiv preprint: 1909.01226.
  * (33) V. Y. Pan, Fast approximate computations with Cauchy matrices, polynomials and rational functions, in Computer Science - Theory and Applications, E. A. Hirsch, S. O. Kuznetsov, J.-É. Pin, and N. K. Vereshchagin, eds., Cham, 2014, Springer International Publishing, pp. 287–299.
  * (34) V. Y. Pan, Transformations of matrix structures work again, Linear Algebra and its Applications, 465 (2015), pp. 107 – 138.
  * (35) B. Peeters, H. Van der Auweraer, P. Guillaume, and J. Leuridan, The PolyMAX frequency-domain method: A new standard for modal parameter estimation?, Shock and Vibration, 11 (2004), pp. 395–409.
  * (36) C. Poussot-Vassal, D. Quero, and P. Vuillemin, Data-driven approximation of a high fidelity gust-oriented flexible aircraft dynamical model, IFAC-PapersOnLine, 51 (2018), pp. 559 – 564.  9th Vienna International Conference on Mathematical Modelling.
  * (37) M. Sahouli and A. Dounavis, Iterative Loewner matrix macromodeling using CUR decomposition for noisy frequency responses, in 2019 IEEE 28th Conference on Electrical Performance of Electronic Packaging and Systems (EPEPS), 2019, pp. 1–3.
  * (38) V. Simoncini and L. Eldén, Inexact Rayleigh quotient-type methods for eigenvalue computations, BIT, 42 (2002), pp. 159–182.
  * (39) V. Simoncini and D. B. Szyld, Theory of inexact Krylov subspace methods and applications to scientific computing, SIAM J. Sci. Comput., 25 (2003), pp. 454–477.
  * (40) M. Stoll, A Krylov-Schur approach to the truncated SVD, Linear Algebra Appl., 436 (2012), pp. 2795–2806.
  * (41) The MORwiki Community, MORwiki - Model Order Reduction Wiki.  http://modelreduction.org.
  * (42) J. van den Eshof and G. L. G. Sleijpen, Inexact Krylov subspace methods for linear systems, SIAM J. Matrix Anal. Appl., 26 (2004), pp. 125–153.
  * (43) R. Vandebril, M. Van Barel, G. Golub, and N. Mastronardi, A bibliography on semiseparable matrices, Calcolo, 42 (2005), pp. 249–270.
  * (44) J. Vogel, J. Xia, S. Cauley, and V. Balakrishnan, Superfast divide-and-conquer method and perturbation analysis for structured eigenvalue solutions, SIAM Journal on Scientific Computing, 38 (2016), pp. A1358–A1382.
  * (45) Y. Xi, J. Xia, S. Cauley, and V. Balakrishnan, Superfast and stable structured solvers for Toeplitz least squares via randomized sampling, SIAM J. Matrix Anal. Appl., 35 (2014), pp. 44–72.
  * (46) J. Xia, S. Chandrasekaran, M. Gu, and X. S. Li, Fast algorithms for hierarchically semiseparable matrices, Numerical Linear Algebra with Applications, 17 (2010), pp. 953–976.
  * (47) J. Xia, S. Chandrasekaran, M. Gu, and X. S. Li, Superfast multifrontal method for large structured linear systems of equations, SIAM Journal on Matrix Analysis and Applications, 31 (2010), pp. 1382–1411.