# Learning Generalized Causal Structure in Time-series

Aditi Kathpalia
Department of Complex Systems  Institute of Computer Science of the Czech
Academy of Sciences
Prague  Czech Republic
<EMAIL_ADDRESS>
Keerti Panchakshari Charantimath
Department of Mathematics  Indian Institute of Technology  Kharagpur
West Bengal  India
<EMAIL_ADDRESS>
Nithin Nagaraj
Consciousness Studies Programme  National Institute of Advanced Studies
Indian Institute of Science Campus  Bengaluru  India
<EMAIL_ADDRESS>
###### Abstract

The science of causality explains/determines ‘cause-effect’ relationship
between the entities of a system by providing mathematical tools for the
purpose. In spite of all the success and widespread applications of machine-
learning (ML) algorithms, these algorithms are based on statistical learning
alone. Currently, they are nowhere close to ’human-like’ intelligence as they
fail to answer and learn based on the important “Why?” questions. Hence,
researchers are attempting to integrate ML with the science of causality.
Among the many causal learning issues encountered by ML, one is that these
algorithms are dumb to the temporal order or structure in data. In this work
we develop a machine learning pipeline based on a recently proposed
‘neurochaos’ feature learning technique (_ChaosFEX_ feature extractor), that
helps us to learn generalized causal-structure in given time-series data.

##  1 Introduction

Current machine learning (ML) algorithms are based on statistical learning,
that is, these algorithms learn mostly by identifying associations in the
given data. As the popularity and range of ML algorithms have increased to
deal with specific tasks of classification and prediction, researchers focus
on the task of bringing these ML algorithms closer to human-like intelligence
so that these algorithms can be used to deal with _higher-level_ problems. One
of the major ways in which human intelligence differs from machine
intelligence is that humans learn through ‘cause-effect’ reasoning. Put more
specifically, humans have the capability to answer _‘what if’_ kind of
questions. _What if I intervene and do something to a given system? What if I
had acted in a different way?_ Machine intelligence is still far away from
answering these kind of questions [1, 2].

Let alone learning based on a causal understanding, machines are not even
capable of making sense of temporal or sequential order in the data. Further,
they lack the ability of _generalization_ which involves transfer of learning
from one problem to another and learning features that categorize together
datasets that are more alike than different. This lack of generalization stems
mainly from an inability of these algorithms to learn causal structures [2].

While a lot of research is being currently pursued at the intersection of
causality and machine learning, we focus on causal learning based on time
series data. Time series data are real valued data points arranged in a
chronological order. Such data is often collected as measurements in different
fields such as finance, medicine, climate science etc. Currently, there are
two main tasks for which causal inference is done using time series data. One
is studying the effect of treatment or interventions on certain variable of a
system. This treatment may be provided as a discrete event or may be
continuously given over time. This task is generally referred to as _causal
treatment effect estimation_ and finds applications in estimating the effects
of policies implemented so as to increase/ reduce the sale of certain goods or
in estimating the effect of drugs given to patients [3, 4, 5]. The second one
is discovering causal relationships between different variables of a system,
for which temporal evolution data is available. This task is generally
referred to as _causal network discovery_ and is useful in domains where we
wish to study the interaction between different variables such as the role of
human and natural factors in climate change [6, 7, 5].

For our work, we take one step backward and ask the question of whether ML
algorithms can even identify whether given time series has an inherent causal
structure associated with it, that is, whether the past values of a time
series effect its present values or if that is not the case. Learning this
structure for time series is essential for making sense of any ordered data
given as input to ML algorithms and is also essential to meaningfully answer
the two causality tasks discussed above. There are of course mathematical
techniques that rely on fitting autoregressive models and testing for time
dependence by estimating serial correlation in given time series. The latter
techniques include the Durbin–Watson test [8], Ljung–Box test [9] and
Breusch–Godfrey test [10]. However, these methods have their limitations: they
are parametric, requiring the choice of maximum order; do not work with
overlapping data or in the case of heteroskedasticity in the error process.
Further, since they are purely statistical, based on autocorrelation, they
cannot handle data with intervention(s)/ perturbation(s) where the cause-
effect relation in time series builds up as a result of external event(s). We
are interested in developing a non-parametric way of identifying causal
structure in time series data that works irrespective of the underlying model
and learns _generalized_ causal-structure in data, unaffected by distribution
or model shifts. For this task, we rely on extracting strong features which
can identify whether or not there is an underlying causal-structure and later
using a simple ML algorithm, classify given time-series data based on having a
causal-structure or not. We have used time series values directly, frequency
domain characteristics of the process, and ‘neurochaos’ based features
extracted from the frequency representation of the process to train the ML
classifiers in separate models. _Neurochaos_ inspired feature learning [11]
has been proposed based on a _ChaosNet_ neural network architecture [12] that
is composed of neurons, which individually mimic chaotic firing properties of
the human brain neurons. We compare the performance of these models and show
that neurochaos based learning is able to extract meaningful features for
learning generalized causal-structure in time series data.

This paper is organized as follows: Section 2 describes the time series data
simulated for the study. Section 3 describes the methods, that is the features
and the classifiers used for distinguishing between causal and non-causal
structure in temporal data. Results are given in Section 4 and we conclude
with Section 5, discussing the results, and hinting at future research
directions.

##  2 Datasets

The following time-series datasets were simulated and used in this work: Time-
series having a causal-structure or causal time-series were generated using
the following models:

  * •

Autoregressive (AR) processes

AR processes are random processes in which value of the time-series at any
time depends linearly on its past values and on a noise term. The general
equation governing an AR ($p$) process, $X$, is given as:

$X(t)=c+\sum_{i=1}^{p}{a_{i}X(t-i)}+\varepsilon_{t},$ (1)

where, $t$ denotes time, $c$ is a constant, $p$ is the order of the AR
process, $a_{i}$ is the coefficient at lag $i$ time step(s) and
$\varepsilon_{t}$ is the noise term at time $t$. Order of an AR process is the
maximum past lag on which the value of the process at any time depends.

AR processes are used to model several processes occurring in nature,
economics etc. [13, 14]

  * •

Autoregressive Moving Average (ARMA) processes

These processes have an AR part and an MA part. Based on the AR part, each
term in the time-series is regressed on its past values and based on the MA
part, the noise term at any time point is modelled as a linear combination of
instantaneous noise term (at that time point) and noise terms at past time
points. In mathematical form, an ARMA($p,q$) process $X$ can be expressed as:

$X(t)=c+\sum_{i=1}^{p}{a_{i}X(t-i)}+\sum_{i=0}^{q}{b_{i}\varepsilon(t-i)},$
(2)

where, $t$ denotes time, $c$ is a constant, $p$ is the order of the AR part
and $q$ is the order of the MA part, $a_{i}$ is the AR coefficient at lag $i$
time step(s) and $b_{i}$ is the MA coefficient at lag $i$ time step(s),
$\varepsilon(t)$ is the noise term at time $t$.

ARMA processes are widely used in the modelling of financial and geological
processes [15, 13].

  * •

Autoregressive Fractionally Integrated Moving Average (ARFIMA) These processes
are used to model long term memory processes. An ARMA($p^{\prime},q$) process
when expressed in terms of the lag or backshift operator can be written as:

$(1-\sum_{i=1}^{p^{\prime}}{a_{i}B^{i}})X(t)=(1+\sum_{i=1}^{q}{b_{i}}B^{i})\varepsilon(t),$
(3)

Now, if the polynomial $(1-\sum_{i=1}^{p^{\prime}}{a_{i}B^{i}})$ has a unit
root of multiplicity $d$, then:

$(1-\sum_{i=1}^{p^{\prime}}{a_{i}B^{i}})=(1-\sum_{i=1}^{p^{\prime}-d}{a_{i}B^{i}})(1-B)^{d},$
(4)

This is the property expressed by Autoregressive Integrated Moving Average or
ARIMA($p,d,q$) processes with $p=p^{\prime}-d$. For the case of ARFIMA
processes, the difference parameter $d$ is allowed to take non-integer values.
They are thus generally expressed as:

$(1-\sum_{i=1}^{p^{\prime}}{a_{i}B^{i}})(1-B)^{d}X(t)=(1+\sum_{i=1}^{q}{b_{i}}B^{i})\varepsilon(t),$
(5)

ARFIMA models are also widely used, for example in the modelling of economic,
geological and physiological time series [16, 17].

Time-series without a causal structure or non-causal time-series were
generated as having time-series value at each time point, as independently and
randomly chosen from a normal distribution, $\mathcal{N}(\mu,\sigma^{2})$
(where $\mu$ and $\sigma$ denote the mean and standard deviation of the
distribution) or uniform distribution, $U(b_{u},b_{l})$ (where $b_{u}$ and
$b_{l}$ denote the upper bound and lower bound of the distribution).

##  3 Methods

###  3.1 Using time-series values directly

In this case, time-series values were directly passed through ML classifiers
to check if they could classify time series with a causal structure and
without one appropriately. For this purpose, Logistic Regression (LR) [18] and
Long-short term memory (LSTM) [19] classifiers were used. LSTM classifiers
have a recurrent neural network architecture and include not only feedforward
but also feedback connections. Along with being applied to single data points,
they have also been shown to work for large sequential datasets such as speech
and text data.

###  3.2 Using frequency domain characteristics of time-series

Any periodic structure in a time series is reflected in its frequency domain
amplitude. Since continuous causal effect within a time series (where its past
values affect present values), can be thought of as imbuing some periodic
nature to the time series, we consider it better to analyze the signals in
frequency domain. The fast fourier transform algorithm was used to obtain the
discrete fourier transform of simulated signals. Frequency domain amplitudes
of the signal were then directly passed through ML classifiers: LR and LSTM.

###  3.3 Using ChaosFEX features of the frequency domain signal

