Linear Simple Cycle Reservoirs at the edge of stability perform Fourier decomposition of the input driving signals

Recurrent Neural Networks (RNNs) are input-driven parametric state-space machine learning models designed to capture temporal dependencies in sequential input data streams. They encode time series data into a latent state space, dynamically storing temporal information within state-space vectors.

Reservoir Computing (RC) models is a subset of RNNs that operate with a fixed, non-trainable input-driven dynamical system (known as the reservoir) and a static trainable readout layer producing model responses based on the reservoir activations. This design uniquely simplifies the training process by concentrating adjustments solely on the readout layer (thus avoiding back-propagating the error information backwards through time), leading to enhanced computational efficiency. The simplest implementations of RC models include Echo State Networks (ESNs) [Jae01, MNM02, TD01, LJ09].

Simple Cycle Reservoirs (SCR) represent a specialized class of RC models characterized by a single degree of freedom in the reservoir construction (modulo the state space dimensionality), structured through uniform ring connectivity and binary input weights with an aperiodic sign pattern. Recently, SCRs were shown to be universal approximators of time-invariant dynamic filters with fading memory over $\mathbb{C}$ and $\mathbb{R}$ respectively in [LFT24, FLT24], making them highly suitable for integration in photonic circuits for high-performance, low-latency processing [BVdSBV22, LBFM⁺17, HVK⁺20].

Understanding the intricacies of SCRs in depth is essential. In this work, we employ the kernel view of linear ESNs introduced in [Tin20], in which the state-space ‘reservoir’ representation of (potentially left-infinite) input sequences is treated as a feature map corresponding to the given reservoir model through the associated reservoir kernel. For linear reservoirs, the canonical dot product of two input sequences’ feature representations is analytically expressible as the (semi-)inner product of the sequences themselves. The corresponding metric tensor reveals the representational structure imposed by the reservoir on the input sequences, in particular in the metric tensor’s eigenspace containing dominant projection axes (time series ’motifs’) and scaling (‘importance’ factors).

To assess the “richness” of linear SCR state-space representations, [Tin20] proposed analyzing the relative area of motifs under the Fast Fourier Transform (FFT). It was observed that the richness of these representations collapses at the edge of stability when the spectral radius $\rho$ of the dynamic coupling matrix equals to 1.

In this paper we theoretically analyze the collapse of motif richness at the edge of stability and show that when $\rho=1$ the SCR kernel motifs correspond to Fourier basis. We begin by reviewing the notion of SCRs [RT10], kernel view of ESNs [Tin20], and Reservoir Motif Machines [TFL24] in Section 2. The contribution of this paper, in the subsequent sections, are outlined as follows:

1. (1)

   In Section 4, we show in $\mathbb{C}$ that motifs of linear SCR at the edge of stability are harmonic functions.

2. (2)

   In Section 5, we show in $\mathbb{R}$ that $n$ dimensional linear SCR has $\lceil\frac{n}{2}\rceil$ symmetric motifs and $\lfloor\frac{n}{2}\rfloor$ skew-symmetric motifs.

3. (3)

   In Section 6, we combine the results of the previous two sections and demonstrate numerically that in $\mathbb{R}$, the motifs of linear SCR at edge of stability are exactly the columns of real Fourier basis matrix.

4. (4)

   Finally in Section 7, we conclude the paper with numerical experiments supporting our findings.

Let $\mathbb{K}=\mathbb{R}$ or $\mathbb{C}$ be a field. We first formally define the principal object of our study - parametrized linear driven dynamical system with a (possibly) non-linear readout.

The idea of simple cycle reservoir was presented in [RT10] as a reservoir system with a very small number of design degrees of freedom, yet retaining performance capabilities of more complex or (unnecessarily) randomized constructions. In fact, it can be shown that even with such a drastically reduced design complexity, SCR models are universal approximators of fading memory filters [LFT24, FLT24].

