  Large-scale linear, time-invariant (LTI) dynamical systems are widely used to
characterize complicated physical phenomena. We propose a two-stage algorithm
to reduce the order of a large-scale LTI system given samples of its transfer
function for a target degree $k$ of the reduced system. In the first stage, a
modified adaptive Antoulas--Anderson (AAA) algorithm is used to construct a
degree $d$ rational approximation of the transfer function that corresponds to
an intermediate system, which can be numerically stably reduced in the second
stage using ideas from the theory on Hankel norm approximation (HNA). We also
study the numerical issues of Glover's HNA algorithm and provide a remedy for
its numerical instabilities. A carefully computed rational approximation of
degree $d$ gives us a numerically stable algorithm for reducing an LTI system,
which is more efficient than SVD-based algorithms and more accurate than
moment-matching algorithms.