_ChaosNet_ is an artificial neural network inspired by the chaotic firing of
neurons in the human brain [12]. Based on ChaosNet, a hybrid learning
architecture was proposed later that uses ‘neurochaos’ or ‘ChaosFEX’ 111The
code for ChaosFEX feature extraction is available at
https://github.com/pranaysy/ChaosFEX and this code was used in our work. This
code is available open source under the license: Apache License, Version 2.0.
Consent of the authors was taken to use the code. features combined with
traditional machine learning approaches for the purpose of classification
[11]. The basis of ChaosNet or the extraction of ChaosFEX features is the
_topological transitivity_ property of Chaos. The ChaosNet architecture is
composed of 1-D chaotic neurons called as Generalized Luröth Series (GLS)
maps. The number of neurons in the architecture depend upon the number of
inputs/features provided to the architecture, with each neuron receiving a
single value as input for each data instance. These neurons are set to have an
initial activity of $q$ units. Input data is first normalized (to lie between
$0$ and $1$) and then passed to these GLS neurons. Each GLS neuron starts
firing and keeps on firing until it reaches the epsilon neighborhood of the
input value (also called stimulus) provided to it. The epsilon neighborhood of
a stimulus $y$ is defined as the interval $(y-\varepsilon,y+\varepsilon)$,
where $\varepsilon>0$.

In this work, we use a GLS neuron $T:[0,1)\rightarrow[0,1)$ defined by the
following equation:

$T(y)=\begin{cases}\frac{y}{b},&0\leq y<b,\\\ \frac{1-y}{1-b},&b\leq
y<1,\end{cases}$ (6)

where, $y\in[0,1)$ and $0<b<1$. Let the trajectory of a GLS neuron for a
stimulus $y$ be given by $A=[q\rightarrow T(q)\rightarrow
T^{2}(q)\ldots\rightarrow T^{N}(q)]$. Thus, it can be seen that the neuron
takes $N$ time steps to reach the epsilon neighborhood of the stimulus. $N$ is
referred to as the firing time. The fraction of the time when the trajectory
is above the discrimination threshold ($b$) is referred to as the _Topological
Transitivity-Symbolic Sequence (TTSS) Feature_ for the stimulus $y$.

To classify time-series as having a causal structure or not, normalized
frequency domain amplitudes of given time series were passed as input to the
ChaosFEX feature extractor and the TTSS features extracted from the same
(there will be one _TTSS_ feature at each frequency) were then passed as
inputs to the ML classifier, LR. This scheme, thus exploited the neurochaos
based hybrid learning architecture for our task. Such an architecture has been
shown to learn from finite sample datasets and can help to use chaos as a
kernel trick which can be useful to make given data linearly separable.

##  4 Results

The results obtained by using each of the methods discussed above are detailed
in the subsections below. We term each of our methods as models, as what the
methods are essentially trying to do is distinguish between a causal and a
non-causal underlying structure. So, even though we are not strictly trying to
fit a model, we are trying to learn a generalized characteristic of the
temporal order in given data.

Each of the models were trained and tested using the AR training and testing
set which consisted of AR series as the causal time-series and time series
with independent entries, randomly chosen from normal distribution,
$\mathcal{N}(0,0.01)$, as the non-causal time series. We term time-series
generated in the latter way as random time-series. 1250 AR and 1250 random
time-series, each having a length of 2000 time points were simulated. Each AR
series was of the form, $X(t)=a_{k}X(t-k)+\varepsilon_{t}$, with the order
$k$, being randomly chosen between 1 and 20. These series were initialized to
random values. The noise term followed the distribution $\mathcal{N}(0,0.01)$
and the AR coefficient for each simulation was randomly chosen such that
$a_{k}\in U(0.8,0.9)$. The training to testing split for this dataset was
70:30.

For further testing of the models, the following datasets were generated.
These testing sets had a shift in the probability distribution of time-series
when compared to the training set. We call them Distribution shift testing set
I and Distribution shift testing set II. For both of these, causal time-series
were generated in exactly the same manner as for the AR training and testing
set. The non-causal time series in Distribution shift testing set I were
generated using $\mathcal{N}(0,0.09)$ and in Distribution shift testing set II
were generated using $U(-0.6,0.6)$. The number of non-causal time series in
both datasets were 1250 and each time-series was simulated with 2000 time
points. For these datasets, models were trained using the AR testing set
described in the previous paragraph and tested based on these datasets.

###  4.1 Time-series values model

Performance metrics for the time-series values model using both LR and LSTM on
the simulated training and testing sets are shown in Table 1. The non-causal
time series were labelled as Class-0 and the causal time series as Class-1.
The Precision, Recall and F1-Score columns are given in the format
$(\cdot,\cdot)$, with the first value being for Class-0 and the second value
for Class-1.

LR was implemented using the _Scikit-Learn_ [20] package and LSTM was
implemented using the _Keras_ [21] framework in Python. For LR, all parameters
were kept as default. For LSTM, the layers were stacked in the following
order: LSTM layer with 10 outputs, dropout layer with a parameter value of
0.5, dense hidden layer with 10 outputs and relu activation function and
finally, a dense output layer with two outputs and soft-max activation
function. This classification model used categorical cross-entropy as the loss
function and adam optimizer for optimization. The optimizer was run for 50
epochs with a batch size of 1250. This classification model was run 15 times
to avoid local minima and the run which gave the best testing accuracy was
used.

Table 1: Prediction for time-series values model. LR: Logistic Regression, LSTM: Long Short-Term Memory. Scores are given in the format $(\cdot,\cdot)$, with the first value being for Class-0 and the second value for Class-1. Dataset |  Classifier |  Precision |  Recall |  $\mathbf{F_{1}Score}$ |  Accuracy
---|---|---|---|---|---
AR training set |  LR |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
|  LSTM |  (0.98, 1.00) |  (1.00, 0.98) |  (0.99, 0.99) |  99%
AR testing set |  LR |  (0.53, 0.54) |  (0.70, 0.37) |  (0.60, 0.44) |  54%
|  LSTM |  (0.98, 0.99) |  (0.99, 0.98) |  (0.99, 0.99) |  99%
Distribution shift |  LR |  (0.49, 0.48) |  (0.58, 0.38) |  (0.53, 0.43) |  48%
testing set I |  LSTM |  (0.15, 0.49) |  (0.01, 0.97) |  (0.01, 0.65) |  49%
Distribution shift |  LR |  (0.48, 0.47) |  (0.56, 0.38) |  (0.52, 0.42) |  47%
testing set II |  LSTM |  (0.00, 0.49) |  (0.00, 0.97) |  (0.00, 0.65) |  48%

It can be seen from Table 1 that the LR classifier works well for the AR
training set but fails for the AR testing set, indicating that it is
overfitting the training set. Tuning the hyperparameters for the classifier
did not help to improve its performance. LSTM gives good performance for both
AR training as well as testing set but fails when there is a distribution
shift in the random time-series. This is probably because the LSTM is
recognizing only the random time-series used for training as non-causal and
anything other than that as causal. Thus, it fails to learn the causal
structure in time-series.

###  4.2 Frequency-domain representation model

Performance metrics for the frequency-domain representation model using both
LR and LSTM are shown in Table 2 The Precision, Recall and F1-Score columns
are given in the format $(\cdot,\cdot)$, with the first value being for
Class-0 (non-causal time series) and the second value for Class-1 (causal time
series). In order to observe the characteristics of frequency amplitude of the
signals, it would be essential to demean the time-series so as to remove any
DC component that is present. However, since all our simulated datasets have a
zero mean, we skipped this step.

The hyperparameters and architecture used for LR and LSTM as well as the
packages used for the implementation remained the same as in Section 4.1.

Table 2: Prediction for frequency domain representation model. LR: Logistic Regression, LSTM: Long Short-Term Memory. Scores are given in the format $(\cdot,\cdot)$, with the first value being for Class-0 and the second value for Class-1. Dataset |  Classifier |  Precision |  Recall |  $\mathbf{F_{1}Score}$ |  Accuracy
---|---|---|---|---|---
AR training set |  LR |  (0.99, 1.00) |  (1.00, 0.99) |  (1.00, 1.00) |  100%
|  LSTM |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
AR testing set |  LR |  (0.99, 1.00) |  (1.00, 0.99) |  (0.99, 0.99) |  99%
|  LSTM |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
Distribution shift |  LR |  (0.00, 0.50) |  (0.00, 0.99) |  (0.00, 0.66) |  50%
testing set I |  LSTM |  (0.00, 0.50) |  (0.00, 1.00) |  (0.00, 0.67) |  50%
Distribution shift |  LR |  (0.00, 0.50) |  (0.00, 0.99) |  (0.00, 0.66) |  50%
testing set II |  LSTM |  (0.00, 0.50) |  (0.00, 1.00) |  (0.00, 0.67) |  50%

Figure 1 shows the amplitude spectrum of the fourier transformed signal for a
realization of AR(15) process. Figures 2(a) and 2(b) show the same for
realizations of random time-series generated using distributions
$\mathcal{N}(0,0.01)$ and $U(-0.6,0.6)$ respectively. It can be seen from
these figures that the frequency domain characteristics of causal and non-
causal time series used are quite different, yet the implemented ML
classifiers are unable to distinguish between the two.

Figure 1: Frequency amplitude spectrum for a realization of AR(15) process

(a) (b)

Figure 2: Frequency amplitude spectrum for a realization of random time-series
generated using $\mathcal{N}(0,0.01)$ (left) and $U(-0.6,0.6)$ (right).

It can be seen from Table 2 that both LR and LSTM accurately classify the AR
training and testing sets, however, fail for testing sets with distribution
shifts in non-causal time-series. This again seems to be because the
classifiers are recognizing only the random time-series used for training as
non-causal and anything other than that as causal and have failed to learn any
causal-structure from the data.

###  4.3 ChaosFEX feature representation model (FT+ChaosFEX)

Performance metrics when TTSS features of the amplitudes at different
frequencies were passed as inputs to the ML classifier logistic regression are
shown in Table 3. The Precision, Recall and F1-Score columns are given in the
format $(\cdot,\cdot)$, with the first value being for Class-0 (non-causal
time series) and the second value for Class-1 (causal time series).