One possibility to understand inner representations of the input-driving time series forming inside the reservoir systems is to view the reservoir state space as a temporal feature space of the associated reservoir kernel [Tin20, GGO22] . Consider a linear reservoir system $R=(\mbox{\bf{W}},\mbox{\bf{w}},h)$ over $\mathbb{K}$ with dimensions $(n,1,d)$ operating on univariate input.

Let $\tau>n$ denote the length of the look back window and consider two sufficiently long time series of length $\tau>n$,

	u	$\displaystyle=\left(u\left(-\tau+1\right),u\left(-\tau+2\right),\ldots,u\left(-1\right),u\left(0\right)\right)$
		$\displaystyle=:\left(u_{1},u_{2},\ldots,u_{\tau}\right)\in\mathbb{K}^{\tau}$

and

	v	$\displaystyle=\left(v\left(-\tau+1\right),v\left(-\tau+2\right),\ldots,v\left(-1\right),v\left(0\right)\right)$
		$\displaystyle=:\left(v_{1},v_{2},\ldots,v_{\tau}\right)\in\mathbb{K}^{\tau}$

we consider the reservoir states reached upon reading them (with zero initial state) their feature space representations [Tin20]:

$$\phi(\mbox{\bf{u}})=\sum_{j=1}^{\tau}u_{j}\mbox{\bf{W}}^{\tau-j}\mbox{\bf{w}},\ \ \ \phi(\mbox{\bf{v}})=\sum_{j=1}^{\tau}v_{j}\mbox{\bf{W}}^{\tau-j}\mbox{\bf{w}}.$$

The canonical dot product (reservoir kernel)

$$K(\mbox{\bf{u}},\mbox{\bf{v}})=\langle\phi(\mbox{\bf{u}}),\phi(\mbox{\bf{v}})\rangle$$

can be written in the original time series space as a semi-inner product $\langle\mbox{\bf{u}},\mbox{\bf{v}}\rangle_{Q}=\mbox{\bf{u}}^{\top}\mbox{\bf{Q}}\mbox{\bf{v}}$, where

(2.2)

$$Q_{i,j}=\mbox{\bf{w}}^{\top}\left({\mbox{\bf{W}}}^{\top}{}\right)^{i-1}\mbox{\bf{W}}^{j-1}\mbox{\bf{w}}.$$

Since the (semi-)metric tensor $\mbox{\bf{Q}}\in\mathbb{M}_{\tau\times\tau}(\mathbb{K})$ is symmetric and positive semi-definite, it admits the following eigen-decomposition:

(2.3)

$\displaystyle\mbox{\bf{Q}}=\mbox{\bf{M}}\Lambda_{Q}\mbox{\bf{M}}^{\top},$

where $\Lambda_{Q}:=\operatorname{diag}\left(\lambda_{1},\lambda_{2},\ldots,\lambda_{N_{m}}\right)$ is a diagonal matrix consisting of non-negative eigenvalues of Q with the corresponding eigenvectors $\mbox{\bf{m}}_{1},\mbox{\bf{m}}_{2},\ldots,\mbox{\bf{m}}_{N_{m}}\in\mathbb{K}^{\tau}$ (columns of M). The $N_{m}:=\operatorname{rank}\left(Q\right)\leq N\leq\tau$ eigenvectors of M with positive eigenvalues are called the motifs of $R$. We have:

$$K(\mbox{\bf{u}},\mbox{\bf{v}})=\left(\Lambda_{Q}^{\frac{1}{2}}\mbox{\bf{M}}^{\top}\mbox{\bf{u}}\right)^{\top}\left(\Lambda_{Q}^{\frac{1}{2}}\mbox{\bf{M}}^{\top}\mbox{\bf{v}}\right).$$

In particular, the reservoir kernel is a canonical dot product of time series projected onto the motif space spanned by $\left\{\mbox{\bf{m}}_{i}\right\}_{i=1}^{N_{m}}$:

