Computing Semilinear Sparse Models for Approximately Eventually Periodic Signals

Given an AEP signal $\Sigma=\{x_{t}\}_{t\geq 1}\subset\mathbb{C}$ and a lag value $L>0$. For the study reported in this document we will consider semilinear signal models of the form:

where $\mathcal{L}(\mathbf{x}_{L}(t))$ denotes a linear operation defined by the expression

the term $\mathcal{G}(\mathbf{x}_{L}(t))$ represents a linear combination of neural networks based on gated recurrent unit (GRU) layers, whose structure is described by the block diagram

and $\mathcal{E}(\mathbf{x}_{L}(t))$ represents some suitable error term. For some given integer $m>0$, each GRU cell $\mathbf{G}_{j}$ in the GRU block in (3) is described for each $j=1,\ldots,m$ by the following equations:

with $\mathbf{h}(t)=\begin{bmatrix}h_{1}(t)&\cdots&h_{m}(t)\end{bmatrix}^{\top}$, and where $\sigma$ denotes the sigmoid function. The configuration considered in (3) is based on a GRU layer in order to prevent vanishing gradients, by taking advantage of the GRU structure presented in Cho et al. (2014). Although LSTM networks can be used to prevent vanishing gradients as well, for the study reported in this document GRU networks were chosen instead of LSTM networks, as they have fewer trainable parameters. The affine layer $\mathbf{A}$ of the neural network described in (3) is determined by the expression

In order to compute models of the form (1), we can combine sparse autoregressive models of the form (2) that can be computed using the methods presented in Vides (2021b), with GRU neural network models of the form (3) that can be computed using the computational tools provided as part of TensorFlow, Keras and PyTorch, that are described as part of Chollet et al. (2015) and Paszke et al. (2019).

An approximate representation

of the linear component of (1) such that

can be computed using some sample $\Sigma_{N}=\{x_{t}\}_{t=1}^{N}$ and a corresponding subsample $\Sigma_{0}=\{x_{t}\}_{t=1}^{N-1}\subset\Sigma_{N}$ for some suitable $N>L$, by approximately solving the matrix equation

using sparse least-squares approximation methods.

Schematically, the semilinear autoregressors considered in this study can be described by a block diagram of the form,

where the block $\mathfrak{L}$ is represented by (2), each block $\mathfrak{G}_{j}$ is represented by (3), and where the block $\mathfrak{M}$ is a mixing block defined by the expression

for some coefficients $w_{j}$ to be determined and some given $N$, with $y_{1}(t)=\mathfrak{L}(\mathbf{x}_{L}(t))$ and $y_{k+1}(t)=\mathfrak{G}_{k}(\mathbf{x}_{L}(t))$ for each $k=1,\ldots,N$ and each $t\geq L$.

The details of the computation of the neural network blocks of model (1) will be omitted for brevity, for details on the theory and computation of the GRU neural network models considered for this study the reader is kindly referred to Cho et al. (2014), Chollet et al. (2015), Paszke et al. (2019) and Vides (2021a).

Several interesting papers have been written on the subject of hybrid time series models that combine ARIMA and ANN models, two important references on this matter are Zhang (2003) and Khandelwal et al. (2015). Besides using sparse AR models instead of ARIMA models, another important distinctive aspect of the modeling approach reported in this document, is that instead of using the recurrent neural network components of (6) represented by $\mathcal{G}$ in (1) to approximate the residuals $r_{t}=x_{t+1}-\mathcal{L}(\mathbf{x}_{L}(t))$. Using some suitable training subsets $\Sigma_{I},\Sigma_{M}$ of a given data sample $\Sigma_{N}$ from an arbitrary AEP signal $\Sigma=\{x_{t}\}_{t\geq 1}$ under consideration, first the parameters of the blocks $\mathfrak{L}$, $\mathfrak{G}_{1}$, $\cdots$, $\mathfrak{G}_{N}$ of (6) are fitted independently using $\Sigma_{I}$, and then the coefficients of the block $\mathfrak{M}$ of (6) are fitted using $\Sigma_{M}$ and some corresponding predicted values generated with $\mathfrak{L}$, $\mathfrak{G}_{1}$, $\cdots$, $\mathfrak{G}_{N}$. Computing sparse representations of some of the matrix parameters of the resulting model along the process.

Given an AEP signal $\Sigma=\{x_{t}\}_{t\geq 1}$ whose behavior can be approximately described by a model of the form (2), that can be computed by approximately solving an equation of the form (5) for some suitable integers $N,L>0$ with $N>L$, and given some sample $\Sigma_{0}=\{x_{t}\}_{t=1}^{N-1}$, if we consider any sample $\tilde{\Sigma}_{0}=\{\tilde{x}_{t}\}_{t=1}^{N-1}\subset\Sigma$, such that for some positive integer $S$ the states in $\tilde{\Sigma}_{0}$ satisfy the conditions $\tilde{x}_{t}=x_{t+S}$, for each $t=1,\ldots,N-1$. We will have that the matrix

will approximately satisfy the condition

Using matrices of the form (7) one can express linear models of the form (2) as follows.