Table 3: Prediction for FT+ChaosFEX model with Logistic Regression classifier. Scores are given in the format $(\cdot,\cdot)$, with the first value being for Class-0 and the second value for Class-1. Dataset |  Precision |  Recall |  $\mathbf{F_{1}Score}$ |  Accuracy
---|---|---|---|---
AR training set |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
AR testing set |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
Distribution shift testing set I |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
Distribution shift testing set II |  (1.00, 1.00) |  (1.00, 1.00) |  (1.00, 1.00) |  100%
AR(100) testing set |  (NA, 1.00) |  (NA, 0.99) |  (NA, 1.00) |  99%
ARMA testing set |  (NA, 1.00) |  (NA, 1.00) |  (NA, 1.00) |  100%
ARFIMA testing set |  (NA, 1.00) |  (NA, 1.00) |  (NA, 1.00) |  100%

It can be seen that this model classifies accurately not just the AR training
and testing sets but also both the testing sets with distribution shifts in
the non-causal time series. To further check the robustness of the model, it
was tested on more testing sets in which there was a shift in the distribution
of causal time-series. Three testing sets generated in this manner included
AR(100), ARMA and ARFIMA as causal time-series and no non-causal time series.
Since these datasets did not include any non-causal time series, the first
value for precision, recall and $F_{1}$ score are marked as NA (Not-
Applicable) in Table 3.

The specific details of the three new testing sets mentioned above are as
follows. Each of these datasets consisted of 1250 causal time-series with 2000
time points each. Each time series was initialized to random values and had an
instantaneous noise term, $\varepsilon_{t}\in\mathcal{N}(0,0.01)$. AR(100)
processes followed the form $X(t)=a_{100}X(t-100)+\varepsilon_{t}$, where, for
each realization, $a_{100}\in U(0.8,0.9)$. For each realization in the set of
ARMA and ARFIMA processes, the AR and MA orders were randomly chosen between 1
and 20 and the AR and MA coefficients were randomly chosen from $U(0.8,0.9)$.
The difference parameter $d$ was randomly chosen from $U(-0.5,0.5)$ for ARFIMA
processes.

It can be seen from the performance metrics in Table 3 that the model worked
extremely well for accurate classification of above discussed causal time-
series with distribution shifts and with the presence of causal influence that
was occurring at very different scales as compared to the influence in the
causal-time series used in the training set.

We also plot figures to illustrate the difference in TTSS features of the
amplitudes at different frequencies for causal and non-causal time series used
in this section. Figures 3(a), 3(b), 3(c) and 3(d) show the figures for causal
time series, AR(15), AR(100), ARMA and ARFIMA respectively. Figures 4(a) and
4(b) show the figures for non-causal time series generated using
$\mathcal{N}(0,0.01)$ and $U(-0.6,0.6)$ respectively. Clearly the number of
peaks and valleys in TTSS feature plots of causal time-series are much lower
than those in TTSS feature plots of non-causal time series.

(a) (b)

(c) (d)

Figure 3: TTSS feature representation for the amplitude at each frequency for
causal time-series: (a)AR(15), (b)AR(100), (c)ARMA and (d)ARFIMA

(a) (b)

Figure 4: TTSS feature representation for the amplitude at each frequency for
non-causal time-series generated from: (a) $\mathcal{N}(0,0.01)$ and (b)
$U(-0.6,0.6)$.

Hyperparameters used for ChaosFEX feature extraction for all datasets in this
section are $q=0.33$, $b=0.499$, $\varepsilon=0.01$ and maximum trajectory
length of firing for each GLS neuron was set to 1000. This means that the GLS
neuron would stop firing after a maximum trajectory length of 1000 if it did
not reach the $\varepsilon$ neighborhood of the stimulus. In [22] and [23],
the authors have found the hyperparameters tuned to $q=0.34$, $b=0.499$ and
$\varepsilon$ in the range $0.005$ to $0.3$, to be the most effective for the
classification of sequential data such as genome sequences and spoken digits.
Thus we fine tuned the hyperparameters by limiting the exploration of the
parameters around this range. Hyperparameter tuning for the LR classifier was
also done in this case. All parameters were set to default values other than
_‘ $max\\_iter$’, ‘tol’_ and _‘C’_ , which were set to $1000,0.001$ and
$0.001$ respectively.

##  5 Discussion, concluding remarks and future research directions

We show that ML classifiers (Logistic regression and LSTM), when used by
themselves directly on time-series measurements are dumb to the temporal/
causal-structure in the data. This fact has also been discussed in existing
literature [1, 2]. When time-series values were directly passed to the
classifiers, LR and LSTM, they failed to learn any causal-structure
characteristics from the data. Even though LSTM has been developed for the
classification of sequential data, it seemed to be learning some statistical
features from the data, giving high classification accuracy when the testing
set followed the same distribution as the training set and failing when there
was a distribution shift in the non-causal time series structure.

Even though characteristics of the fourier transformed time series of causal
and non-causal type are quite different, classifiers LR and LSTM, when
directly used on the frequency domain amplitudes fail to learn the causal
structure. Though they could accurately classify testing set which was exactly
the same as training set, the classification accuracy dropped to about $50\%$,
when there was a change in the distribution of non-causal time series. Almost
all of these time series were being wrongly classified as causal.

TTSS feature representation of the amplitude values at different frequencies
proved to be robust in learning the causal structure. We show that it works in
case of distribution shifts in both causal and non-causal time series and
could learn a generalized structure, identifying the processes correctly
irrespective of the scale at which causal-influence existed in the processes.
Good performance of the developed pipeline on our simple problem illustrates
the power of chaos and of ‘neurochaos’ based hybrid learning architectures for
developing more sophisticated causality based ML architectures as well as in
using ML for causal inference tasks. Since, even a simple linear classifier,
LR, is able to distinguish between TTSS features from causal and non-causal
time series, the strength of this feature seems to lie in the ability to
transform the given data to make it linearly separable. Hence, the developed
pipeline seems to be helpful in strong feature extraction for causal-structure
learning purposes.

Future work would involve more rigorous testing of the TTSS based pipeline
with other datasets having different levels of noise, strength of causal
coefficient and different causal/ non-causal structures. Simulated data with
‘interventional perturbations’ will also be provided to the algorithms to test
their performance on such data. Understanding and classification of these kind
of cause-effect treatment datasets will benefit the most from the application
of the proposed approach. Using classifiers other than LR, for example Support
Vector Machine with a linear kernel, to classify based on the TTSS model will
also be done. Other features based on the chaotic firing trajectory of GLS
neurons in the Chaosnet architecture have also been utilized in hybrid ML
based learning tasks [11, 22]. We would also like to check the performance of
the ChaosFEX based model using features other than the TTSS feature (or GLS
firing rate) for classification of causal and non-causal time-series.
Experiments and analysis will be done to get a better theoretical
understanding of how and why the ChaosFEX model works well for causal-
structure learning tasks. Finally, we would like to apply the developed
technique to recognize real time-series data with causal structure. We would
like to check the performance of the method for cause-effect treatment
estimation tasks for real data, both when the treatment/intervention being
provided is a discrete event or a continuous event. There are still not many
well developed techniques available for the latter purpose [5] and the
proposed pipeline seems like a promising approach.

## Acknowledgment

The authors are thankful to Harikrishnan N.B., National Institute of Advanced
Studies, for providing help with the use of ChaosFEX toolbox. N. Nagaraj
gratefully acknowledges the financial support of Tata Trusts and Dept. of
Science and Tech., Govt. of India (grant no. DST/CSRI/2017/54). A. Kathpalia
acknowledges the financial support of the Czech Science Foundation, Project
No. GA19-16066S and the Czech Academy of Sciences, Praemium Academiae awarded
to M. Paluš.

## References

  * [1] Judea Pearl and Dana Mackenzie.  The book of why: the new science of cause and effect.  Basic Books, 2018.
  * [2] Bernhard Schölkopf.  Causality for machine learning.  arXiv preprint arXiv:1911.10500, 2019.
  * [3] Erica EM Moodie, Thomas S Richardson, and David A Stephens.  Demystifying optimal dynamic treatment regimes.  Biometrics, 63(2):447–455, 2007.
  * [4] Susan Athey and Guido W Imbens.  The state of applied econometrics: Causality and policy evaluation.  Journal of Economic Perspectives, 31(2):3–32, 2017.
  * [5] Raha Moraffah, Paras Sheth, Mansooreh Karami, Anchit Bhattacharya, Qianru Wang, Anique Tahir, Adrienne Raglin, and Huan Liu.  Causal inference for time series analysis: Problems, methods and evaluation.  arXiv preprint arXiv:2102.05829, 2021.
  * [6] Siyang Leng, Ziwei Xu, and Huanfei Ma.  Reconstructing directional causal networks with random forest: Causality meeting machine learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(9):093130, 2019.
  * [7] Yu Huang, Zuntao Fu, and Christian LE Franzke.  Detecting causality from time series in a machine learning framework.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(6):063116, 2020.
  * [8] James Durbin and Geoffrey S Watson.  Testing for serial correlation in least squares regression: I.  Biometrika, 37(3/4):409–428, 1950.
  * [9] Greta M Ljung and George EP Box.  On a measure of lack of fit in time series models.  Biometrika, 65(2):297–303, 1978.
  * [10] Trevor S Breusch.  Testing for autocorrelation in dynamic linear models.  Australian Economic Papers, 17(31):334–355, 1978.
  * [11] NB Harikrishnan and Nithin Nagaraj.  Neurochaos inspired hybrid machine learning architecture for classification.  In 2020 International Conference on Signal Processing and Communications (SPCOM), pages 1–5. IEEE, 2020.
  * [12] Harikrishnan Nellippallil Balakrishnan, Aditi Kathpalia, Snehanshu Saha, and Nithin Nagaraj.  Chaosnet: A chaos based artificial neural network architecture for classification.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(11):113125, 2019.
  * [13] Robert H Shumway, David S Stoffer, and David S Stoffer.  Time series analysis and its applications, volume 3.  Springer, 2000.
  * [14] Hans Von Storch and Francis W Zwiers.  Statistical analysis in climate research.  Cambridge university press, 2001.
  * [15] George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung.  Time series analysis: forecasting and control.  John Wiley & Sons, 2015.
  * [16] Clive WJ Granger and Roselyne Joyeux.  An introduction to long-memory time series models and fractional differencing.  Journal of time series analysis, 1(1):15–29, 1980.
  * [17] Timothy Graves, Robert Gramacy, Nicholas Watkins, and Christian Franzke.  A brief history of long memory: Hurst, mandelbrot and the road to arfima, 1951–1980.  Entropy, 19(9):437, 2017.
  * [18] Peter McCullagh and John A Nelder.  Generalized linear models.  Routledge, 2019.
  * [19] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [20] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  the Journal of machine Learning research, 12:2825–2830, 2011.
  * [21] Francois Chollet et al.  Keras, 2015.
  * [22] Harikrishnan NB, Pranay SY, and Nithin Nagaraj.  A neurochaos learning architecture for genome classification.  arXiv preprint arXiv:2010.10995, 2020.
  * [23] NB Harikrishnan and Nithin Nagaraj.  When noise meets chaos: Stochastic resonance in neurochaos learning.  Neural Networks, 143:425–435, 2021.