$$K(\mbox{\bf{u}},\mbox{\bf{v}})=\langle\tilde{}\mbox{\bf{u}},\tilde{}\mbox{\bf{v}}\rangle,$$

where

(2.4)		$\displaystyle\tilde{}\mbox{\bf{u}}$	$\displaystyle=\Lambda_{Q}^{\frac{1}{2}}\mbox{\bf{M}}^{\top}\mbox{\bf{u}}$
		$\displaystyle=\begin{bmatrix}\lambda^{\frac{1}{2}}_{1}\cdot\langle\mbox{\bf{m}}_{1},\mbox{\bf{u}}\rangle\\ \vdots\\ \lambda^{\frac{1}{2}}_{N_{m}}\cdot\langle\mbox{\bf{m}}_{N_{m}},\mbox{\bf{u}}\rangle\end{bmatrix}$
		$\displaystyle=\left(\lambda_{i}^{\frac{1}{2}}\cdot\langle\mbox{\bf{m}}_{i},\mbox{\bf{u}}\rangle\right)_{i=1}^{N_{m}}\in\mathbb{M}_{N_{m}}(\mathbb{K}).$

Reservoir Motif Machine (RMM)[TFL24] is a predictive model motivated by the kernel view of linear Echo State Networks described above. By projecting the $\tau$-blocks of input time series onto the reservoir motif space given by $\operatorname{span}(\{\mbox{\bf{m}}_{i}\}_{i=1}^{N_{m}})$, RMM captures temporal and structural dynamics in a computationally efficient feature map.

In particular, rather than relying on the motif weights $\lambda_{i}^{\frac{1}{2}}$ determined by the reservoir $R$ in Equation (2.4), RMM introduces a set of adaptable motif coefficients, denoted as $C:=\{c_{i}\in\mathbb{R}\}_{i=1}^{N_{m}}$, to define its feature map as follows:

	$\displaystyle\varphi\left(\mbox{\bf{u}};C\right)$	$\displaystyle=\begin{bmatrix}c_{1}\cdot\langle\mbox{\bf{m}}_{1},\mbox{\bf{u}}\rangle\\ \vdots\\ c_{N_{m}}\cdot\langle\mbox{\bf{m}}_{N_{m}},\mbox{\bf{u}}\rangle\end{bmatrix}$
(2.5)			$\displaystyle=\left(c_{i}\cdot\langle\mbox{\bf{m}}_{i},\mbox{\bf{u}}\rangle\right)_{i=1}^{N_{m}}\in\mathbb{R}^{N_{m}}.$

This feature map is used to train a predictive model, such as linear regression or kernel-based methods, directly in the motif space.

Having defined the reservoir temporal kernel, one can ask how “representationally rich” is the associated feature space (span of the reservoir motifs). To quantify the ‘richness’ of the reservoir feature space of a linear reservoir system over $\mathbb{R}$, [Tin20] proposed the following procedure:

Consider a linear reservoir system $R=(\mbox{\bf{W}},\mbox{\bf{w}},h)$ with dimensions $(n,1,d)$ over $\mathbb{R}$. Suppose W has spectral radius $\rho$. Recall from Remark 2.4 that: to study the SCR motif structure of, it is sufficient to consider the past horizon $\tau=n$.The motif matrix $\mbox{\bf{M}}\in\mathbb{M}_{n\times n}\left(\mathbb{R}\right)$ is constructed according to Equation (2.3) from the matrix Q (metric tensor of the inner product of the reservoir kernel).

First, Fast Fourier Transform (FFT) is applied to the kernel motifs (columns of M), considering only those with motif weights upto a threshold of $10^{-2}$ of the highest motif weight. This yields an $n\times n^{\prime}$ matrix of Fourier coefficients over $\mathbb{C}$ with $n^{\prime}\leq n$. These Fourier coefficients are then collected in a (multi)set $Z:=\left\{z_{k}\right\}_{k=1}^{n\cdot n^{\prime}}$.