One can observe that to each model of the form (2), there corresponds a matrix of the form (7). From here on, a matrix that satisfies the previous conditions will be called the matrix form of a linear model $\mathcal{L}$ of the form (2).

Given $\delta>0$, and two matrices $A\in\mathbb{C}^{m\times n}$ and $Y\in\mathbb{C}^{m\times p}$, let us write $AX\approx_{\delta}Y$ to represent the problem of finding $X\in\mathbb{C}^{n\times p}$, $\alpha,\beta\geq 0$ and an orthogonal projector $Q$ such that $\|AX-Y\|_{F}\leq\alpha\delta+\beta\|(I_{m}-Q)Y\|_{F}$. The matrix $X$ will be called a solution to the problem $AX\approx_{\delta}Y$.

Since $\mathrm{rk}_{\delta}(\mathcal{H}_{L}(\Sigma_{N})^{\top})=\mathrm{rk}_{\delta}(\mathcal{H}_{L}(\Sigma_{N}))>0$ by (Vides, 2021b, Lemma 3.2). This result is a consequence of the application of (Vides, 2021b, Theorem 3.6) to the problem $\mathcal{H}_{L}(\Sigma_{N})^{\top}\hat{A}\approx_{\delta}\mathcal{H}_{L}(\Sigma_{N})^{\top}A$.

Given an AEP signal $\Sigma=\{x_{t}\}_{t\geq 1}$ with approximate period $T$ and a signal model $\mathcal{S}$ for $\Sigma$ of the form (1), if for the linear component $\mathcal{L}$ of $\mathcal{S}$ the corresponding residuals $r_{t}=|x_{t+1}-\mathcal{L}(\mathbf{x}_{L}(t))|$ are small, then the significative contribution of $\mathcal{L}$ to the modeling process of $\Sigma$, will be beneficial for interpretability purposes.

Let us consider an AEP signal $\Sigma$ with approximate period $T$ and approximately periodic tail $\tilde{\Sigma}=\{x_{t}\}_{t\geq S}$. Given some suitable integer lag value $L>0$ such that there is an approximate linear model $\mathcal{L}$ for $\Sigma$ of the form (2) with small residuals $r_{t}=|x_{t+1}-\mathcal{L}(\mathbf{x}_{L}(t))|$. Let $C_{L}$ be the matrix form of $\mathcal{L}$ and let $s=S+L-1$. By applying a Krylov subspace approach along the lines presented in (Saad, 2011, §6.1), and as a consequence of (Vides, 2021b, Theorem 4.3.), one can find a matrix $W_{k}\in\mathbb{C}^{L\times k}$ whose columns form an orthonormal basis of the subspace $\mathcal{K}_{T}\subset\mathbb{C}^{L}$ determined by the expression

such that each $z\in\sigma(W_{k}^{\ast}C_{L}W_{k})$ satisfies the relation $|z^{T}-1|\leq\varepsilon$, for some $\varepsilon>0$. The matrix $W_{k}^{\ast}C_{L}W_{k}$ will be called an approximately periodic (AP) $\Sigma$-section of $C_{L}$ and will be denoted by $C_{L}|_{\Sigma}^{AP}$.

In this section, algorithm 1 is applied to compute a SpARS model for the signal data sample recorded in the csv file VortexSheddingSignal.csv in the DataSets folder in Vides (2021a). The graphic representations of the results produced by the command line

>>> NumericalExperiment(1)

are shown in figures 1 and 2, respectively.

In every figure like figure 2, the red dots represent the points in each considered spectrum, the blue lines represent $\mathbf{S}^{1}$, and the number $T$ represents the estimated approximate period for each signal considered.

In this section, algorithm 1 is applied to compute a SpARS model for the signal data sample recorded in the csv file NLOscillatorSignal.csv in the DataSets folder in Vides (2021a). The graphic representations of the results produced by the command line

>>> NumericalExperiment(2)

are shown in figures 3 and 4, respectively.

In this section, algorithm 1 is applied to compute a SpARS model for the signal data sample recorded in the csv files:

• •

  art_daily_no_noise.csv

• •

  art_daily_small_noise.csv

that are included as part of the datasets described in Ahmad et al. (2017). The graphic representations of the results produced by the command line

>>> NumericalExperiment(3.1)

for the periodic signal without noise are shown in figures 5 and 6, respectively. The graphic representations of the results produced by the command line

>>> NumericalExperiment(3.2)

for the periodic signal with noise are shown in figures 7 and 8.

The approximation root mean square errors (RMSE) are summarized in table 1.

It is appropriate to mention that the root mean square errors can present little fluctuations as one performs several numerical simulations, due primarily to the nature of the neural-network models, as the linear components tend to present very low or no variability from simulation to simulation.

The Python programs that support the findings of this study are openly available in the SPAAR repository, reference number Vides (2021a). The time series data used for the experiments 1 and 2 documented in §5.1 and §5.2, respectively, are available as part of Vides (2021a). The time series data used for experiment 3 in §5.3 are available as part of the Numenta Anomaly Benchmark (NAB) described in Ahmad et al. (2017).