# Frequency Limited $\mathcal{H}_{2}$ Optimal Model Reduction of Large-Scale
Sparse Dynamical Systems

Xin Du  School of Mechatronic Engineering and Automation, Shanghai University,
Shanghai-200072, China and Key Laboratory of Modern Power System Simulation
and Control & Renewable Energy Technology, Ministry of Education(Northeast
Electric Power University), Jilin-132012, China<EMAIL_ADDRESS>M. Monir
Uddin Department of Mathematics and Physics, North south University,
Dhaka-1229, Bangladesh<EMAIL_ADDRESS>A. Mostakim Fony
Department of Mathematics, Chittagong University, Chittagong, Bangladesh,
<EMAIL_ADDRESS>Md. Tanzim Hossain Department of Electrical and
Computer Engineering, North South University, Dhaka-1229, Bangladesh,
<EMAIL_ADDRESS>Mohammaed Sahadat-Hossain Department of
Mathematics and Physics, North south University, Dhaka-1229, Bangladesh,
<EMAIL_ADDRESS>
###### Abstract

We mainly consider the frequency limited $\mathcal{H}_{2}$ optimal model order
reduction of large-scale sparse generalized systems. For this purpose we need
to solve two Sylvester equations. This paper proposes efficient algorithm to
solve them efficiently. The ideas are also generalized to index-1 descriptor
systems. Numerical experiments are carried out using Python Programming
Language and the results are presented to demonstrate the approximation
accuracy and computational efficiency of the proposed techniques.

keywords Frequency limited model reduction, $\mathcal{H}_{2}$ optimal
condition, frequency limited Graminas, Sylvester and Lyapunov equations

##  1 Introduction

Model order reduction (MOR) is a process to approximate a high-order dynamical
system by a substantially low-order system with a maximum accuracy. This tool
is now widely used in different disciplines of science, engineering and
technology to reduce the complexity of the model. In general, the reduced
order models are used in controller design, simulation and optimization. For
motivation, applications and techniques of MOR see, e.g,. [1, 2].

The commonly useful methods for the model reduction of large-scale linear time
invariant dynamical systems are the balanced truncation and $\mathcal{H}_{2}$
optimal model reduction [1]. Both the methods are well established and
successfully investigated to find the model reduction of large-scale sparse
dynamical systems. In the recent times frequency and time limited model
reduction methods have taken a lot of attentions due to its demand in real-
life applications. In many applications, a specific frequency interval is more
important, i.e., the ROM should maintain a superior accuracy within that
desired frequency interval. Balanced truncation based frequency limited model
reduction was discussed by Gawronski and Juang in [3]. The computational
techniques for the time and frequency limited balanced truncation were
discussed in [4].

The optimal $\mathcal{H}_{2}$ model reduction methods have been studied and
investigated in [5, 6, 7, 8, 9]. See, also the references cited therein. In
all these papers the proposed technique is based on either Gramian assistance
first-order optimality conditions [5] or the tangential interpolaton [6] of
the transfer function. In fact both the conditions are coincided which is
shown in [9]. These papers only discuss the model reduction on the infinite
frequency interval. For the time limited case the we refer the readers to
[10]. Although, in [11] authors briefly introduced optimal $\mathcal{H}_{2}$
model reduction problem of standard state space systems considering a
restricted frequency interval, there the implementation details were not
given.

This paper focuses on the computational techniques of the frequency limited
optimal $\mathcal{H}_{2}$ model reduction method of large-scale sparse
systems. We mainly generalized the idea as in [12, 9] in which the proposed
algorithm is called two sided iteration algorithm (TSIA). Moreover, to
implement the frequency limited TSIA we need to solve two frequency limited
Sylvester equations. This paper also discusses how to solve the Sylvester
equations efficiently by preserving the sparsity of the system. Besides the
generalized systems the idea is also extended for index-1 descriptor systems.
The benefits of the algorithmic improvements presented in this paper are
illustrated by several numerical examples. We have generated the the results
by using Python Programming Language. Rest of this paper is organized as
follows.

Section 2 overview the TSIA and the optimal $\mathcal{H}_{2}$ model reduction
of generalized system. Then the ideas of this section are discussed in the
next sections for the frequency-limited model order reduction. Section 5
presents the algorithm for solving frequency limited Sylvesters equations
which provides projectors to carry out the FLMOR. The results of the numerical
experiments are depicted in Section 6 which show the efficiency and capability
of the proposed methods.

##  2 The TSIA and $\mathcal{H}_{2}$ optimal model order reduction of
generalized Systems

The goal of this section is to review the basic idea of $\mathcal{H}_{2}$
optimal model order reduction of generalized Systems. Let us consider a linear
time invariant continuous-time system of the form

$\displaystyle E\dot{x}(t)$ $\displaystyle=Ax(t)+Bu(t),$ (1) $\displaystyle
y(t)$ $\displaystyle=Cx(t)+Du(t),$

where $E\in\mathbb{R}^{n\times n}$ is non-singular, and
$A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times p}$,
$C\in\mathbb{R}^{m\times n}$ and ${D}\in\mathbb{R}^{m\times p}$. The transfer-
function matrix of this system is defined by $\text{G}(s)=C(sE-A)^{-1}B+D$,
where $s\in\mathbb{C}$. The controllability and the observability Gramians of
the system on the infinite frequency range can be defined as [13]

$\displaystyle P$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}(\imath\omega
E-A)^{-1}BB^{T}(\imath\omega E^{T}-A^{T})\mathop{d}\omega\quad\text{and}$
$\displaystyle Q$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}(\imath\omega
E^{T}-A^{T})^{-1}C^{T}C(\imath\omega E-A)\mathop{d}\omega,$

and they are the solutions of the continuous-time algebraic Lyapunov equations

$\displaystyle APE^{T}+EPA^{T}+BB^{T}=0\quad\text{and}$ (2) $\displaystyle
A^{T}QE+E^{T}QA+C^{T}C=0,$

respectively. We want to construct a substantially reduced-order model

$\displaystyle\dot{\hat{x}}(t)$ $\displaystyle=\hat{A}\hat{x}(t)+\hat{B}u(t),$
(3) $\displaystyle\hat{y}(t)$ $\displaystyle=\hat{C}\hat{x}(t)+\hat{D}u(t),$

where $\hat{A}\in\mathbb{R}^{r\times r}$, $\hat{B}\in\mathbb{R}^{r\times p}$,
$\hat{C}\in\mathbb{R}^{m\times r}$ and $\hat{D}\in\mathbb{R}^{m\times p}$. The
goal is to minimize the error

$\displaystyle\Xi=\|\text{G}-\hat{\text{G}}\|_{\mathcal{H}_{2}},$ (4)

where $\|.\|_{\mathcal{H}_{2}}$ denotes the system’s $\mathcal{H}_{2}$-norm
[1] and $\hat{G}(s)=\hat{C}(sI-\hat{A})^{-1}\hat{B}+\hat{D}$ is the transfer-
function matrix of the reduced system. The $\mathcal{H}_{2}$-norm of the error
system as defined in (4) can be measured by

$\displaystyle\Xi=:\sqrt{\mathrm{Tr}\left(C_{\Xi}P_{\Xi}C_{\Xi}^{T}\right)}\quad\text{or}\quad\sqrt{\mathrm{Tr}\left(B_{\Xi}^{T}Q_{\Xi}B_{\Xi}\right)},$
(5)

where $P_{\Xi}$ and $Q_{\Xi}$ are the solutions of the Lyapunov equations

$\displaystyle
A_{\Xi}P_{\Xi}E_{\Xi}^{T}+E_{\Xi}P_{\Xi}A_{\Xi}^{T}+B_{\Xi}B_{\Xi}^{T}=0\quad\text{and}$
(6a) $\displaystyle
A_{\Xi}^{T}Q_{\Xi}E_{\Xi}+E_{\Xi}^{T}Q_{\Xi}A_{\Xi}+C_{\Xi}^{T}C_{\Xi}=0,$
(6b)

in which

$\displaystyle E_{\Xi}=\begin{bmatrix}E&\\\ &\hat{I}\end{bmatrix},\quad
A_{\Xi}=\begin{bmatrix}A&\\\ &\hat{A}\end{bmatrix},\quad
B_{\Xi}=\begin{bmatrix}B\\\ \hat{B}\end{bmatrix}\quad\text{and}\quad
C_{\Xi}=\begin{bmatrix}C&-\hat{C}\end{bmatrix}.$

Now, partitioning $P_{\Xi}$ and $Q_{\Xi}$ as

$\displaystyle P_{\Xi}=\begin{bmatrix}P&M\\\
M^{T}&\hat{P}\end{bmatrix}\quad\text{and}\quad Q_{\Xi}=\begin{bmatrix}Q&N\\\
N^{T}&\hat{Q}\end{bmatrix},\,\text{respectively}$

and plugging into (6) we obtain the following algebraic matrix equations