To evaluate the diversity and spread of the Fourier coefficients in the complex plane, [Tin20] proposed calculating the coarse-grained area occupied by $Z\subseteq\mathbb{C}$. In particular, the box $\left[-7,7\right]^{2}$ in the complex plane is partitioned into a grid of cells (following [Tin20] we use side length $0.05$). The relative area covered by $Z\subseteq\mathbb{C}$ is defined as the ratio of the number of cells visited by the coefficients $z_{k}\in Z$ to the total number of cells in the grid. An example of the distribution of Fourier coefficients $Z$ of linear SCR with $n=97$ at $\rho=0.9813,0.999,1$ and w being the first $n$ digits of binary expansion of $\pi$ is presented in Figure 1 ³³3In particular, $\rho=0.9812798473475446$ is where the relative area peaks at Figure 2..

Replicating the experiment in [Tin20] in Figure 2, we observe that the “richness” of the motif space of SCR, as measured by relative area under the current setup, increases as the spectral radius approaches approximately $\rho\approx 0.9813$. Beyond this point, the measure sharply declines, aligning with the results for a randomly generated reservoir. This decline is also observed in is also shown in Figure 2, as $\rho$ increases from approximately $0.98$ (the peak of Figure 2) to the edge of stability at $\rho=1$.

In this paper, we aim to investigate this collapse of representational richness of SCR models at the edge of stability $\rho=1$. In particular, we will show that at the edge of stability, the motifs of linear SCR are (sampled) harmonic functions.

The next two section will be a two-part exposition of the properties of the motifs of linear SCR:

1. (1)

   In Section 4, we show in $\mathbb{C}$ that the motifs are harmonic functions. Following the approach of [LFT24], we begin with unitary dynamical coupling and progress to cyclic permutations.

2. (2)

   In Section 5, we demonstrate in $\mathbb{R}$ that the number of symmetric motifs is $\lceil\frac{n}{2}\rceil$, while the number of skew-symmetric motifs is $\lfloor\frac{n}{2}\rfloor$. In line with [FLT24], we start with orthogonal coupling and then move onto cyclic permutation.

Combining these two results, we conclude that at $\rho=1$, the motifs alternate between real and imaginary components of the first Fourier basis, which correspond to cosines and sines, respectively. We supplement our theoretical findings with numerical simulations of the motif space of SCR in the next section, which then lead to the numerical experiments in the final section.

We first show that, in the complex domain $\mathbb{C}$, the motifs of SCR can be derived explicitly. In particular, in this section we set $\mathbb{K}=\mathbb{C}$ and show that when the spectral radius $\rho=1$, the motifs of linear SCR are harmonic, i.e. they are precisely the Fourier basis (columns of the Fourier matrix).

Consider a linear reservoir system $R=(\mbox{\bf{W}},\mbox{\bf{w}},h)$ over $\mathbb{C}$ with dimensions $(n,1,d)$ and W of spectral radius $\rho$. Let $\mbox{\bf{Q}}_{\rho}$ denote the metric tensor of the reservoir kernel (Eq. (2.2).

In the spirit of [LFT24], we begin by considering a more general setting of $\mbox{\bf{W}}=\rho\mbox{\bf{U}}$, where U is a unitary matrix in $\mathbb{M}_{n\times n}(\mathbb{C})$ (i.e. $\mbox{\bf{U}}\mbox{\bf{U}}^{*}=\mbox{\bf{U}}^{*}\mbox{\bf{U}}={\mbox{\bf{I}}}$) and $\rho\in(0,1]$. We then move to the special case where $\mbox{\bf{U}}=\mbox{\bf{C}}$ is a cyclic permutation. Since U is unitary, its eigenvalues all have norm $1$. We let its eigenvalues be $\{s_{j}\in\mathbb{C}:1\leq j\leq n\}$ with corresponding eigenvectors $\{\xi_{j}:1\leq j\leq n\}$, and we know each $|s_{j}|=1$.

By construction:

$$\mbox{\bf{Q}}_{\rho}=[Q_{\rho\ (i,j)}]_{1\leq i,j\leq\tau}=[\mathbf{w}^{*}\mbox{\bf{W}}^{*(i-1)}\mbox{\bf{W}}^{j-1}\mathbf{w}]_{1\leq i,j\leq\tau}.$$

The matrix $\mbox{\bf{Q}}_{\rho}$ is a $\tau\times\tau$ matrix. Denote

$$\mbox{\bf{X}}_{\rho}=[\mbox{\bf{W}}^{*(i-1)}\mbox{\bf{W}}^{j-1}]_{1\leq i,j\leq\tau},$$

and,

$$\hat{\mathbf{w}}=\begin{bmatrix}\mathbf{w}&&&\\ &\mathbf{w}&&\\ &&\ddots&\\ &&&\mathbf{w}\end{bmatrix}.$$

Notice that $\mbox{\bf{X}}_{\rho}$ is an $\tau n\times\tau n$ matrix while $\hat{\mathbf{w}}$ is an $\tau n\times\tau$ matrix. Then by construction, $\mbox{\bf{Q}}_{\rho}=\hat{\mathbf{w}}^{*}\mbox{\bf{X}}_{\rho}\hat{\mathbf{w}}$. Notice that since $\mbox{\bf{W}}=\rho\mbox{\bf{U}}$ and U is unitary, we can rewrite $\mbox{\bf{X}}_{\rho}$ as follows:

$$\mbox{\bf{X}}_{\rho}=[\rho^{i+j-2}\mbox{\bf{U}}^{j-i}]_{1\leq i,j\leq\tau}=\begin{bmatrix}{\mbox{\bf{I}}}&\rho\mbox{\bf{U}}&\rho^{2}\mbox{\bf{U}}^{2}&\cdots&\rho^{\tau-1}\mbox{\bf{U}}^{\tau-1}\\ \rho\mbox{\bf{U}}^{*}&\rho^{2}{\mbox{\bf{I}}}&\rho^{3}\mbox{\bf{U}}&&\vdots\\ \rho^{2}\mbox{\bf{U}}^{*2}&\rho^{3}\mbox{\bf{U}}^{*}&\rho^{4}{\mbox{\bf{I}}}&&\\ \vdots&&&\ddots&\\ \rho^{\tau-1}\mbox{\bf{U}}^{*(\tau-1)}&&&&\rho^{2(\tau-1)}{\mbox{\bf{I}}}\end{bmatrix}$$

Finally, we let $\lambda=1+\rho^{2}+\rho^{4}+\cdots+\rho^{2(\tau-1)}=\begin{cases}\tau,&\text{ if }\rho=1\\ \frac{1-\rho^{2\tau}}{1-\rho^{2}},&\text{ if }\rho<1\end{cases}.$

With the ultimate goal of computing the eigen-decomposition of $\mbox{\bf{Q}}_{\rho}$, we first characterize the eigenvalues and eigenvectors of $\mbox{\bf{X}}_{\rho}$. For the eigenvalues we first observe the following:

As a result, we can now fully characterize the eigenvalues of $\mbox{\bf{X}}_{\rho}$.

We now turn to characterize the non-zero eigenvectors of $\mbox{\bf{X}}_{\rho}$. We will then use these vectors to compute the motif matrix of $\mbox{\bf{Q}}_{\rho}$. Recall $\xi_{j}$ are the eigenvectors of U with the corresponding eigenvalues $s_{j}$. For each $j$, define

$$\hat{\xi}_{j}=\begin{bmatrix}\xi_{j}\\ \rho s_{j}^{-1}\xi_{j}\\ \vdots\\ \rho^{\tau-1}s_{j}^{-(\tau-1)}\xi_{j}\end{bmatrix}\in\mathbb{C}^{n\tau}$$

Since all the other eigenvalues of $\mbox{\bf{X}}_{\rho}$ are $0$, we immediately have the desired eigen-decomposition of the matrix $\mbox{\bf{X}}_{\rho}$ as:

$$\mbox{\bf{X}}_{\rho}=\frac{1}{\lambda}\sum_{j=1}^{n}\hat{\xi}_{j}\hat{\xi}_{j}^{*}$$

Recall an $n\times n$ full-cycle permutation matrix is given by:

$$\mbox{\bf{C}}:=\begin{bmatrix}0&\cdots&\cdots&\cdots&1\\ 1&0&&\cdots&0\\ 0&1&0&\vdots&\vdots\\ \vdots&0&\ddots&0&\vdots\\ 0&\vdots&0&1&0\end{bmatrix}$$

For the rest of the section, we set $\rho=1$ and $\mbox{\bf{W}}={\mbox{\bf{C}}}$ is a full-cycle permutation and we will now compute eigen-decomposition of the metric tensor $\mbox{\bf{Q}}=\mbox{\bf{Q}}_{1}$ of the corresponding reservoir kernel.

From elementary matrix analysis, we know the eigenvalues of a full-cycle permutation C are precisely the $n$-th root of unities $\{\omega_{j}=e^{\frac{2\pi ij}{n}},1\leq j\leq n\}$. Its normalized eigenvectors for each eigenvalue $\omega_{j}$ is given by the Fourier basis:

(4.1)

$$\xi_{j}=\frac{1}{\sqrt{n}}\begin{bmatrix}1\\ \omega_{j}\\ \omega_{j}^{2}\\ \vdots\\ \omega_{j}^{n-1}\end{bmatrix}\in\mathbb{C}^{n},\quad j=1,\ldots,n.$$

We will now compute the eigenvalues and eigenvectors of Q explicitly. First, we see that

$$\mbox{\bf{Q}}=\hat{\mathbf{w}}^{*}\mbox{\bf{X}}_{1}\hat{\mathbf{w}},$$

where $\mbox{\bf{X}}_{1}=\frac{1}{n}\sum_{j=1}^{n}\hat{\xi}_{j}\hat{\xi}_{j}^{*}$. Define $\xi_{j}^{*}\mathbf{w}=\langle\mathbf{w},\xi_{j}\rangle=d_{j}$.

While, following Remark 2.4, we focused on the case where $\tau=n$, the results of this section up to this point (Lemma 4.1, Corollary 4.2, and Lemma 4.3) hold for general $\tau$.

As a result, in the case when $\rho=1$ and the matrix $\mbox{\bf{W}}={\mbox{\bf{C}}}$ is a full-cycle permutation, the motif matrix for Q consists of $\xi_{j}$, which is precisely the discrete Fourier transform matrix. In other words, the columns of the motif matrix are precisely the Fourier basis.

In this section we return to the real domain, i.e. $\mathbb{K}=\mathbb{R}$ and show that at unit spectral radius, the motifs of SCR consists of a fixed number of symmetric and skew symmetric vectors. As in [FLT24], we begin by deriving properties of the motif space of linear reservoir systems with orthogonal dynamical coupling and then move on to the special case of cyclic permutation.

Let $\mbox{\bf{W}}\in\mathbb{M}_{n\times n}\left(\mathbb{R}\right)$ be the dynamical coupling matrix of a reservoir system $R=(\mbox{\bf{W}},\mbox{\bf{w}},h)$ over $\mathbb{R}$. Suppose W is orthogonal with spectral radius $\rho=1$. We now show that the matrix corresponding to the reservoir kernel $\mbox{\bf{Q}}:=\mbox{\bf{Q}}_{\rho=1}$ is Toeplitz.

Q is Toeplitz if and only if $\mbox{\bf{Q}}_{ij}=\mbox{\bf{Q}}_{i+1,j+1}$ for all $i,j=1,\ldots,n-1$. By construction of Q, this is satisfied if and only if:

$$\mbox{\bf{w}}^{\top}\left(\mbox{\bf{W}}^{\top}\right)^{i-1}\mbox{\bf{W}}^{j-1}\mbox{\bf{w}}=\mbox{\bf{w}}^{\top}\left(\mbox{\bf{W}}^{\top}\right)^{i}\mbox{\bf{W}}^{j}\mbox{\bf{w}}.$$

Now, orthogonality of W implies $\mbox{\bf{W}}^{\top}\mbox{\bf{W}}=\mbox{\bf{I}}$. Without loss of generality assume that $i\geq j$. Then

$$\left(\mbox{\bf{W}}^{\top}\right)^{i-1}\mbox{\bf{W}}^{j-1}=\left(\mbox{\bf{W}}^{\top}\right)^{i-j}=\left(\mbox{\bf{W}}^{\top}\right)^{i}\mbox{\bf{W}}^{j},$$

showing that Q is indeed Toeplitz.

Let J denote the exchange matrix with $1$ on the antidiagonal and $0$ everywhere else:

$$\mbox{\bf{J}}:=\begin{bmatrix}0&0&0&\cdots&0&1\\ 0&0&0&\cdots&1&0\\ 0&0&\cdots&1&0&0\\ \vdots&\vdots&\iddots&\iddots&\vdots&\vdots\\ 0&1&0&\cdots&0&0\\ 1&0&0&\cdots&0&0\end{bmatrix}$$

By [CB76]: symmetric centrosymmetric matrices of order n admits an orthonormal basis of eigenvectors with:

1. (1)

   $\lceil\frac{n}{2}\rceil$ symmetric eigenvectors,

2. (2)

   $\lfloor\frac{n}{2}\rfloor$ skew-symmetric eigenvectors

This section bridges the theoretical results and numerical experiments by explicitly constructing the Fourier basis matrix corresponding to a linear SCR over $\mathbb{R}$. We present numerical simulations to validate our theoretical findings and outline the components of the numerical experiments in the final section.

Combining the results of the previous sections, we conclude that: At unit spectral radius, the motifs of linear SCR over $\mathbb{R}$ with $n$ neurons of look back window $\tau=n$ are:

1. R.1

   Harmonic, i.e. they correspond to the Fourier basis (Theorem 4.4).

2. R.2

   There are $\lceil\frac{n}{2}\rceil$ symmetric vectors (cosines) and $\lfloor\frac{n}{2}\rfloor$ skew-symmetric vectors (sines) with positive eigenvalues. (Theorem 5.4).

We now demonstrate the explicit construction of the Fourier basis matrix corresponding to the motifs of a linear SCR at unit spectral radius.

For a linear reservoir system over $\mathbb{R}$, the motifs must also be real. By conditions R.1 and R.2, each motif must then correspond to either the real or imaginary parts of the Fourier basis. More precisely, for $k=0,\ldots,\lceil\frac{n}{2}\rceil$ and $i=0,\ldots,\tau-1$ (The ceiling function accounts for Theorem 5.4.), the real Fourier basis matrix $\mbox{\bf{F}}=\mbox{\bf{F}}[i,j]\in\mathbb{M}_{\tau\times n}(\mathbb{R})$ is defined as

(6.1)

$\displaystyle\begin{split}\mbox{\bf{F}}[i,2k]&=\sqrt{\frac{2}{\tau}}\cos{\frac{2\pi ki}{\tau}}\quad\ldots\quad\text{Even columns}\\ \mbox{\bf{F}}[i,2k+1]&=\sqrt{\frac{2}{\tau}}\sin{\frac{2\pi ki}{\tau}}\quad\ldots\quad\text{Odd columns}\end{split}$

We claim that the motifs of a linear SCR are exactly the first $n$ columns of F.

In Figure 3 and Figure 4, we present the Fourier analysis of motifs of the SCR at unit spectral radius when $\tau=n=97$.

In Figure 4, we present examples of $4$ randomly chosen motifs and their corresponding Fourier basis in F. While some motifs and their corresponding Fourier basis are shifted by a fixed phase of $\pi$ or reflected over the $x$-axis, their Fourier spectra align as illustrated in Figure 3. Furthermore, Figure 3, we observe that the eigenvalues come in pairs, as discussed in 4.6. The motifs may not be strictly harmonic in the classical sense due to the coarseness of the discretization grid (non-integer division of period) (See Remark 6.1 below).

We conclude the paper with numerical experiments illustrating our findings. We compare time series forecasting results using the so-called Reservoir Motif Machines (RMM) [TFL24] and the linear SCR. RMM is a simple time series forecasting method based on the feature space representation of time series derived from linear reservoirs and has demonstrated remarkable predictive performance, even surpassing that of more complex transformer models on univariate time series forecasting tasks. A brief exposition of RMM is presented in Section 2.

We use RMM because it allows us to explicitly define the feature space representation by imposing a set of motifs, rather than having the feature space representation implicitly defined as in classical ESNs, such as SCR. This provides an ideal platform to showcase our theoretical findings: that the feature space representation of the motif space of a linear SCR is the same as that defined by the Fourier basis matrix F in Equation 6.1.

For reproducibility of the experiments, all experiments are CPU-based and are performed on Apple M3 Max with 128GB of RAM. The source code and data of the numerical analysis is openly available at https://github.com/Lampertos/motif_Fourier.

We compare the prediction results on univariate time series forecasting across the following models:

1. (1)

   Lin-RMM with SCR motifs and unit spectral radius,

2. (2)

   Lin-RMM with Fourier basis motifs as discussed in Section 6, and

3. (3)

   Linear SCR with unit spectral radius.

It is essential to emphasize that this experiment is not intended to showcase the predictive capability of the models, but rather to highlight the similarities in feature space representation between the linear SCR model over $\mathbb{R}$ and the Fourier basis motifs introduced in Section 6. Consequently, hyperparameters are kept constant throughout the experiments to underscore this feature space comparison.

The fixed set of hyperparameters for all experiments are as follows: $r_{\text{in}}=0.05$ for input weights, $n=97$ for the number of reservoir neurons, and $\tau=2n=194$ for the look-back window length. All models are trained using ridge regression with a ridge coefficient of $10^{-3}$. The prediction horizon is set to $168$.

Linear recurrent neural networks (RNN), such as Echo State Networks (ESN) can be thought of as providing feature representations of the input-driving time series in their state space [Tin20, GGO22]. By endowing the feature (state) space with the canonical dot-product, one can “reverse engineer” the corresponding inner product in the space of time series (time series kernel) that is defined through the the dot product in the RNN state space [Tin20]. This in turn helps to shed light on the inner representational schemes employed by the RNN to process the input-driving time series. In particular, the induced (semi-)inner product in the time series space can be theoretically analyzed through eigen-decomposition of the corresponding metric tensor. The eigenvectors (time series motifs) define the projection basis of the induced feature space and the (decay of) eigenvalues its dominant subspace and effective dimensionality. The induced time series kernels by the Simple Cycle Reservoir (SCR) models were shown to be superior (in terms of dimensionality, motif variability, and memory) to several alternative ESN constructions [Tin20].

In this paper we have shown a rather surprising result: When SCR is constructed at the edge of stability, the basis of its induced time series feature space correspond to the well known and widely used basis for signal decomposition - namely the Fourier basis.

This insight also explains the reduction in relative area covered by Fourier representations of SCR motifs observed by [Tin20] at the edge of stability. Our results imply that the feature space representation of a linear SCR at unit spectral radius effectively performs a weighted projection onto the Fourier basis.

This observation is supported by numerical experiments, in which we compared the time series forecasting accuracy of Lin-RMM with the motif space defined by a linear SCR at unit spectral radius and the Fourier basis, respectively.