$\displaystyle APE^{T}+EPA^{T}+BB^{T}$ $\displaystyle=0,$ (7a)
$\displaystyle\hat{A}\hat{P}+\hat{P}\hat{A}^{T}+\hat{B}\hat{B}^{T}$
$\displaystyle=0,$ (7b) $\displaystyle AM+EM\hat{A}^{T}+B\hat{B}^{T}$
$\displaystyle=0,$ (7c) $\displaystyle A^{T}QE+E^{T}QA+C^{T}C$
$\displaystyle=0,$ (7d)
$\displaystyle\hat{A}^{T}\hat{Q}+\hat{Q}\hat{A}+\hat{C}^{T}\hat{C}$
$\displaystyle=0,$ (7e) $\displaystyle A^{T}N+E^{T}N\hat{A}-C^{T}\hat{C}$
$\displaystyle=0,$ (7f)

where $\hat{P}$ and $\hat{Q}$ are respectively known as the controllability
and observability Gramians of the reduced systems. Therefore, the
$\mathcal{H}_{2}$ norm of the error system in (5) can be measured by

$\Xi=:\begin{cases}\sqrt{{\mathrm{Tr}\left(CPC^{T}\right)}+{\mathrm{Tr}\left(\hat{C}\hat{P}\hat{C}^{T}\right)}+2{\mathrm{Tr}\left(CM\hat{C}^{T}\right)}}\\\
\quad\text{or}\quad\\\
\sqrt{{\mathrm{Tr}\left(B^{T}QB\right)}+{\mathrm{Tr}\left(\hat{B}^{T}\hat{Q}\hat{B}\right)}+2{\mathrm{Tr}\left(B^{T}N\hat{B}\right)}}.\end{cases}$
(8)

The first-order optimality conditions for the optimal $\mathcal{H}_{2}$ model
reduction was given in [14], which is known as Wilson conditions. Wilson
conditions are based on the first derivatives of (4) with respect to
$\hat{A}$, $\hat{B}$ and $\hat{C}$ as follows:

$\displaystyle\nabla\Xi_{\hat{A}}=2(\hat{Q}\hat{P}+W^{T}EV),\,\nabla\Xi_{\hat{B}}=2(\hat{Q}\hat{B}+W^{T}E^{-1}B),\,\nabla\Xi_{\hat{C}}=2(\hat{C}\hat{P}-CV).$

Setting these three derivatives to zero leads to the _Wilson conditions_ ,

$\displaystyle\hat{Q}\hat{P}+N^{T}EM=0,$ (9)
$\displaystyle\hat{Q}\hat{B}+N^{T}E^{-1}B=0,$ (10)
$\displaystyle\hat{C}\hat{P}-CM=0.$ (11)

These three conditions in fact yield the left and and right projection
matrices to compute an optimal reduced order system (3) and in the optimal
reduced order system the reduced matrices are formed as

$\displaystyle\hat{A}=W^{T}E^{-1}AV,\quad\hat{B}=W^{T}B,\quad\,\text{and}\quad\hat{C}=CV,$
(12)

where $V=M\hat{P}^{-1}$ and $W^{T}=-\hat{Q}^{-1}N^{T}$ and hence it can be
proved that $W^{T}EV=I$. However, we can not guarantee that $\hat{P}$ and
$\hat{Q}$ are invertible, since to assure this the reduced model should be
completely controllable and observable [1]. In the case that they are
invertible, the multiplication from the right is only a transformation of
bases and does not change the subspace. The idea by Xu and Zeng [9] was to
satisfy the Wilson conditions by setting

$\displaystyle W=N\quad\text{and}\quad V=M.$ (13)

Note that $V$ and $W$ can be computed by solving the matrix equations (7c) and
(7f), respectively which are known as Sylvesters equations. Another important
observation is that if we want to compute the optimal projection subspace we
already need the optimal solution $\hat{A}$, $\hat{B}$ and $\hat{C}$. However,
this is not known prior. A possible solution is to start with a reduced model,
which emerged from an arbitrary projection of the original model, solve matrix
equations (7c) and (7c), compute the projectors, and restart the process with
the newly obtained reduced model until we are satisfied. In this way we get a
kind of a fixed point iteration. This procedure is called two sided iteration
algorithm (TSIA) by Xu and Zeng in [9].

The Wilson conditions are Gramian-based conditions since they are related to
Gramians of the system. The Hyland-Bernstein conditions [7] are another
gramian-based first-order optimal conditions, which were shown to be
equivalent to the Wilson conditions [15]. Van Dooren et al., were
characterized the tangential interpolation based $\mathcal{H}_{2}$ optimal
conditions in [12]. One drawback of interpolation based model reduction is to
selection of interpolation points. However in [15] authors proposed the
Iterative Rational Krylov Algorithms (IRKA) to resolve this problem. On the
other hand Xu and Zeng showed in [9] that both the Gramian and interpolation
based optimality conditions are the same. In [9] authors also presented the
two sided iterative algorithm (TSIA) for a standard system which is slightly
modified as in [12]. For our convenient the TSIA for the $\mathcal{H}_{2}$
optimal model reduction for the generalized system (1) is summarized in
Algorithm 1

Input : $E,A,B,C,D$.

Output : $\hat{A},\hat{B},\hat{C}$, $\hat{D}:=D$.

1 Choose matrices $W_{0}\in\mathbb{R}^{n\times r}$ and
$V_{0}\in\mathbb{R}^{n\times r}$ such that $W_{0}^{T}V_{0}=I$.

2 Construct the reduced-order matrices:
$\hat{A}=W_{0}^{T}E^{-1}AV_{0},\,\hat{B}=W_{0}^{T}E^{-1}B\quad\text{and}\quad\hat{C}=CV_{0}$
.

3 while  _$i\leq N-1$_ do

4 Compute $V_{i}$ and $W_{i}$ by solving Sylvesters  $\displaystyle
AV+EV\hat{A}^{T}+B\hat{B}^{T}=0$ (14a) $\displaystyle
A^{T}W+E^{T}W\hat{A}-C^{T}\hat{C}=0,$ (14b)

5 Compute $W_{i+1}=W_{i}(V_{i}^{T}W_{i})^{-1}$ and $V_{i+1}=V_{i}$

6 Construct the reduced-order matrices
$\hat{A}=W_{i+1}^{T}E^{-1}AV_{i+1},\hat{B}=W_{i+1}^{T}E^{-1}B$ and
$\hat{C}=CV_{i+1}$.

7 $i=i+1$.

8 end while

Algorithm 1 Two-sided iteration algorithm (TSIA).

##  3 FL $\mathcal{H}_{2}$ optimal MOR of generalized systems

This section moves to the frequency limited $\mathcal{H}_{2}$ optimal model
reduction of the system (1). For this purpose we first define frequency
limited Gramians. In the previous section the system Gramians have been
defined on the infinite frequency interval. If we replace the infinite
interval (-$\infty$, $\infty$) into a finite interval $\omega=[\omega_{1}$,
$\omega_{2}]$) then the controlability and the observablity Gramians can be
defined as

$\displaystyle P_{\omega}$
$\displaystyle=\frac{1}{2\pi}\int_{\omega_{1}}^{\omega_{2}}(\imath\nu
E-A)^{-1}BB^{T}(\imath\nu E^{T}-A^{T})\mathop{d}\nu$ (15) $\displaystyle
Q_{\omega}$
$\displaystyle=\frac{1}{2\pi}\int_{\omega_{1}}^{\omega_{2}}(\imath\nu
E^{T}-A^{T})^{-1}C^{T}C(\imath\nu E-A)\mathop{d}\nu,$ (16)

which satisfies the frequency limited Lyapunov equations

$\displaystyle
AP_{\omega}E^{T}+EP_{\omega}A^{T}+B_{\omega}BB^{T}+BB^{T}B_{\omega}^{\ast}=0,$
(17) $\displaystyle
A^{T}Q_{\omega}E+E^{T}Q_{\omega}A+C_{\omega}^{\ast}C^{T}C+C^{T}CC_{\omega}=0,$

with

$\displaystyle
B_{\omega}=\frac{\imath}{2\pi}\ln(A+\imath\omega_{2}E)(A+\imath\omega_{1}E)^{-1}\,\,\text{and}$
$\displaystyle
C_{\omega}=\frac{\imath}{2\pi}\ln(A+\imath\omega_{1}E)^{-1}(A+\imath\omega_{2}E).$

The goal in this paper is to constructed a reduced model
$\hat{G}=(\hat{A},\hat{B},\hat{C})$ from a given model $G=(E,A,B,C)$ that
minimizes the error

$\displaystyle\Xi_{\omega}=\|\text{G}-\hat{\text{G}}\|_{\mathcal{H}_{2,\omega}},$
(18)

where $\|.\|_{\mathcal{H}_{2,\omega}}$ denotes the $\mathcal{H}_{2}$-norm on
the prescribed frequency range $\omega$. The transfer-function matrix of the
reduced system is The $\mathcal{H}_{2}$-norm of the error system as defined in
(4) can be measured efficiently by

$\displaystyle\Xi_{\omega}=:\sqrt{\mathrm{Tr}\left(C_{\Xi}P_{\Xi,\omega}C_{\Xi}^{T}\right)}\quad\text{or}\quad\sqrt{\mathrm{Tr}\left(B_{\Xi}^{T}Q_{\Xi,\omega}B_{\Xi}\right)},$
(19)

where $P_{\Xi,\omega}$ and $Q_{\Xi,\omega}$ are the solutions of the Lyapunov
equations

$\displaystyle
A_{\Xi}P_{\Xi,\omega}E_{\Xi}^{T}+E_{\Xi}P_{\Xi,\omega}A_{\Xi}^{T}+B_{\Xi,\omega}B_{\Xi}B_{\Xi}^{T}+B_{\Xi}B_{\Xi}^{T}B_{\Xi,\omega}^{\ast}=0\quad\text{and}$
(20a) $\displaystyle
A_{\Xi}^{T}Q_{\Xi,\omega}E_{\Xi}+E_{\Xi}^{T}Q_{\Xi,\omega}A_{\Xi}+C_{\Xi,\omega}^{\ast}C_{\Xi}^{T}C_{\Xi}+C_{\Xi}^{T}C_{\Xi}C_{\Xi,\omega}=0,$
(20b)

where

$\displaystyle
B_{\Xi,\omega}=\frac{\imath}{2\pi}\ln(A_{\Xi}+\imath\omega_{2}E_{\Xi})(A_{\Xi}+\imath\omega_{1}E_{\Xi})^{-1}\,\,\text{and}$
$\displaystyle
C_{\Xi,\omega}=\frac{\imath}{2\pi}\ln(A_{\Xi}+\imath\omega_{1}E_{\Xi})^{-1}(A_{\Xi}+\imath\omega_{2}E_{\Xi}).$

Due to the structure of $E_{\Xi}$, $A_{\Xi}$, $B_{\Xi}$ and $C_{\Xi}$ we can
Partition $P_{\Xi,\omega}$, $Q_{\Xi,\omega}$, $B_{\Xi,\omega}$ and
$C_{\Xi,\omega}$ as follows

$\displaystyle P_{\Xi,\omega}=\begin{bmatrix}P_{\omega}&M_{\omega}\\\
M_{\omega}^{T}&\hat{P}_{\omega}\end{bmatrix},\quad
Q_{\Xi,\omega}=\begin{bmatrix}Q_{\omega}&N_{\omega}\\\
N_{\omega}^{T}&\hat{Q}_{\omega}\end{bmatrix},$ $\displaystyle
B_{\Xi,\omega}=\begin{bmatrix}B_{\omega}&0\\\
0&\hat{B}_{\omega}\end{bmatrix},\qquad
C_{\Xi,\omega}=\begin{bmatrix}C_{\omega}&0\\\
0&\hat{C}_{\omega}\end{bmatrix}.$

Therefore (20) yields

$\displaystyle
AP_{\omega}E^{T}+EP_{\omega}A^{T}+B_{\omega}BB^{T}+BB^{T}B_{\omega}^{\ast}=0,$
(21a)
$\displaystyle\hat{A}\hat{P}_{\omega}+\hat{P}_{\omega}\hat{A}^{T}+\hat{B}_{\omega}\hat{B}\hat{B}^{T}+\hat{B}\hat{B}^{T}\hat{B}_{\omega}^{\ast}=0,$
(21b) $\displaystyle
AM_{\omega}+EM_{\omega}\hat{A}^{T}+B_{\omega}B\hat{B}^{T}+B\hat{B}^{T}\hat{B}_{\omega}^{\ast}=0,$
(21c) $\displaystyle
AQ_{\omega}E^{T}+EQ_{\omega}A^{T}+C_{\omega}^{\ast}C^{T}C+C^{T}CC_{\omega}=0,$
(21d)
$\displaystyle\hat{A}\hat{Q}_{\omega}+\hat{Q}_{\omega}\hat{A}^{T}+\hat{C}_{\omega}^{\ast}\hat{C}^{T}\hat{C}+\hat{C}^{T}\hat{C}\hat{C}_{\omega}=0,$
(21e) $\displaystyle
AN_{\omega}+EN_{\omega}\hat{A}^{T}+C_{\omega}^{\ast}C^{T}\hat{C}-C^{T}\hat{C}\hat{C}_{\omega}=0,$
(21f)

with

$\displaystyle\hat{B}_{\omega}=\frac{\imath}{2\pi}\ln(\hat{A}+\imath\omega_{2}I)(\hat{A}+\imath\omega_{1}I)^{-1}\,\,\text{and}$
$\displaystyle\hat{C}_{\omega}=\frac{\imath}{2\pi}\ln(\hat{A}+\imath\omega_{1}I)^{-1}(\hat{A}+\imath\omega_{2}I).$

Following the discussion in the above section here we also construct reduced
order model by constructing the reduced matrices as in (12). We solve the the
Sylvesters equations (21c) and (21f) to construct $V=M_{\omega}$ and
$W=N_{\omega}$. The constructed reduced order model is $\mathcal{H}_{2}$
optimal on the limited frequency range and satisfies Wilson’s first-order
optimality conditions. The whole procedure is summarized in Algorithm 2. The
main computation tasks in this algorithm is to solve sparse-dense Sylvesters
equations (22a) and (22b). Following section will presents how to solve them
efficiently.

Input : $E,A,B,C$.

Output : $\hat{A},\hat{B},\hat{C}$, $\hat{D}:=D$.

1 Choose matrices $W_{0}\in\mathbb{R}^{n\times r}$ and
$V_{0}\in\mathbb{R}^{n\times r}$ such that $W_{0}^{T}V_{0}=I$.

2 Construct the reduced-order matrices
$\hat{A}=W_{0}^{T}E^{-1}AV_{0},\hat{B}=W_{0}^{T}E^{-1}B$ and $\hat{C}=CV_{0}$.

3 while  _( $i\leq N-1$)_ do

4 Compute $V_{i}=M_{\omega}$ and $W_{i}=N_{\omega}$ by solving the Sylvester
$\displaystyle
AM_{\omega}+EM_{\omega}\hat{A}^{T}+B_{\omega}B\hat{B}^{T}+B\hat{B}^{T}\hat{B}_{\omega}^{\ast}=0,$
(22a) $\displaystyle
AN_{\omega}+EN_{\omega}\hat{A}^{T}+C_{\omega}^{\ast}C^{T}\hat{C}-C^{T}\hat{C}\hat{C}_{\omega}=0,$
(22b) Compute $W_{i+1}=W_{i}(V_{i}^{T}W_{i})^{-1}$ and $V_{i+1}=V_{i}$.

5 Construct the reduced-order matrices
$\hat{A}=W_{i+1}^{T}E^{-1}AV_{i+1},\hat{B}=W_{i+1}^{T}E^{-1}B$ and
$\hat{C}=CV_{i+1}$.

6 $i=i+1$.

7 end while

Algorithm 2 Two-sided iteration algorithm (TSIA).

##  4 Solution of semi-generalized Sylvester equations

Above section shows that to perform the frequency limited model reduction of
system (1) we need to solve two frequency limited matrix equations namely
Sylvester equations (22a) and (22b). This section discusses how to solve them
efficiently. Since the Sylvester equations (22a) and (22b) are duel of each
other we only interested to elaborate the solution of (22a). Another one can
be solved by applying the same procedure. For our convenient we rewrite the
equation (22a) as

$\displaystyle AX+EX\hat{A}^{T}+F=0,$ (23)

where $F=B_{\omega}BB^{T}+BB^{T}B_{\omega}^{\ast}$ and $X=M_{\omega}$. The
technique that we have followed here was presented in [16] where $E=I$ is an
identity matrix. In [17] authors generalized the idea of [16] for the equation
like (23) where $F=B\hat{B}$.

Considering the _Schur decomposition_ of $\hat{A}$ as $QSQ^{\ast}$ such that
$QQ^{\ast}=Q^{\ast}Q=I$ and inserting this into (23) we get

$\displaystyle AX+EXQSQ^{\ast}+F=0.$ (24)

By multiplying this from the right with $Q$ we obtain

$\displaystyle
A\underbrace{XQ}_{\tilde{X}}+E\underbrace{XQ}_{\tilde{X}}S+\underbrace{FQ}_{\tilde{F}}=0,$
(25)

Observing that $S$ is a upper triangular matrices leads to a formula for the
first column of $\tilde{X}$:

$\displaystyle A\tilde{X}_{1}+E\tilde{X}_{1}S_{1,1}+\tilde{F}_{1}=0,$ (26)
$\displaystyle\Leftrightarrow\quad$
$\displaystyle(A+S_{11}E)\tilde{X}_{1}=-\tilde{F}_{1}.$ (27)

For all other columns we have to take care of the linear combination of $E$
matrix. If we consider the second column of the solution

$\displaystyle(A+S_{22}E)\tilde{X}_{2}=-\tilde{F}_{2}-S_{12}E\tilde{X}_{1}.$
(28)

In this way the arbitrary column $j$ of $\tilde{X}$ we find

$\displaystyle(A+S_{jj}E)\tilde{X}_{j}=-\tilde{F}_{j}-E\sum_{i=1}^{j-1}S_{ij}\tilde{X}_{i}.$
(29)

To obtain the solution of the original system we multiply $\tilde{X}$ by
$U^{\ast}$ from the right.

Input : $E,A,\hat{A},F$ from (23).

Output : $X\in\mathbb{R}^{n\times r}$ solution of (23).

1 Compute the Schur decomposition $\hat{A}=QSQ^{\ast}$ and Define
$\tilde{F}=FQ$

2 for  _$i=1,\cdots,r$_ do

3 Compute $\hat{F}=-\tilde{F}_{j}-E\sum_{i=1}^{j-1}S_{i,j}\tilde{X}_{i}$

4 Solve $(A+S_{jj}E)\tilde{X}_{j}=\hat{F}$

5 end for

$X=\tilde{X}Z^{\ast}$.

Algorithm 3 Solution of semi-generalized Sylvester equations.

##  5 FLMOR of structured index-1 systems

The index 1 descriptor system that we consider in the section has the
following form.

$\begin{array}[]{rcl}E_{1}\dot{x}(t)&=J_{1}x(t)+J_{2}z(t)+B_{1}u(t)\\\
0&=J_{3}x(t)+J_{4}z(t)+B_{2}u(t)\\\
y(t)&=C_{1}x(t)+C_{2}z(t)+D_{a}u(t),\end{array}$ (30)

Where $x(t)\in\mathbb{R}^{n_{1}}$ is the vector with differential variables
and $z(t)\in\mathbb{R}^{n_{2}}$ is the vector with algebraic variables. Model
reduction of such descriptor system has been discussed in a couple of previous
research articles, e.g., [18, 19, 20, 21, 22, 20, 23, 24] on an unrestricted
frequency limit. In [18] author uses spectral projector to split the system
into finite and infinite parts and balancing based method was applied to the
finite part. Other papers implemented MOR without computing the spectral
projector, rather eliminating the algebraic part the system was converted into
an ODE system. However, practical implementation was carried out without
computing the ODE system explicitly. This paper generalizes this idea for of
the FLMOR discussed in the above section.

By eliminating the algebraic variables i.e., $z(t)\in\mathbb{R}^{n_{2}}$ of
the system we obtain a generalized system (1) where the coefficients matrices
are defined as

$\displaystyle E:=E_{1},\quad A:=J_{1}-J_{2}{J_{4}}^{-1}J_{3},\quad
B:=B_{1}-J_{2}{J_{4}}^{-1}B_{2},$ (31) $\displaystyle
C:=C_{1}-C_{2}{J_{4}}^{-1}J_{3},\quad D:=D_{a}-C_{2}{J_{4}}^{-1}B_{2}.$

The index-1 and generalized systems are equivalent since the responses of the
systems are same and their finite eigenvalues are coincided. Such structured
system are arising in power system model [25]. For the FLMOR of index-1 system
we define $V$ and $W$ by solving the corresponding Sylvesters equations of the
generalized system as discussed in Section 3. Now, applying these
transformations the reduced system matrices can be constructed as:

$\displaystyle\hat{A}:=\hat{J}_{1}-\hat{J}_{2}{J}_{4}^{-1}\hat{J}_{3},\quad\hat{B}:=\hat{B}_{1}-\hat{J}_{2}\hat{J}_{4}^{-1}B_{2},$
(32)
$\displaystyle\hat{C}:=\hat{C}_{1}-C_{2}J_{4}^{-1}\hat{J}_{3},\quad\hat{D}:=D_{a}-C_{2}{J}_{4}^{-1}B_{2},$

where $\hat{J}_{1}=W^{T}E_{1}^{-1}J_{1}V$, $\hat{J}_{2}=W^{T}J_{2}$,
$\hat{J}_{3}=J_{3}V$, $\hat{B}_{1}=W^{T}E_{1}^{-1}B_{1}$ and
$\hat{C}_{1}=C_{1}V$. To compute the the projectors $V$ and $W$ by solving the
corresponding Sylvesters equations is a challenging task since the input
matrices in (32) are highly dense. In the following text we discuss how to
solve the Sylvesters equations related to the index-1 system (30) efficiently.

To solve the Sylvesters equations we can use Algorithm 3. In this algorithm
the main expensive task is to solve a linear system at each iteration step. At
Step 4 of the algorithm we need to solve the linear system

$\displaystyle(A+S_{ii}E)\tilde{X}_{j}=\hat{F}.$

Plugging $A$ and $E$ from (31) we obtain

$\displaystyle(J_{1}-J_{2}J_{4}^{-1}J_{3}+S_{ii}E_{1})\tilde{X}_{j}=\hat{F},$

which can be rewritten as

$\displaystyle(J_{1}+S_{ii}E_{1}-J_{2}J_{4}^{-1}J_{3})\tilde{X}_{j}=\hat{F}.$
(33)

A close observation reveal that instead of this we can solve the following
linear system

$\displaystyle\begin{bmatrix}J_{1}+S_{ii}E_{1}&J_{2}\\\
J_{3}&J_{4}\end{bmatrix}\begin{bmatrix}\tilde{X}_{j}\\\
\Gamma\end{bmatrix}=\begin{bmatrix}\hat{F}\\\ 0\end{bmatrix}$ (34)

for $\tilde{X}_{j}$. Although, linear system in (34) is larger than the system
in (33) it is sparse and hence can be solved by any sparse solvers (e.g.,
direct [26, 27] or iterative e.g., [28, 29]) efficiently.

Input : $E_{1},J_{1},J_{2},J_{3},J_{4},B_{1},B_{2},\hat{A}$.

Output : $X\in\mathbb{R}^{n\times r}$.

1 Form $F=B_{\omega}B\hat{B}^{T}+B\hat{B}^{T}\hat{B}_{\omega}^{\ast}$, where
$\displaystyle
B_{\omega}=\frac{\imath}{2\pi}\ln(J_{1}+\imath\omega_{2}E_{1}-J_{2}J_{4}^{-1}J_{3})(J_{1}+\imath\omega_{1}E_{1}-J_{2}J_{4}^{-1}J_{3})^{-1},$
$\displaystyle\hat{B}_{\omega}=\frac{\imath}{2\pi}\ln(\hat{A}+\imath\omega_{1}I)^{-1}(\hat{A}+\imath\omega_{2}I),\quad\text{and}\quad
B=B_{1}-J_{2}J_{4}^{-1}B_{2}$ Compute the Schur decomposition
$\hat{A}=QSQ^{\ast}$ and Define $\tilde{F}=FQ$

2 for  _$i=1,\cdots,r$_ do

3 Compute $\hat{F}=-\tilde{F}_{j}-E_{1}\sum_{i=1}^{j-1}S_{i,j}\tilde{X}_{i}$

4 Solve  $\displaystyle\begin{bmatrix}J_{1}+S_{ii}E_{1}&J_{2}\\\
J_{3}&J_{4}\end{bmatrix}\begin{bmatrix}\tilde{X}_{j}\\\
\Gamma\end{bmatrix}=\begin{bmatrix}\hat{F}\\\ 0\end{bmatrix}$ for
$\tilde{X}_{j}$

5 end for

$X=\tilde{X}Z^{\ast}$.

Algorithm 4 Sylvester equation for index-1 system.

##  6 Numerical results

To asses the efficiency of our proposed techniques in this section we discuss
the numerical results. For our convenient we have splitted the section into
several subsections.

###  6.1 Model examples

We consider the following model examples for the numerical experiments.

######  Example 6.1 (International Space Station (ISS)).

This is a model of stage 1R (Russian Service Module) of the ISS. It has n=270
states, p=3 inputs and m=3 outputs. The details of the model can be obtain in
[30].

######  Example 6.2 (Clamped beam model (CBM)).

This structural model was obtained by spatial discretization of an appropriate
partial differential equation (see [31]). The dimension of the model is n=348
and it is a single input single output i.e., SISO system. The input represents
the force applied to the structure at the free end, and the output is the
resulting displacement.

######  Example 6.3 (Triple chain oscillator (TCO) model).

Although this example was originated in [32], the setup was described in [33]
which resulted in second-order system. We convert into firs-order form in
which the dimension of the system is n=10000. It is also an SISO system and
input, output matrices are transpose of each other.

######  Example 6.4 (Power system model).

We consider several Brazilian Power System (BPS) models from [21] which are in
index-1 descriptor form. Table 1 shows number of differential ($n_{1}$) and
algebraic $n_{2}$ variables and inputs/outputs of several models which are
used for the numerical experiments here.

model name | $n_{1}$ | $n_{2}$ | $m/p$
---|---|---|---
BPS-606 | 606 | 1142 |
BPS-1142 | 1142 | 8593 | 4/4
BPS-1450 | 1450 | 9855 |
BPS-1693 | 1693 | b11582 |
Table 1: Dimension of differential and algebraic variables of different
Brazilian power system models.

###  6.2 Setup Hardware and Software

The experiments are carried out with Python 3.7.9 on a board with AMD
$\text{Ryzen}^{\text{TM}}$ $\text{Threadripper}^{\text{TM}}$ 1920X 12-Core
Processor with a 3.5 GHz clock speed and 128 GB RAM.

###  6.3 Error analysis of reduce-order model

The proposed techniques was applied to the all model examples mentioned above.
For the ISS, CBM, TCO we apply Algorithm 2 to obtain frequency limited reduced
order model. On the other for the BPS models we have applied the techniques
discussed in Section 5. We have computed different dimensional reduced order
models for different model examples which are mentioned in Table 2. The table
also shows $\mathcal{H}_{2}$ norm of the error systems in both frequency
restricted and unrestricted cases. For all the models frequency restricted
reduced order model show much more better accuracy than the frequency
frequency underused ones within the assigned frequency intervals.

Model |  $r$ |  $\Xi_{\omega}$ |  $\Xi$
---|---|---|---
ISS |  30 |  7.3244$\times 10^{-8}$ |  2.8000$\times 10^{-3}$
TCO |  30 |  3.4850$\times 10^{-4}$ |  9.6000$\times 10^{-3}$
BPS-606 |  30 |  7.8992$\times 10^{-4}$ |  1.0600$\times 10^{-2}$
BPS-1142 |  35 |  3.4421$\times 10^{-5}$ |  5.3000$\times 10^{-3}$
BPS-1450 |  35 |  9.3818$\times 10^{-8}$ |  1.5063$\times 10^{-5}$
BPS-1693 |  45 |  4.2995$\times 10^{-6}$ |  1.3000 $\times 10^{-3}$
Table 2: Dimension of reduced order models and $\mathcal{H}_{2}$ norms of the
error systems with and without limited frequency intervals.

We also have investigated the frequency domain analysis including the errors
of the original and reduced order models by using sigma plots. Exemplary, we
have depicted the frequency responses of TCO and BPS models only. Figures 1
and 2 show comparisons of the frequency responses of the of the original and
the reduced order models of TCO and BPS models, respectively . In both the
figures from absolute and the relative errors we observe that frequency
restricted reduced order models approximate with the original models with
higher accuracy within on the prescribed frequency ranges.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$\displaystyle\omega$$\sigma{}_{\text{max}}\text{(G(j}\omega\text{))}$Full
modelFrequency unrestrictedFrequency restricted (a) Sigma plot

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$10^{-4}$$10^{-1}$$10^{2}$$\displaystyle\omega$$\sigma{}_{\text{max}}\text{(G(j}\omega\text{)-}\hat{\text{G}}\text{(j}\omega\text{))}$
(b) Absolute error

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\displaystyle\omega$$\frac{\sigma{}_{\text{max}}\text{(G(j}\omega\text{)-}{\hat{\text{G}}}\text{(j}\omega\text{))}}{\sigma{}_{\text{max}}\text{(G(j}\omega\text{))}}$
(c) Relative error

Figure 1: Comparison of original and 30 dimensional reduced systems on the
frequency range [1,2] for TCO.

$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$16$$10^{1.5}$$10^{2}$$\displaystyle\omega$$\sigma{}_{\text{max}}\text{(G(j}\omega\text{))}$Full
modelFrequency unrestrictedFrequency restricted (a) Sigma plot

$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$16$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$\displaystyle\omega$$\sigma{}_{\text{max}}\text{(G(j}\omega\text{)-}\hat{\text{G}}\text{(j}\omega\text{))}$
(b) Absolute error

$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$16$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\displaystyle\omega$$\frac{\sigma{}_{\text{max}}\text{(G(j}\omega\text{)-}{\hat{\text{G}}}\text{(j}\omega\text{))}}{\sigma{}_{\text{max}}\text{(G(j}\omega\text{))}}$
(c) Relative error

Figure 2: Comparison of original and 45 dimensional reduced systems on the
frequency range [6,10] for BPS-1693.

###  6.4 Comparison of sparse and dense system

We already mentioned that in the TSIA the main computational task is solving
two Sylvester equations. This paper presents Algorithm 4 to solve the
Sylvester equations efficiently for index-1 system. We know that by converting
index-1 system into a generalized system we can solve the Sylvester equation
using Algorithm 3. Figure 3 shows the time comparisons of Algorithms 3 and 4
for solving Sylvester equations of index-1 system. We see that if the
dimension of the system is increased then the computational times for
Algorithms 3 is increased rapidly. On the other hand the computational time
with Algorithm 4 is nominal in compare to Algorithms 3.

BPS-606BPS-1142BPS-1450BPS-1693$0$$1$$2$$3$$4$$5$$6$$7$$8$0.621.943.197.740.530.730.871.16time(sec)Algorithm
3Algorithm 4 Figure 3: Comparison of computational time to solve the Sylvester
equation for different dimensional BPS model using Algorithms 3 and 4.

##  7 Conclusions

In this paper we have discussed the frequency limited $\mathcal{H}_{2}$
optimal model order reduction of large-scale sparse dynamical systems. For
this purpose we have solved two mixed (Sparse and dense) Sylvester equations
which are dual of each other. We have shown that how to solve the Sylvester
equations efficiently without loosing the sparsity of large sparse system
matrices. The ideas are also generalized to index-1 descriptor systems.
Index-1 system can be converted into a generalized system by eliminating the
algebraic equations, which however convert the system from sparse to dense. We
have discussed how to perform the model reduction without converting the dense
system explicitly. Numerical experiments has been carried out to demonstrate
the approximation accuracy and computational efficiency of the proposed
algorithm using Python Programming Language.

#### Acknowledgments.

This research work was funded by NSU-CTRG research grant under the project
No.: CTRG-19/SEPS/05. It was also supported by National Natural Science
Foundation of China under Grant No. (61873336, 61873335), the Fundamental
Research Funds for the Central Universities under Grant (FRF-BD-19-002A), and
the High-end foreign expert program of Shanghai University,

## References

  * [1] A. Antoulas, _Approximation of Large-Scale Dynamical Systems_ , ser. Advances in Design and Control. Philadelphia, PA: SIAM Publications, 2005, vol. 6.
  * [2] M. M. Uddin, _Computational Methods for Approximation of Large-Scale Dynamical Systems_. New York, USA: Chapman and Hall/CRC, 2019.
  * [3] W. Gawronski and J. Juang, “Model reduction in limited time and frequency intervals,” _Int. J. Syst. Sci._ , vol. 21, no. 2, pp. 349–376, 1990.
  * [4] P. Benner, P. Kürschner, and J. Saak, “Frequency-limited balanced truncation with low-rank approximations,” _SIAM J. Sci. Comput._ , vol. 38, no. 1, pp. A471–A499, Feb. 2016.
  * [5] P. Van Dooren, K. Gallivan, and P.-A. Absil, “$\mathcal{H}_{2}$-optimal model reduction of MIMO systems,” _Appl. Math. Lett._ , vol. 21, pp. 1267–1273, 2008.
  * [6] S. Gugercin, A. C. Antoulas, and C. A. Beattie, “$\mathcal{H}_{2}$ model reduction for large-scale dynamical systems,” _SIAM J. Matrix Anal. Appl._ , vol. 30, no. 2, pp. 609–638, 2008.
  * [7] D. Hyland and D. Bernstein, “The optimal projection equations for model reduction and the relationships among the methods of wilson, skelton, and moore.”
  * [8] W. Y. Yan and J. Lam, “An approximate approach to $\mathcal{H}_{2}$ optimal model reduction,” _IEEE Trans. Autom. Control_ , vol. 44, no. 7, pp. 1341–1358, 1999.
  * [9] Y. Xu and T. Zeng, “Optimal $\mathcal{H}_{2}$ model reduction for large scale MIMO systems via tangential interpolation,” _International Journal of Numerical Analysis and Modeling_ , vol. 8, no. 1, pp. 174–188, 2011.
  * [10] P. Goyal and M. Redmann, “Time-limited $\mathcal{H}_{2}$-optimal model order reduction,” _Applied Mathematics and Computation_ , vol. 355, pp. 184–197, 2019.
  * [11] D. Petersson and J. Löfberg, “Model reduction using a frequency-limited $\mathcal{H}_{2}$-cost,” _Systems & Control Letters_, vol. 67, pp. 32–39, 2014.
  * [12] P. Van Dooren, K. A. Gallivan, and P. A. Absil, “H2-optimal model reduction of mimo systems,” _Appl. Math. Lett._ , vol. 21, no. 12, pp. 1267–1273, 2008\.
  * [13] K. Zhou, G. Salomon, and E. Wu, “Balanced realization and model reduction for unstable systems,” _Internat. J. Robust and Nonlinear Cont._ , vol. 9, no. 3, pp. 183–198, 1999.
  * [14] D. A. Wilson, “Optimum solution of model-reduction problem,” _Proceedings of the Institution of Electrical Engineers_ , vol. 117, no. 6, pp. 1161–1165, 1970.
  * [15] S. Gugercin, A. C. Antoulas, and C. Beattie, “H_2 model reduction for large-scale linear dynamical systems,” _SIAM journal on matrix analysis and applications_ , vol. 30, no. 2, pp. 609–638, 2008.
  * [16] D. C. Sorensen and A. C. Antoulas, “The sylvester equation and approximate balanced reduction,” _Numer. Lin. Alg. Appl._ , vol. 351–352, pp. 671–700, 2002.
  * [17] P. Benner, M. Köhler, and J. Saak, “Sparse-dense sylvester equations in $h_{2}$-model order reduction,” 2011.
  * [18] T. Stykel, “Analysis and numerical solution of generalized Lyapunov equations,” Ph.D. Thesis, Technische Universität Berlin, Berlin, 2002.
  * [19] S. Gugercin, T. Stykel, and S. Wyatt, “Model reduction of descriptor systems by interpolatory projection methods,” _SIAM J. Sci. Comput._ , vol. 35, no. 5, pp. B1010–B1033, 2013.
  * [20] M. M. Uddin, J. Saak, B. Kranz, and P. Benner, “Computation of a compact state space model for an adaptive spindle head configuration with piezo actuators using balanced truncation,” _Production Engineering_ , vol. 6, pp. 577–586, 2012.
  * [21] F. Freitas, J. Rommes, and N. Martins, “Gramian-based reduction method applied to large sparse power system descriptor models,” _IEEE Trans. Power Syst._ , vol. 23, no. 3, pp. 1258–1270, Aug. 2008.
  * [22] M. M. Uddin, “Model reduction for piezo-mechanical systems using Balanced Truncation,” Master’s thesis, Stockholm University, Stockholm, Sweden, 2011\. [Online]. Available: http://www.qucosa.de/fileadmin/data/qucosa/documents/7822/Master\\_Thesis\\_Uddin.pdf
  * [23] M. S. Hossain and M. M. Uddin, “Iterative methods for solving large sparse lyapunov equations and application to model reduction of index 1 differential-algebraic-equations,” _Numerical Algebra, Control & Optimization_, vol. 9, no. 2, p. 173, 2019.
  * [24] M. M. Uddin, M. S. Hossain, and M. F. Uddin, “Rational krylov subspace method (rksm) for solving the lyapunov equations of index-1 descriptor systems and application to balancing based model reduction,” in _9th International Conference on Electrical & Computer Engineering (ICECE) 2016_. IEEE, 2016, pp. 451–454.
  * [25] R. W. Freund, “Structure-preserving model order reduction of RCL circuit equations,” in _Model Order Reduction: Theory, Research Aspects and Applications_ , W. H. A. Schilders, H. A. van der Vorst, and J. Rommes, Eds. Springer-Verlag, 2008, pp. 49–73.
  * [26] T. A. Davis, _Direct Methods for Sparse Linear Systems_ , ser. Fundamentals of Algorithms. Philadelphia, PA, USA: SIAM, 2006, no. 2.
  * [27] I. S. Duff, A. M. Erisman, and J. K. Reid, _Direct methods for sparse matrices_. Oxford, UK: Clarendon Press, 1989.
  * [28] H. A. Van der Vorst, _Iterative Krylov Methods for Large Linear Systems_. Cambridge: Cambridge University Press, 2003.
  * [29] Y. Saad, _Iterative Methods for Sparse Linear Systems_. Philadelphia, PA, USA: SIAM, 2003.
  * [30] S. Gugercin, A. C. Antoulas, and N. Bedrossian, “Approximation of the international space station $1$r and $12$a models,” in _Proc. of the 40th IEEE Conferences on Decision and Control_ , 2001, pp. 1515–1516.
  * [31] A. C. Antoulas, D. C. Sorensen, and S. Gugercin, “A survey of model reduction methods for large-scale systems,” _Contemp. Math._ , vol. 280, pp. 193–219, 2001.
  * [32] N. Truhar and K. Veselić, “An efficient method for estimating the optimal dampers’ viscosity for linear vibrating systems using Lyapunov equation,” _SIAM J. Matrix Anal. Appl._ , vol. 31, no. 1, pp. 18–39, 2009.
  * [33] J. Saak, “Efficient numerical solution of large scale algebraic matrix equations in pde control and model order reduction,” Dissertation, TU Chemnitz, July 2009, available from http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642.