# Quality of internal representation shapes learning performance in feedback
neural networks

Lee Susman Equal contribution, listed in random order Interdisciplinary
program in Applied Mathematics, Technion Israel Institute of Technology,
Haifa, 32000, Israel  Network Biology Research Laboratories, Technion Israel
Institute of Technology, Haifa, 32000, Israel,
<EMAIL_ADDRESS>Francesca Mastrogiuseppe ††footnotemark:  Gatsby
Computational Neuroscience Unit, University College London, London
<EMAIL_ADDRESS>Naama Brenner Equal contribution Network
Biology Research Laboratories, Technion Israel Institute of Technology, Haifa,
32000, Israel,
<EMAIL_ADDRESS>Dept. of Chemical Engineering, Technion Israel Institute
of Technology, Haifa, 32000, Israel,
<EMAIL_ADDRESS>Omri Barak ††footnotemark:  Network Biology Research
Laboratories, Technion Israel Institute of Technology, Haifa, 32000, Israel,
<EMAIL_ADDRESS>Rappaport Faculty of Medicine, Technion Israel Institute
of Technology, Haifa, 32000, Israel,
<EMAIL_ADDRESS>
A fundamental feature of complex biological systems is the ability to form
feedback interactions with their environment. A prominent model for studying
such interactions is reservoir computing, where learning acts on low-
dimensional bottlenecks. Despite the simplicity of this learning scheme, the
factors contributing to or hindering the success of training in reservoir
networks are in general not well understood. In this work, we study non-linear
feedback networks trained to generate a sinusoidal signal, and analyze how
learning performance is shaped by the interplay between internal network
dynamics and target properties. By performing exact mathematical analysis of
linearized networks, we predict that learning performance is maximized when
the target is characterized by an optimal, intermediate frequency which
monotonically decreases with the strength of the internal reservoir
connectivity. At the optimal frequency, the reservoir representation of the
target signal is high-dimensional, de-synchronized, and thus maximally robust
to noise. We show that our predictions successfully capture the qualitative
behaviour of performance in non-linear networks. Moreover, we find that the
relationship between internal representations and performance can be further
exploited in trained non-linear networks to explain behaviours which do not
have a linear counterpart. Our results indicate that a major determinant of
learning success is the quality of the internal representation of the target,
which in turn is shaped by an interplay between parameters controlling the
internal network and those defining the task.

##  1 Introduction

A fundamental feature of the brain, and biological networks in general, is the
ability to form closed-loop interactions with their environment. Such
interactions are often implemented through a dimensionality bottleneck: while
networks typically consist of large numbers of units, signals exchanged with
the environment are low-dimensional. In fact, external stimuli can often be
represented in terms of a few scalar variables (e.g. the angle and speed of a
tennis ball approaching); these low-dimensional variables are encoded in the
high-dimensional activity of a large population of neurons [1, 2] before being
again transformed into low-dimensional decision variables and motor outputs
(e.g. the angle and speed of the hand holding the racket).

Simple but effective models for studying closed-loop interactions are feedback
networks. These models implement a simple form of closed-loop interaction: the
output (or readout) signal, which is extracted from a reservoir of randomly
connected units as a linear combination of unit activities, is directly
injected back into the reservoir as external input [3, 4]. By adjusting the
weights which specify how reservoir activity is mapped to the output, feedback
networks can be trained to produce the desired readout signal. In the most
common training algorithms [5, 6, 7], readout weights are updated through
least-squares (LS) regression; this can be performed only once, by using a
complete batch of activity samples [5], or in an online fashion, by
recursively integrating activity samples as they are simulated [7, 8].

What kind of closed-loop dynamics can feedback networks implement? Despite
some theoretical advancement [4, 9, 10, 11, 12], computational properties of
feedback networks are still poorly understood. Early theoretical work has
indicated that most feedback models are expected to be able to approximate
readout signals characterized by arbitrarily complex dynamics [4]. However, it
has been reported that not all feedback architectures and target dynamics
result in the same performance: trained networks can experience dynamical
instabilities [10, 11], and converge to fragile solutions for certain choices
of the feedback architecture and parameters [7, 13].

For a fixed task, several studies have reported that training performance is
strongly influenced by the overall strength of recurrent connections in the
reservoir [14, 7, 6]. Specifically, performance is high when recurrent
connections are strong, but not strong enough to lead to the appearance of
chaotic activity [15] – a parameter region named _edge-of-chaos_ [16].
Intuitively, the edge-of-chaos defines an optimal tradeoff point where the
internal reservoir dynamics are rich but stable.

Reservoir activity, however, is not determined by connectivity alone: because
the system is coupled to the environment, activity depends also on the
statistics and dynamics of the target signal, which specify the task. How the
internal reservoir dynamics interact with the target in determining trained
networks performance is a fundamental question in feedback systems which is
still not well understood [17]. In particular: are there specific target
features which optimize performance, and how do they depend on internal
properties of the reservoir network? For given values of the target
parameters, what are the properties of reservoir activity that support optimal
training? How sensitive is the optimal performance to the learning algorithm?
To the current date, these questions remain largely unsolved.

In this work, we consider a simple setup consisting of a non-linear reservoir
of rate units which is trained to sustain a sinusoidal output with given
frequency $\omega$. Consistently across three different training techniques,
we find that learning performance is maximized at a finite “preferred”
frequency $\bar{\omega}$, which in turn depends on reservoir connectivity: as
the connectivity strength is increased towards the edge of chaos,
$\bar{\omega}$ decreases towards zero. This nontrivial dependence of
performance, even in a simple task, provides a test case to study the
interplay between reservoir and target properties and its effects on learning.

To gain analytical insight into this phenomenon, we consider a simplified
setup where reservoir dynamics are linearized, and perform exact mathematical
analysis. By averaging over the ensemble of random reservoir networks, we
characterize reservoir activity in response to the target signal, and show
that a “resonance” frequency $\omega^{*}$ emerges, which decreases with the
connectivity strength. Under this frequency, dimensionality of neural activity
is maximal and synchrony across different units in the reservoir is minimal.
When training the network to output the target signal, feedback interactions
are most robust in the vicinity of the resonance frequency, thus resulting in
optimal performance. Moreover, this behaviour is predicted to be qualitatively
consistent across different training algorithms, even if performance itself is
sensitive to the algorithm used. We show that our theoretical predictions
correctly capture the qualitative behaviour of learning performance observed
numerically in non-linear network models. Overall, our results shed light on
the learning capacity of recurrent network architectures by quantifying how
learning precision is determined by the interaction between internal reservoir
connectivity and target dynamics.

##  2 Results

###  2.1 Emergence of a preferred frequency in trained feedback networks

We consider a reservoir network consisting of $N$ units characterized by the
evolution dynamics:

$\dot{\mathbf{x}}(t)=-\mathbf{x}(t)+\mathbf{J}\Phi(\mathbf{x}(t))+\mathbf{m}u(t)$
(1)

where $\Phi(x)=\tanh(x)$ is applied to the activation vector $\mathbf{x}$
element-wise. Recurrent weights $\mathbf{J}$ are fixed, and are generated
independently from the normalized Gaussian distribution
$\mathcal{N}(0,g^{2}/N)$ [15, 7], so that the parameter $g$ controls the
strength of reservoir connectivity. The one-dimensional external signal $u(t)$
acts as a forcing on the reservoir through input weights $\mathbf{m}$, which
are fixed and drawn as independent standard Gaussian variables.

The output of the reservoir network is a one-dimensional readout signal,
defined as:

$z(t)=\mathbf{n}^{\top}\Phi(\mathbf{x}(t))$ (2)

through a set of decoding weights $\mathbf{n}$ that are assumed to be plastic.
The feedback is realized by using the output signal as input: $u(t)=z(t)$ (see
Fig. 1A for an illustration), which yields the final autonomous dynamics

$\dot{\mathbf{x}}(t)=-\mathbf{x}(t)+(\mathbf{J}+\mathbf{m}\mathbf{n}^{\top})\Phi(\mathbf{x}(t)).$
(3)

During training, the vector $\mathbf{n}$ is updated until the output $z(t)$
best matches the desired target $f(t)$. The target function that we consider
is a simple sinusoidal wave of frequency $\omega$, i.e. $f(t)=A\cos(\omega
t)$.

We trained multiple instances of this feedback architecture and analyzed how
performance depends on the frequency of the target signal $\omega$ and on
internal coupling strength $g$ (Fig. 1). Three common training algorithms
(least-squares (LS) regression, ridge regression [18] and recursive least-
squares (RLS) [19, 7]) were used (training details are reported in Appendix
4.1). We quantified the error as the mismatch between the target $f(t)$ and
the readout $z(t)$ averaged over a finite number of target cycles in the post-
training activity.

Figure 1:  Emergence of preferred frequency in non-linear feedback networks
trained to sustain a sinusoidal output. A. Illustration of network
architecture used in open (yellow) and closed (purple) loop. B. Example
readout signal (dark grey; target shown in yellow) for three learning trials
corresponding to the frequency values indicated in green in C ($w=$ 0.1, 0.7,
2.1; $g=1$). Other parameters as in C, except (for illustration purposes)
training is performed on a smaller number of target cycles ($N^{\text{tot}}=4$
and $N^{\text{tr}}=2$, see Appendix 4.1). LS regression was used for training;
examples trials for Ridge and RLS are reported in Supp. Fig. 8. C. Readout
error as a function of $\omega$, for a range of $g$ values (blue shades), for
networks trained via LS (left), ridge regression (middle) and RLS (right). We
take here $A=1$. Training details and parameters are reported in Appendix 4.1.
D. Error-minimizing frequency $\bar{\omega}$ as a function of $g$, for the
three learning algorithms as in C. Three different target amplitudes $A$ were
tested (grey shades).

We observe that, for fixed reservoir connectivity $g$, the accuracy of signal
reconstruction by the output strongly depends on the target frequency: while
training on one frequency results in highly precise readout for many cycles,
others result in a runaway from the target signal (Fig. 1B). For every value
of $g$, the error has a non-monotonous dependence on $\omega$, and reaches a
minimum at a finite frequency that we name $\bar{\omega}$ (Fig. 1C). Each
curve, corresponding to a different value of $g$, has a different optimal
frequency: specifically, $\bar{\omega}$ decreases as the strength of reservoir
connectivity $g$ increases from zero towards the edge-of-chaos (Fig. 1D; see
Appendix 4.11 for a characterization of the edge-of-chaos in our framework).
Although the exact value of the preferred frequency $\bar{\omega}$ is found to
be algorithm-dependent, the same qualitative behaviour is observed
consistently across the three different algorithms we used for training. It is
also observed for both small and large amplitudes of the target signal $A$,
which are expected to elicit, respectively, weakly or strongly non-linear
activity in the reservoir.

The observations from Fig. 1 provide a striking example of the non-trivial
interplay between reservoir features (the connectivity parameter $g$) and
external task parameters (the target frequency $\omega$) in determining
learning performance. Because the network is completely random, one might
naively think that its dynamics do not exhibit a typical timescale, and are
thus blind to the signal frequency; instead, the network appears to have its
preference even for a simple task. In the rest of this paper, we aim to
understand this observation in detail through mathematical analysis.

To this end, we consider a simplified model which greatly eases the analysis:
the case of linear reservoir dynamics ($\Phi(x)=x$). The analysis strategy we
use consists of two steps [10]. To begin with, we examine the feedback network
in an open-loop setup (Fig. 1A, yellow), where the encoding of the input and
the decoding of the output signals can be analyzed separately. In the encoding
phase, we take the input to the reservoir network to be identical to the
target function: $u(t)=f(t)$, and characterize analytically the reservoir
response $\mathbf{x}(t)$ both at the level of single units and the population
as a whole (Section 2.2). In the decoding phase, we use the reservoir response
to pick a readout $\mathbf{n}$ which allows the network to reconstruct the
correct output: $\mathbf{n}^{\top}\mathbf{x}(t)=f(t)$ (Section 2.3). At that
point, our feedback architecture admits the desired target as a solution; to
investigate success of such solutions in performing the task, in Section 2.4
we close the loop (Fig. 1A, purple), and characterize dynamics stability.
Taken together, the open- and closed-loop descriptions fully characterize
trained feedback architectures, and thus allow us to make predictions about
dynamical mechanisms and training performance. These predictions are shown to
compare favorably to numerical simulations obtained by training linear
feedback networks (Section 2.5). Finally, in Section 2.6, we show that they
qualitatively carry over to the case of non-linear dynamics.

###  2.2 Open loop: encoding the target signal

We begin our analysis by examining encoding in the open-loop framework: this
corresponds to a random reservoir with linear dynamics driven by the target
signal. The time evolution is described by

$\dot{\mathbf{x}}(t)=-\mathbf{x}(t)+\mathbf{J}\mathbf{x}(t)+\mathbf{m}f(t);$
(4)

here $\mathbf{J}$ is a Gaussian random matrix as defined above; to avoid
dynamic instabilities, we consider $g<1$ [20]. The linear dynamics is
indifferent to the amplitude of the input, so we set $A=1$; in response to the
periodic input $f(t)=\text{cos}(\omega t)$, the stationary solution for
$t\to\infty$ is

$\mathbf{x}(t)=\frac{1}{2}\left(\mathbf{x}_{+}e^{i\omega
t}+\mathbf{x}_{-}e^{-i\omega t}\right),$ (5)

where

$\mathbf{x}_{\pm}:=[(1\pm i\omega)\mathbf{I}-\mathbf{J}]^{-1}\mathbf{m}$ (6)

are complex conjugate vectors representing the reservoir activity in Fourier
space (see Appendix 4.2).

#### A geometric description

The stationary solution may be written as

$\mathbf{x}(t)=\mathbf{v}_{+}\text{cos}(\omega
t)+\mathbf{v}_{-}\text{sin}(\omega t),$ (7)

showing that activity occupies the plane spanned by two vectors
$\mathbf{v}_{\pm}$, which are the real and imaginary parts of
$\mathbf{x}_{\pm}$: $\mathcal{R}(\mathbf{x}_{+})=\mathbf{v}_{+}$ and
$\mathcal{I}(\mathbf{x}_{+})=-\mathbf{v}_{-}$. In this plane, the state-space
trajectory is a closed, elliptic curve (Fig. 2A), with geometry determined by
the spanning vectors.

Figure 2:  Encoding of the target signal: geometry of activity trajectories.
A. Projection of one example trajectory $\mathbf{x}(t)$ (Eq. (7)) on the plane
spanned by vectors $\mathbf{v}_{\pm}$. B. Norm of the two vectors
$\mathbf{v}_{\pm}$ as a function of the target frequency $\omega$. C. Angle
between the two spanning vectors $\mathbf{v}_{\pm}$. In A-B-C we used $g=0.5$.
In B-C, the black vertical line indicates the resonance frequency $\omega^{*}$
where the two norms are equal ($r=1$, B) and the angle is maximized (C). D.
Linear dimensionality: participation ratio computed from the principal
components of activity. We plot results for five increasing values of $g$
(blue shades); black stars indicate the position of $\omega^{*}$ for every
value of $g$. E. Resonance frequency $\omega^{*}$. In all panels, continuous
lines indicate the analytical predictions. In D, dots show results averaged
over $20$ simulations of finite size networks, $N=2000$.

The spanning vectors $\mathbf{v}_{\pm}$, in turn, depend both on the recurrent
connectivity $\mathbf{J}$ and on the driving frequency $\omega$ (Eq. (6)).
Their geometry is self-averaging in the limit of large networks, and can be
computed by averaging over the ensemble of randomly connected reservoir
networks (see Appendix 4.3). Fig. 2B shows the dependence upon $\omega$ of the
norms $\|\mathbf{v}_{\pm}\|$. For very small frequencies, the trajectory
follows the drive adiabatically and $\mathbf{v_{-}}\approx 0$; there is
practically only one spanning vector. As frequency increases, the response
acquires a phase shift and the second spanning vector $\mathbf{v_{-}}$ becomes
non-negligible. At high frequencies, both norms decrease due to the filtering
property of the network; the second spanning vector thus obtains a maximal
norm at an intermediate frequency.

We quantify the elliptical trajectory by its linear dimensionality, i.e. the
participation ratio computed from the principal components of reservoir
activity [21, 22]. Denoting the activity cross-correlation matrix by
$C:=\frac{1}{T}\int_{0}^{T}\mathbf{x}(t)\mathbf{x}^{\top}\>\mathop{}\\!\mathrm{d}t$
and its eigenvalues by $\nu_{i}$, the trajectory dimensionality $d$ is defined
as

$d:=\frac{\left(\sum_{i=1}^{N}\nu_{i}\right)^{2}}{\sum_{i=1}^{N}\nu_{i}^{2}}.$
(8)

By using Eq. (7), and by integrating out time, we find that $C$ is a rank-two
matrix,
$C=\frac{1}{2}\left(\mathbf{v}_{+}\mathbf{v}_{+}^{\top}+\mathbf{v}_{-}\mathbf{v}_{-}^{\top}\right)$,
whose non-zero eigenvalues (which we take to be $\nu_{1},\nu_{2}$) are
identical to those of the $2\times 2$ reduced cross-correlation matrix [23]

${C}^{R}=\begin{pmatrix}\|\mathbf{v}_{+}\|^{2}&\mathbf{v}_{+}\cdot\mathbf{v}_{-}\\\
\mathbf{v}_{+}\cdot\mathbf{v}_{-}&\|\mathbf{v}_{-}\|^{2}\end{pmatrix}.$ (9)

Explicitly computing the eigenvalues of ${C}^{R}$ yields the expression

$d=\frac{\left(\|\mathbf{v}_{+}\|^{2}+\|\mathbf{v}_{-}\|^{2}\right)^{2}}{\|\mathbf{v}_{+}\|^{4}+2\left(\mathbf{v}_{+}\cdot\mathbf{v}_{-}\right)^{2}+\|\mathbf{v}_{-}\|^{4}}.$
(10)

We observe that the linear dimensionality, which is bounded between 1 and 2,
is insensitive to the overall trajectory magnitude, but depends on the ratio
of norms $r=\|\mathbf{v}_{-}\|/\|\mathbf{v}_{+}\|$ and on the angle $\theta$
between the spanning vectors:

$d=\frac{1}{1-\frac{2r^{2}}{(1+r^{2})^{2}}\,\sin^{2}(\theta)}.$ (11)

The ratio $r$ indicates how much the curve is squeezed along a single
direction, with both extremes ($r$ very small or very large) resulting in
trajectories squeezed along the dominant spanning vector. For a fixed angle,
as the ratio passes through $r=1$, the trajectory goes through a shape which
is most similar to a circle and has maximal dimensionality. For a fixed $r$,
the angle $\theta$ determines to what degree the curve is skewed relative to a
perfect ellipse; the dimensionality increases monotonically as $\theta$ opens
up from zero to $\pi/2$.

Examination of the vector norms $\|\mathbf{v}_{\pm}\|$ in Fig. 2B indicates
that they intersect at a frequency value that we name $\omega^{*}$, where
$r=1$. Fig. 2C shows how the angle $\theta$ varies as a function of frequency;
surprisingly, we find that it displays a maximum at $\omega^{*}$. These
dependencies are reflected in the behaviour of the dimensionality (Fig. 2D),
which itself attains a maximum at frequency $\omega^{*}$. Our mathematical
analysis reveals that (see Appendix 4.4)

$\omega^{*}=\sqrt{1-g^{2}},$ (12)

i.e. the resonance frequency $\omega^{*}$ decreases to zero as $g$ increases
towards the instability boundary ($g=1$). This analytic result is in excellent
agreement with finite network simulations, as shown in Fig. 2E.

#### A single-unit description

The analysis above considered the geometry of trajectories spanned by the
reservoir population in its high-dimensional activity space, and revealed that
trajectory dimensionality is maximized at the resonance frequency
$\omega^{*}$. An alternative viewpoint is obtained by considering the
statistics of single-unit activity profiles across the population. As we shall
see, this alternative perspective reveals that the optimal frequency
$\omega^{*}$ has a second natural interpretation in terms of population
synchrony.

To do so, we derive a self-consistent expression for $\mathbf{x}_{+}$ by
inserting Eq. (5) into the evolution equations (Eq. (4)):

$\mathbf{x}_{+}=\frac{1}{1+i\omega}\left({\mathbf{m}}+\mathbf{J}\mathbf{x}_{+}\right).$
(13)

This form highlights that vector $\mathbf{x}_{+}$ is given by the sum of two
contributions: one associated with the external forcing via the input vector
$\mathbf{m}$, and one associated with the reservoir response via the recurrent
input $\mathbf{Jx}_{+}$. Since $\mathbf{J}$ is random, the direction of the
latter contribution is random (i.e., it varies across realizations of
$\mathbf{J}$), but its amplitude is self-averaging and depends on the strength
of recurrent connectivity $g$ [15].

We use Eq. (13) to gain intuition about how the network encodes the external
oscillatory signal at the level of single-unit activity. To this end, we
visualize the entries of the $\mathbf{x}_{+}$ vector as points in the complex
plane: $(\mathbf{x}_{+})_{i}={R_{i}}e^{i{\phi_{i}}}$, where $R_{i}$ and
$\phi_{i}$ represent the amplitude and phase with which a single unit responds
to the forcing input (Fig. 3A). How are points corresponding to different
units distributed on the complex plane? When recurrent connections are very
weak ($g\simeq 0$), different units behave as uncoupled filters of the input;
we have $\mathbf{x}_{+}\simeq{\mathbf{m}}/(1+i\omega)$, implying that the real
and imaginary part of $(\mathbf{x}_{+})_{i}$ for different $i$ are
proportional one to each other. As a consequence, points on the complex plane
are collinear (Fig. 3A left), and phases are identical: $\phi_{i}=\phi$.
Responses of different units are thus synchronized (Fig. 3B left). As $g$
grows from 0, the second term in Eq. (13), which originates from recurrent
interactions, starts spreading the real and imaginary parts of
$(\mathbf{x}_{+})_{i}$ away from the line $\phi_{i}=\phi$ (Fig. 3A right), and
introduces variability in response phases (Fig. 3B right).

Figure 3:  Encoding of the target signal: single-unit description. A. Entries
of vector $\mathbf{x}_{+}$ in the complex plane. Left: $g=0.1$, right:
$g=0.5$; $w=0.6$ for both panels. Continuous lines are contour lines of the
bivariate Gaussian distribution predicted by the theory. The most external
contour indicates a probability of 0.01. Grey points are results of a finite
network simulation ($N=400$). B. Sample of activity from four randomly
selected units chosen from the corresponding panel in A. C. Spread of response
phases across the population (Eq. (15)) for increasing values of $g$ and as a
function of $\omega$. Black stars indicate the maximum value. D. Value of the
frequency which maximizes the spread of phases from C.

For fixed values of $g$ and $\omega$, the distribution of dots on the complex
plane is a bivariate Gaussian (Fig. 3A); a _narrow_ distribution corresponds
to highly synchronized units, and its _broadening_ at stronger coupling
indicates their desynchronization. As both $\mathbf{m}$ and $\mathbf{J}$ are
generated from a centered Gaussian distribution, the mean of the distribution
vanishes. The covariance is given by:

$\frac{1}{N}\begin{pmatrix}\|\mathbf{v}_{+}\|^{2}&-\mathbf{v}_{+}\cdot\mathbf{v}_{-}\\\
-\mathbf{v}_{+}\cdot\mathbf{v}_{-}&\|\mathbf{v}_{-}\|^{2}\end{pmatrix}$ (14)

implying that the shape distribution is controlled by the statistics of the
spanning vectors $\mathbf{v}_{+}$ and $\mathbf{v}_{-}$. The similarity between
the covariance matrix and the reduced cross-correlation matrix $C^{R}$ (Eq.
(9)) analyzed in the previous paragraph suggests that synchrony in single-unit
response and dimensionality of state-space trajectories are deeply related
properties of reservoir activity. To formalize this relationship, we compute
the spread of phases $\phi_{i}$ across the reservoir population, i.e.

$\Sigma^{2}=\int_{\bar{\phi}-\frac{\pi}{2}}^{\bar{\phi}+\frac{\pi}{2}}\mathop{}\\!\mathrm{d}\phi\>p(\phi)(\phi-\bar{\phi})^{2},$
(15)

where $p(\phi)$ is the probability distribution of phases for a bivariate
Gaussian distribution ([24], see Appendix 4.5). The phase spread for different
values of recurrent strength $g$ and frequency $\omega$ is plotted in Fig. 3C.
These results show that it monotonically increases with $g$; for any fixed
$g$, it reaches a maximum at a finite frequency value, given again by
$\omega^{*}=\sqrt{1-g^{2}}$ (Fig. 3D).

To conclude, we have examined the behaviour of single-unit activity in
response to a sinusoidal forcing input. In line with classical mean-field
studies, we have analyzed the statistical distribution of single-unit activity
profiles across the reservoir population [15, 25, 9]. This approach has
revealed that, for fixed $g$, $\omega^{*}$ corresponds to the frequency at
which single-unit activity is maximally desynchronized. Note that
historically, desynchronized single-unit profiles have been pointed out as a
desirable feature of reservoir activity, as temporally heterogeneous profiles
form a rich set of basis functions from which complex target functions can be
reconstructed [26].

###  2.3 Open-loop setup: decoding the internal representation

After having characterized the reservoir activity during stimulus encoding, we
turn to the decoding step of the open-loop analysis. Decoding corresponds to
finding a readout vector $\mathbf{n}\in\mathbb{R}^{N}$ which satisfies:

$z(t)=\mathbf{n}^{\top}\mathbf{x}(t)=\cos(\omega t);$ (16)

the projection of driven reservoir activity along $\mathbf{n}$ needs thus to
match the target $f(t)$ (Fig. 4A, yellow). In terms of the Fourier-space
representation (Eq. (5)), $\mathbf{n}$ is a solution to the set of two linear
equations given by

$\mathbf{n}^{\top}\mathbf{x}_{\pm}=\mathbf{n}^{\top}[(1\pm
i\omega)\mathbf{I}-\mathbf{J}]^{-1}\mathbf{m}=1.$ (17)

When $g=0$, interactions vanish and the equations above read
$\mathbf{n}^{\top}\mathbf{m}=1\pm i\omega$, which cannot be satisfied by any
$\mathbf{n}$. This scenario corresponds to completely synchronized reservoir
activity, or equivalently, activity spanning one-dimensional state-space
trajectories. For any $g>0$, on the other hand, this system of equations is
under-determined, since it fixes only 2 among the $N$ degrees of freedom in
$\mathbf{n}$.

We explore the effect of these degrees of freedom by defining a family of
readout vectors $\mathbf{n}$ parametrized by an integer $k$, where
$k=2,\ldots,N$; $k$ indicates the number of reservoir units from which the
readout signal is reconstructed. We term such solutions _from-_ $k$
regression. To obtain such a solution, we set all elements except for the
first $k$ of $\mathbf{n}$ to zero, and then solve Eq. (17) by considering the
least-squares (LS) solution of minimal norm, which can be computed through the
pseudo-inverse (see Appendix 4.6). When $k=N$, we obtain the full LS solution,
which reads:

$\mathbf{n}_{LS}=\begin{pmatrix}\mathbf{x}_{+}&\mathbf{x}_{-}\end{pmatrix}{\begin{pmatrix}\|\mathbf{x}_{+}\|^{2}&\mathbf{x}_{+}\cdot\mathbf{x}_{-}\\\
\mathbf{x}_{+}\cdot\mathbf{x}_{-}&\|\mathbf{x}_{-}\|^{2}\end{pmatrix}}^{-1}\begin{pmatrix}1\\\
1\end{pmatrix}$ (18)

or, in terms of $\mathbf{v_{\pm}}$ vectors:

$\displaystyle\mathbf{n}_{LS}=\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}P\begin{pmatrix}1\\\
0\end{pmatrix},$ (19)

where we defined the short-hand notation $P=\left(C^{R}\right)^{-1}$.

All the readouts within the _from-k_ family exactly solve the task in the
open-loop setup. However, it is not clear a-priori whether all of them are
equivalent when closing the loop, i.e. when the feedback network is required
to autonomously generate the target signal (Eq. (3)). In the following, we
assess dynamics and stability of closed-loop networks corresponding to the
different choices of the readout $\mathbf{n}$.

###  2.4 Closing the loop: autonomous signal generation

In the previous two sections, we have analyzed how random networks encode a
one-dimensional periodic signal, and how the network response can be used to
reconstruct the same signal as output. Ultimately, we want the encoding and
the decoding steps to be self-consistent, i.e. we require

$z(t)=\mathbf{n}^{\top}\mathbf{x}(t)=u(t)=f(t),$ (20)

which is equivalent to transforming our problem from an open-loop to a closed-
loop setup, where dynamics are autonomous and follow Eq. (3) with linear
interactions:

$\dot{\mathbf{x}}(t)=-\mathbf{x}(t)+(\mathbf{J}+\mathbf{m}\mathbf{n}^{\top})\mathbf{x}(t)$
(21)

and $\mathbf{n}$ satisfies Eq. (16). This step is illustrated in Fig. 4A by
the purple feedback arrow connecting the reservoir output to the input. If
closing the loop does not perturb activity in the reservoir by changing its
stability properties, then at every time point the readout
$\mathbf{n}^{\top}\mathbf{x}(t)=\cos(\omega t)$ is fed back into the system,
and the solution obtained through the open-loop setup is self-consistent.

The solutions to Eq. (21) and their stability are fully characterized by the
eigenspectrum of $\mathbf{\bar{J}}=\mathbf{J}+\mathbf{m}\mathbf{n}^{\top}$
(the leak term in the dynamics contributes by uniformly shifting the spectrum
by $-1$). For $N$ sufficiently large, the eigenvalues of $\mathbf{J}$ are
distributed uniformly in a disk of radius $g<1$ [20]. The position of some or
all of the eigenvalues can, however, be modified by the rank-one perturbation
$\mathbf{mn}^{\top}$; we refer to these as _outliers_. In order for the
closed-loop system to stably sustain the periodic activity we found in the
encoding step, the eigenspectrum of $\mathbf{\bar{J}}$ must satisfy two key
requirements: (i) a pair of complex outlier eigenvalues with value:
$\lambda_{\pm}=1\pm i\omega$ (which ensures that a periodic trajectory of
frequency $\omega$ is realized), and (ii) a stable bulk of remaining
eigenvalues: $\mathcal{R}(\lambda)<1$ $\forall\lambda\neq\lambda_{\pm}$ (which
ensures that no runaway activity is generated along other directions).

Figure 4:  Closing the loop. A. Transforming the open-loop encoding/decoding
setup (yellow) into a closed-loop system (purple). B. Sample networks trained
through full LS (_from-N_ , left) or _from-2_ (right) regression. The top
panels show the eigenspectra of the closed-loop connectivity matrix
$\mathbf{\bar{J}}$ (red dots); small black dots indicate the unperturbed
eigenspectrum of $\mathbf{{J}}$. The bottom panels show the output generated
by the corresponding networks. Here we used parameters $g=0.8$, $\omega=0.6$
and $N=400$. C. Overlap between readout $\mathbf{n}$ and the principal
components (PC) of driven reservoir activity (Eq. (7)), of which the first two
span the $\mathbf{v}_{\pm}$ plane. The same parameters as in B were used. D.
Fraction of unstable closed-loop systems (over 2000 sample networks) as a
function of connectivity strength $g$, for several values of $k$, measured
over 1000 different realizations. We used $N=1000$ and $\omega=0.6$.

All eigenvalues of $\mathbf{\bar{J}}$ are roots of the characteristic
polynomial

$\det\left((\mathbf{J}+\mathbf{m}\mathbf{n}^{\top})-\lambda\mathbf{I}\right)=0.$
(22)

The Matrix Determinant Lemma [13, 12] allows us to decompose this polynomial
into two factors, corresponding to the two sets of eigenvalues:

$\det\left(1+\mathbf{n}^{\top}(\mathbf{J}-\lambda\mathbf{I})^{-1}\mathbf{m}\right)\det\left(\mathbf{J}-\lambda\mathbf{I}\right)=0.$
(23)

It is seen that the second term vanishes on the spectrum of $\mathbf{J}$,
whereas the first term vanishes for the outlier eigenvalues. The outliers
therefore satisfy

$1=\mathbf{n}^{\top}(\lambda\mathbf{I}-\mathbf{J})^{-1}\mathbf{m}.$ (24)

If $\mathbf{n}$ satisfies Eq. (16), then $\lambda_{\pm}=1\pm i\omega$ are
indeed solutions, implying that condition (i) is satisfied. Note that the
eigenvectors corresponding to $\lambda_{\pm}$ are identical to the vectors
$\textbf{x}_{\pm}$, as
$(\mathbf{J}+\mathbf{m}\mathbf{n}^{\top})\mathbf{u}_{\pm}=\lambda_{\pm}\mathbf{u}_{\pm}$
implies

$\mathbf{u}_{\pm}\propto(\mathbf{J}-\lambda_{\pm}\mathbf{I})^{-1}\mathbf{m}.$
(25)

Thus fixing $\mathbf{n}$ in the open-loop framework is equivalent to directly
controlling the value of the target-relevant eigenvalues in the eigenspectrum
of the closed-loop network.

We next examine whether this pair of eigenvalues are the only outliers
generated by closing the loop: while Eq. (24) is guaranteed to have
$\lambda_{\pm}=1\pm i\omega$ as solutions, other solutions might be admitted
which could violate requirement (ii). Such potential solutions depend on the
overlap between the vectors $\mathbf{n}$ and
$\mathbf{x}_{\lambda}=(\lambda\mathbf{I}-\mathbf{J})^{-1}\mathbf{m}$. Note
that if the readout $\mathbf{n}$ was random (and thus orthogonal to
$\mathbf{J}$ and $\mathbf{m}$), this overlap would vanish and no additional
outliers would be generated.

In the case of the full LS solution ($k=N$, Eqs. (18)-(19)), the readout
vector $\mathbf{n}_{LS}$ is contained in the plane spanned by vectors
$\mathbf{v}_{\pm}$. As a consequence, the overlap between $\mathbf{n}$ and
$\mathbf{x}_{\lambda}$ can be expanded in terms of
$\mathbf{x}_{\pm}^{\top}\mathbf{x}_{\lambda}$. In the limit
$N\rightarrow\infty$, these terms have a simple form which can be evaluated
analytically (see Appendix 4.8), yielding an equation in $\lambda$ which
reads:

$(1+\omega^{2}){\lambda}^{2}-\left[2g^{2}+N(P_{11}+\omega
P_{21})\right]{\lambda}+g^{4}+Ng^{2}P_{11}=0$ (26)

where $P_{11}$ and $P_{21}$ are the elements of the first column of
$P=\left(C^{R}\right)^{-1}$ and depend on $g$ and $\omega$. As the equation
above is quadratic, it admits $\lambda=\lambda_{\pm}$ as unique solutions.
Therefore, in large networks, LS training is guaranteed to result in stable
dynamics, as no additional outliers are generated in the eigenspectrum other
than the task-relevant ones. This is confirmed by numerical simulation in the
left panels of Fig. 4B.

In the more general case of _from-k_ LS regressors with $k<N$, readout vectors
might contain extra components that are correlated with $\mathbf{J}$ and
$\mathbf{m}$ and are not fully contained within the spanning plane; as a
consequence, more than two outlier eigenvalues and unstable dynamics can be
expected. The right panels of Fig. 4B show an example of such a situation (in
the simulation, $k=2$). One outlier eigenvalue with $\mathcal{R}(\lambda)>1$
is seen in the top panel, which induces the dynamic instability seen in the
bottom panel.

As $k$ decreases from the maximal value ($N$) to the minimal one ($2$), the
component of the readout vector $\mathbf{n}$ outside of the $\mathbf{v}_{\pm}$
plane becomes larger (Fig. 4C). Numerical analysis indicates that,
correspondingly, the fraction of networks with unstable dynamics increases
(Fig. 4D).

In summary, we have shown that – although the open-loop setup admits multiple
exact solutions – different solutions are not equivalent in terms of dynamical
stability when the loop is closed. Stability properties are related to the
orientation of the readout vector relative to the driven open-loop trajectory.
In the case of the full LS solution (Eq. (18)), the readout $\mathbf{n}_{LS}$
is completely aligned with the trajectory plane, and closed-loop dynamics are
guaranteed to be stable. Other solutions generally contain components outside
of this plane, which can cause activity to diverge.

###  2.5 Predicting performance of trained linear networks

We now turn back to the problem of understanding performance in trained
feedback networks and its dependence on the target frequency $\omega$. We
start by considering linear feedback networks which are trained (as in Fig. 1C
left) via LS regression. The analysis of the previous sections has shown that
linear networks trained via LS regression can exactly implement the feedback
task with stable dynamics. Due to noise, however, real-life LS regression
never converges to this ideal solution. Noise arises in training from multiple
sources, such as finite sampling of training data, variability due to
different initial conditions or regularization noise. In order to characterize
training performance, we thus use our theoretical framework to analyze the
effect of noise on the dynamics of feedback linear networks trained via LS
regression.

Consider first the encoding phase of learning (Section 2.2), where reservoir
activity is stimulated. Because of noise, learning algorithms may not have
access to the true spanning vectors $\mathbf{v}_{\pm}$. Rather, we assume that
corrupted versions $\tilde{\mathbf{v}}_{\pm}=\mathbf{v}_{\pm}+\bm{\xi}_{\pm}$
(where the entries of $\bm{\xi}_{\pm}$ are independent Gaussian noise) are
measured. The estimated LS readout then reads:

$\tilde{\mathbf{n}}_{LS}=\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}\tilde{P}\begin{pmatrix}1\\\
0\end{pmatrix}\>+\>\begin{pmatrix}\bm{\xi}_{+}&\bm{\xi}_{-}\end{pmatrix}\tilde{P}\begin{pmatrix}1\\\
0\end{pmatrix}$ (27)

where $\tilde{P}$ is the inverse reduced cross-correlation matrix which
includes the noise disturbance.

As in Section 2.4, we can characterize closed-loop dynamics by computing the
outlier eigenvalues of
$\mathbf{J}+\mathbf{m}\>{\tilde{\mathbf{n}}_{LS}}^{\top}$. The second term in
the r.h.s. of Eq. (27) is random and orthogonal to $\mathbf{m}$ and
$\mathbf{J}$, and therefore does not affect the position of outlier
eigenvalues. In contrast, the first term is a vector fully aligned with the
noise-free spanning vectors $\mathbf{v}_{\pm}$, which generates two outlier
eigenvalues $\tilde{\lambda}_{\pm}$. Because of the noise, their values
deviate from the target eigenvalues $\lambda_{\pm}$; they are solutions of an
equation identical to Eq. (26), but with $P$ replaced by $\tilde{P}$. For
every noise realization, the inverse reduced cross-correlation matrix $P$ is
perturbed in a random direction, yielding random modifications to the target
eigenvalues ${\lambda}_{\pm}$. We can estimate the average mismatch between
$\tilde{\lambda}_{\pm}$ and ${\lambda}_{\pm}$ from the sensitivity of matrix
$P$ to perturbations, which is quantified by the condition number of the
reduced correlation matrix $C^{R}$, i.e. the ratio between the largest and
smallest eigenvalue [27, 17]:

$c=\frac{\nu_{2}}{\nu_{1}}.$ (28)

The value of $c$ and its dependence on $\omega$ and $g$ can be computed by
taking the limit $N\rightarrow\infty$ and averaging over the networks
ensemble. Fig. 5A shows that, for fixed connectivity strength $g$, the
condition number is a non-monotonic function of the forcing frequency
$\omega$, and attains a minimum at the resonance frequency
$\omega^{*}=\sqrt{1-g^{2}}$ (see Appendix 4.9). Thus, when training linear
feedback networks through noisy LS regression we expect that, for fixed $g$,
the readout would be closest to the desired one at $\bar{\omega}=\omega^{*}$,
where the reduced cross-correlation matrix is most robust to noise. This
robustness directly reflects the properties of the internal representation of
the target signal within the reservoir, which is characterized by maximal
dimensionality and minimal synchrony at $\omega=\omega^{*}$.

We tested this prediction on finite-size trained networks. Examples from Fig.
5B confirm that the task-related eigenvalue pair $\tilde{\lambda}_{\pm}$
deviate from the target ones. As in the case of non-linear networks (Fig.
1B-C), we find that the error is frequency dependent (Fig. 5C, left);
furthermore, for fixed strength of the internal connectivity $g$, we observe
that the error is minimized at a frequency $\bar{\omega}$ which is very close
to $\omega^{*}$ (Fig. 5D, left).

Figure 5:  Training performance in linear networks is maximized at
$\omega^{*}$. A. Condition number of the reduced cross-correlation matrix
$C^{R}$ computed analytically for different values of $g$ (blue shades in C).
Stars denote minimal condition number. B. Closed-loop outlier eigenvalues
$\tilde{\lambda}_{\pm}$ for example networks from three learning trials
corresponding to the frequencies marked by green triangles in C ($w=$ 0.1,
0.4, 1.5; $g=0.9$). C. Spectrum error as a function of $\omega$, for a range
of $g$ values (blue shades), for networks trained via noisy LS (left), ridge
regression (middle) and RLS (right). Error is measured from the imaginary part
of outlier eigenvalues as
$|\mathcal{I}(\tilde{\lambda}_{\pm})-\mathcal{I}(\lambda_{\pm})|/\omega$, and
similarly for the real part; the spectrum error is an average of the two.
Training details and parameters are reported in Appendix 4.1. D. Frequency
$\bar{\omega}$ minimizing the error in the real and imaginary part of
$\tilde{\lambda}_{\pm}$ as a function of $g$, for the three training
algorithms. Solid orange lines in the middle panel show theoretical prediction
for ridge regression (see Appendix 4.10).

As a second way to characterize performance, like in Fig. 1, we considered
feedback networks trained via Ridge regression [18]. In this case, the readout
vector is deterministic; in the Fourier space, it can be expressed as (see
Appendix 4.10):

$\tilde{\mathbf{n}}_{R}=\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}\tilde{P}\begin{pmatrix}1\\\
0\end{pmatrix}$ (29)

where $\tilde{P}=\left(C^{R}\right)^{-1}+N\sigma^{2}\mathbf{I}$ and
$\mathbf{I}$ is the $2\times 2$ identity matrix. As in the case of noisy LS
regression, also this readout vector generates only two outlier eigenvalues
$\tilde{\lambda}_{\pm}$, whose values can be again computed through Eq. (26);
in this case a closed-form expression can be computed (Appendix 4.10).

For moderate values of the regularization parameter $\sigma$, the resulting
$\tilde{\lambda}_{\pm}$ are complex conjugates that deviate somewhat from
${\lambda}_{\pm}$ (see Supp. Fig. 10 for the full bifurcation diagram).
Specifically, their real part is always smaller than 1, implying that the
resulting autonomous dynamics are always stable (see Appendix 4.10). The
amplitude of the mismatch between the real and the imaginary parts of
$\tilde{\lambda}_{\pm}$ and the target eigenvalues ${\lambda}_{\pm}$ depends
both on $g$ and $\omega$ (see Appendix 4.10), and is minimized at a finite
frequency $\bar{\omega}$ which monotonically decreases with increasing $g$
(Fig. 5D center, solid lines). Importantly, the value of $\bar{\omega}$ is
predicted to behave similarly (although not identically) to $\omega^{*}$. Fig.
5D (middle) shows an excellent match between these predictions and simulation
results.

As a third and final example, we considered linear networks trained via the
RLS algorithm [19, 7]. In this case, an analytical description of the closed-
loop spectrum and resulting dynamics is much harder to obtain; we thus
computed the value of the preferred frequency $\bar{\omega}$ from simulations.
We found that the mismatch between $\tilde{\lambda}_{\pm}$ and
${\lambda}_{\pm}$ displays a strong, non-monotonic dependence on the target
frequency (Fig. 5C, right); the preferred frequency $\bar{\omega}$ is, again,
quite close to $\omega^{*}$ (Fig. 5D, right).

To conclude, we analyzed performance in linear feedback networks; as for non-
linear networks (Fig. 1), we found that performance is maximized for a
preferred frequency $\bar{\omega}$ which decreases with the connectivity
strength $g$. Analysing how the simple LS readout solution interacts with
noise, we predicted that the preferred frequency $\bar{\omega}$ is expected to
lay close to $\omega^{*}$, i.e. the resonance frequency where encoding
dynamics has maximal dimensionality and is minimally synchronized. This
prediction is exactly verified in networks trained via LS regression, but also
carry over in a qualitative fashion to networks trained via different training
algorithms. In fact, we showed that different algorithms are affected by
different kinds of biases, whose effect is to shift the value of the preferred
frequency $\bar{\omega}$ away from $\omega^{*}$ without changing its overall
qualitative behaviour.

###  2.6 Internal representation in non-linear networks

We finally turn back to the original problem of analyzing training performance
in non-linear feedback networks (Fig. 1). Our analysis of linear networks
revealed that a key feature which determines training performance is the
quality of representation of the target signal within the reservoir. This
representation can be characterised by its dimensionality or, equivalently, by
the synchrony of activity across units in the reservoir.

Guided by these insights, we examined the properties of open-loop dynamics
(Eq. (1)) in non-linear networks. Because of the non-linearity, the neural
trajectory $\mathbf{x}(t)$ is in this case not planar, but curved along many
dimensions (Fig. 6A); most of its variance, however, is still explained by two
directions (Fig. 6B). We investigated numerically the properties of non-linear
target representations by using the same measures as for linear networks,
namely the dimensionality and the spread of phases across units. Although in
non-linear systems these are not equivalent measures, we find that their
behaviour is qualitatively similar to one another, and to the behaviour of
their analogues in linear systems (Fig. 6C-D left). First, both measures
increase monotonically with the connectivity strength $g$. Second, for any
fixed value of $g$, both measures display a maximum at an intermediate
frequency $\omega^{*}$.

In the middle panels of Fig. 6C-D, we display the resonance frequency
$\omega^{*}$ computed from both measures of non-linear representations (left
panels) across various values of $g$ and for three target amplitudes (legend).
As in the linear case, we find that the value of $\omega^{*}$ decreases with
the connectivity strength $g$; unlike the linear case, however, it depends on
the target amplitude $A$. For small target amplitudes (light gray), both
measures of non-linear representations yield values of $\omega^{*}$ which are
quantitatively very close to the values predicted by the linear theory, i.e.
$\sqrt{1-g^{2}}$ (yellow line). This is expected, as for low-amplitude driving
the reservoir activity mostly remains in the vicinity of the origin, a region
where the non-linear dynamics are approximately linear. In the nonlinear case,
however, as the target amplitude $A$ increases (darker shades of gray, see
legend), the resonance frequency $\omega^{*}$ also increases. The decrease of
$\omega^{*}$ with $g$ is retained, but to a lesser extent.

In the right panels of Fig. 6C-D, we compare the resonance frequency
$\omega^{*}$ predicted from analyzing non-linear representations to the
preferred frequency $\bar{\omega}$ which minimizes training performance (Fig.
1). Although the two quantities do not exactly coincide, they display
significant correlations. Remarkably, the value of $\omega^{*}$ correctly
captures the behaviour of the preferred frequency $\bar{\omega}$ with the
target amplitude $A$: like $\omega^{*}$, $\bar{\omega}$ increases with $A$, as
can be seen by the clustering of different shades of grey in Fig. 6C-D, right
panels.

Importantly, this observation is not sensitive to training details. It is
consistent across the three training algorithms we used (Fig. 6C-D, right
panels) and across a broad range of training hyper-parameters (Supp. Fig. 7).
These results suggest that the properties of the internal representation –
here measured by the dimensionality and synchronization of the open-loop
dynamics – play a crucial role in determining training performance of non-
linear networks, as they do for linear networks.

Figure 6:  Using open-loop reservoir dynamics to predict training performance
in non-linear networks. A. Reservoir trajectories $\mathbf{x}(t)$ in driven
non-linear networks: example trajectory projected onto the first three PCs of
network activity (note scale of PC3 axis). B. Variance explained by projecting
trajectories on the first two PC axes (note scale of variance). C. Using
representation dimensionality to predict training performance. Left:
participation ratio of driven trajectory as a function of $\omega$ for a range
of $g$ values (blue shades). Results are averages over 20 simulations of
networks of size $N=2000$, with $A=1$. Center: resonance frequency
$\omega^{*}$, measured as the position of maximum dimensionality, as a
function of $g$ for three different values of $A$ (grey shades). Right: error-
minimizing frequency $\bar{\omega}$ (from Fig. 1) plotted against $\omega^{*}$
(from center panel). Results are shown for three training algorithms (legend).
D. Using phase spread of driven trajectories to predict performance. Unit
activities $x_{i}(t)$ were fitted with sinusoidal functions of the driving
frequency $\omega$, and the variance of the phase distribution was measured.
Left, center and right panels are the same as in C.

##  3 Discussion

Ubiquitously across biology, complex high-dimensional systems interact with
their environment through low-dimensional channels. The computational
modelling of such setups has advanced considerably in the past two decades
with the emergence of reservoir computing techniques [3, 26], where learning
acts on such low-dimensional bottlenecks. Despite the simplicity of this
learning scheme, the factors contributing to or hindering the success of
training in reservoir networks are in general not well understood [17]. In
particular, a theory is lacking for predicting – based on the characteristics
of the reservoir and the target function – dynamics and performance of trained
feedback networks.

In this work, we studied learning performance of feedback networks trained to
self-sustain a sinusoidal readout signal. Through mathematical analysis, we
showed that learning performance is mostly controlled by the quality of the
internal representation of the target signal. This quality can be quantified
by analyzing the open-loop dynamics and measuring the condition number of
their cross-correlation matrix, a number that characterizes to what extent the
network dynamics is robust to training noise. We found that the condition
number displays, like training performance, a complex dependence on the
parameters controlling the reservoir internal properties (strength of
reservoir connectivity $g$) and the readout target function (frequency
$\omega$). The parameter values where the condition number is minimized,
$\omega^{*}=\sqrt{1-g^{2}}$, define an optimal spot for learning. At this
optimal point, internal representations are characterized by maximal
dimensionality and minimal synchrony, which are two ways of quantifying the
richness of the dynamic repertoire available to the learning algorithm. Our
insights were derived by studying linearized dynamics and were later tested on
non-linear networks, where they successfully capture non-trivial aspects of
training performance.

The condition number of the cross-correlation matrix has been pointed out in
several studies as a key quantity in determining performance [17, 27]. Our
work analytically quantifies those empirical observations in the framework of
networks trained on a simple task via common LS-based algorithms. We have
shown, however, that performance might depend on other features, such as
closed-loop stability, for other non-standard algorithms (see Fig. 4).

Importantly, our analysis differentiates between two properties that might
hinder network performance: high-norm readouts and non-normality. In a number
of classic studies [7, 17, 12], large norms of the readout vector have been
associated with impaired performance. In addition, recent observations
indicate that training performance is low in parameter regions where the open-
loop dynamics is highly non-normal [13], and link low performance to large
readout vectors. In our framework, the two properties can be analyzed
separately. Non-normality can be measured from the angle $\theta$ between the
two activity eigenvectors $\mathbf{v}_{\pm}$; Fig. 2C indicates that non-
normality is minimal at the resonance frequency $\omega^{*}$. The norm of the
readout vector $\mathbf{n}_{LS}$ can be instead derived from Eq. (16) (see
Appendix 4.7); we show in Supp. Fig. 9 that, for every value of connectivity
$g$, the norm of the readout vector is monotonic in the target frequency
$\omega$. We conclude that these two quantities are not equivalent predictors
of learning performance; in our setting, training performance is optimal close
to $\omega^{*}$, so that non-normality is identified as the dominating factor
in controlling performance.

Several studies have supported the hypothesis that learning capability is
maximized in the parameter region where dynamics is close to the boundary
between ordered and chaotic activity, i.e. the _edge-of-chaos_ [14, 7, 6]. Our
findings are consistent with this hypothesis: we have shown that the condition
number (and, consequently, the training error) monotonically decreases as the
strength of reservoir connectivity $g$ is increased from 0 towards its
critical value. However, our analysis has shown that, together with the
strength of internal connectivity, learning performance is crucially shaped by
the properties of the target function. By analysing non-linear networks,
furthermore, we have found that the parameter region characterized by
maximally high-dimensional and de-synchronized internal representations does
not necessarily coincide with the edge-of-chaos; the two regions in fact
diverge as the target amplitude $A$ is increased and activity becomes strongly
non-linear (Figs. 6 and 11). Specifically, as $A$ increases, the critical
frequency where activity becomes chaotic moves to very high values [25] (Fig.
11), while the resonance frequency $\omega^{*}$ (which measures activity
dimensionality and synchrony) remains close to the training-preferred
frequency $\bar{\omega}$ (Fig. 6). This result suggests that future research
should focus on characterizing the properties of driven non-linear activity
rather than analysing the transition to chaos per-se.

The numerical analysis of non-linear networks (Fig. 6), which was led by the
insights gained from the linear theory, suggests that representation quality
is a major determinant of closed-loop performance also in the case of non-
linear networks. Exploiting the link between the two, we were able to predict
the dependence of the preferred frequency $\bar{\omega}$ on both the internal
connectivity $g$ and the target amplitude $A$, which plays no role in the
linear counterpart. This is despite the fact that the non-linearity of the
dynamics introduces, in trained networks, new qualitative behaviours which do
not exist in linear networks. In particular, we observe that the training
error (and, consequently, the value of $\bar{\omega}$) strongly depends on the
hyper-parameters controlling the stability of the limit cycle which
constitutes the internal representation (see Appendix 4.1 and Supp. Fig. 8).
In this respect, a more detailed analysis is called for; we hope that future
work would extend our analytic framework to cover non-linear reservoirs.

## References

  * [1] R. Chaudhuri, B. Gercek, B. Pandey, A. Peyrache, and I. Fiete.  The intrinsic attractor manifold and population dynamics of a canonical cognitive circuit across waking and sleep.  Nat. Neurosci., 22(9):1512–1520, 2019.
  * [2] R. J Low, S. Lewallen, D. Aronov, R. Nevers, and D. Tank.  Probing variability in a cognitive map using manifold inference from neural dynamics.  bioRxiv, 10.1101/418939, 2018.
  * [3] H. Jaeger and H. Haas.  Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication.  Science, 304(5667):78–80, 2004.
  * [4] W. Maass, P. Joshi, and E. D Sontag.  Computational aspects of feedback in neural circuits.  PLOS Comput. Biol., 3(1):1–20, 01 2007.
  * [5] H. Jaeger.  The “echo state” approach to analysing and training recurrent neural networks - with an erratum note.  GMD Report, 2001.
  * [6] M. Lukoševičius and H. Jaeger.  Reservoir computing approaches to recurrent neural network training.  Comput. Sci. Rev., 3(3):127 – 149, 2009.
  * [7] D. Sussillo and L.F. Abbott.  Generating coherent patterns of activity from chaotic neural networks.  Neuron, 63(4):544 – 557, 2009.
  * [8] H. Jaeger.  Adaptive nonlinear system identification with echo state networks.  NIPS, pages 609–616, 2002.
  * [9] M. Massar and S. Massar.  Mean-field theory of echo state networks.  Phys. Rev. E, 87:042809, 2013.
  * [10] A. Rivkind and O. Barak.  Local dynamics in trained recurrent neural networks.  Phys. Rev. Lett., 118:258101, 2017.
  * [11] F. Mastrogiuseppe and S. Ostojic.  A geometrical analysis of global stability in trained feedback networks.  Neural Computat., 31(6):1139–1182, 2019.
  * [12] F. Schuessler, A. Dubreuil, F. Mastrogiuseppe, S. Ostojic, and O. Barak.  Dynamics of random recurrent networks with correlated low-rank structure.  Phys. Rev. Research, 2:013111, 2020.
  * [13] L. Logiaco, L. Abbott, and S. Escola.  A model of flexible motor sequencing through thalamic control of cortical dynamics.  bioRxiv, 2019.
  * [14] R. Legenstein and W. Maass.  Edge of chaos and prediction of computational performance for neural circuit models.  Neural Networks, 20(3):323 – 334, 2007.
  * [15] H. Sompolinsky, A. Crisanti, and H. J. Sommers.  Chaos in random neural networks.  Phys. Rev. Lett., 61:259–262, 1988.
  * [16] N. Bertschinger and T. Natschläger.  Real-time computation at the edge of chaos in recurrent neural networks.  Neural Computat., 16(7):1413–1436, 2004.
  * [17] H. Jaeger.  Reservoir riddles: Suggestions for echo state network research (extended abstract).  volume 3, pages 1460 – 1462 vol. 3, 01 2005.
  * [18] Arthur E. Hoerl and Robert W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics, 12(1):55–67, 1970.
  * [19] W. Liu, J. Principe, and S. Haykin.  Kernel Adaptive Filtering: A Comprehensive Introduction.  Wiley Publishing, 2010.
  * [20] V. L. Girko.  Circular law.  Theory Probab. Appl., 29(4):694–706, 1985.
  * [21] P. Gao and S. Ganguli.  On simplicity and complexity in the brave new world of large-scale neuroscience.  Curr. Opin. Neurobiol., 32:148–55, 2015.
  * [22] A. Litwin-Kumar, K. Harris, R. Axel, H. Sompolinsky, and L. Abbott.  Optimal degrees of synaptic connectivity.  Neuron, 93(5):1153 – 1164.e7, 2017.
  * [23] Y Nakatsukasa.  The low-rank eigenvalue problem.  arXiv, 1905.11490, 2019.
  * [24] V. Aalo, G. Efthymoglou, and C. Chayawan.  On the envelope and phase distributions for correlated gaussian quadratures.  IEEE Comm. Lett., 11(12):985–987, 2007.
  * [25] K. Rajan, L. F. Abbott, and H. Sompolinsky.  Stimulus-dependent suppression of chaos in recurrent neural networks.  Phys. Rev. E, 82:011903, 2010.
  * [26] W. Maass, T. Natschläger, and H. Markram.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  Neural Computat., 14(11):2531–2560, 2002.
  * [27] M. Lukosevicius.  A Practical Guide to Applying Echo State Networks. Neural Networks: Tricks of the Trade. Lecture Notes in Computer Science. Springer, 2012.
  * [28] J. Eckmann and D. Ruelle.  Ergodic theory of chaos and strange attractors.  In The theory of chaotic attractors, pages 273–312. Springer, 1985\.
  * [29] K. Geist, U. Parlitz, and W. Lauterborn.  Comparison of different methods for computing lyapunov exponents.  Prog. Theor. Phys., 83(5):875–893, 1990.

### Acknowledgments

This work was supported in part by the Israeli Science Foundation (grant
number 346/16, OB; and grant number 155/18, NB). FM would like to thank Haim
Sompolinsky and Ran Rubin for their supervision on an early version of the
project. FM would also like to thank Laureline Logiaco for useful discussions,
and Srdjan Ostojic and Manuel Beiran for their feedback on a previous version
of the manuscript. LS would like to thank Friedrich Schuessler for helpful
discussions.

##  4 Appendix

###  4.1 Training of feedback networks

In the following, we report the procedures used to train feedback
architectures (Figs. 1 and 5). Procedures are detailed for the general case of
non-linear networks; the case of linear networks corresponds to taking
$\Phi(x)=x$. Results are averages across 1000 different network and training
realizations.

#### LS regression training

Training is performed in the open-loop setup. In a first phase, open-loop
activity (where we enforce $u(t)=A\cos(\omega t)$ in Eq. (1)) is simulated by
using the Scipy odeint routine from $t=0$ to $t=T^{\text{tot}}$, with
$T=N^{\text{tot}}2\pi/\omega$. Activity is stored in a $L\times N$ matrix
$\bm{\Phi}$, where $L$ indicates the number of time points used for
integration. The $L$-dimensional vector $\bm{F}$ is constructed by computing
the target function $f(t)$ at the same time points. Activity and target
function from $t=0$ to $T^{\text{tr}}=N^{\text{tr}}2\pi/\omega$ are later
discarded, resulting in $L^{\prime}\times N$ and $L^{\prime}\times 1$ matrices
$\bm{\Phi}$ and $\bm{F}$, where $L^{\prime}$ indicates the number of time
points kept after discarding the transient. We used $N^{\text{tot}}=20$ and
$N^{\text{tr}}=8$. In order to regularize the cross-correlation matrix and to
ease local stability in non-linear networks, white noise is then added on top
of activity: $\bm{\tilde{\Phi}}=\bm{\Phi}+\sigma^{\text{LS}}\bm{\xi}$, where
$\bm{\xi}$ is a $L^{\prime}\times N$ matrix of standard Gaussian variables.
The trained readout vector $\mathbf{n}$ is finally computed as:

$\mathbf{n}=(\bm{\tilde{\Phi}}^{\top}\bm{\tilde{\Phi}})^{-1}\bm{\tilde{\Phi}}^{\top}\bm{F}.$
(30)

In linear networks, training performance is measured in terms of the mismatch
between the target outlier eigenvalues $\lambda_{\pm}$ (see Section 2.4) and
the outlier eigenvalues $\tilde{\lambda}_{\pm}$, defined as the pair of
complex conjugate eigenvalues of
$\bar{\mathbf{J}}=\mathbf{J}+\mathbf{m}\mathbf{n}^{\top}$ whose real part is
maximally close to one. In non-linear networks, performance is measured on
closed-loop activity. To this end, closed loop dynamics (Eq. (3)) is simulated
from $t=0$ to $t=T^{\text{tot}}$. The initial condition is taken to be equal
to activity in the last time step of the open-loop simulation; on top of it,
an $N$-dimensional vector of white noise of amplitude $\sigma^{\text{pert}}A$
is added. The latter perturbation was used to take into account training error
generated by unstable local dynamics; we take $\sigma^{\text{pert}}=0$ in
linear networks. To measure the test error we fitted a sinusoidal function
$F(t)$ of fixed amplitude $A$ and frequency $\omega$ to the readout signal
$\mathbf{z}=\mathbf{n}^{\top}\bm{\Phi}^{\top}$ obtained in the closed-loop
simulation, yielding a novel $L$-dimensional vector $\bm{F}$. Readout error is
finally measured as $\langle\>|z_{k}-F_{k}|\>\rangle_{k}$, where the average
is taken over all the integration time points from $t=0$ to
$t=T^{\text{tot}}$. If the fit fails, we set the readout error to 1.
Parameters used in Fig. 1 are $N=400$, $\sigma^{\text{LS}}=0.01$ and
$\sigma^{\text{pert}}=0.1$. Parameters used in Fig. 5 are $N=400$,
$\sigma^{\text{LS}}=0.01$ and $\sigma^{\text{pert}}=0$.

#### Ridge regression training

As in the LS case, training is performed in the open-loop setup. The
$L^{\prime}\times N$ open-loop activity $\bm{\Phi}$ matrix is obtained as
above. The trained readout vector $\mathbf{n}$ is then computed as:

$\mathbf{n}=(\bm{{\Phi}}^{\top}\bm{{\Phi}}+(\sigma^{\text{R}})^{2}\mathbf{I})^{-1}\bm{{\Phi}}^{\top}\bm{F}$
(31)

where $\mathbf{I}$ indicates the $N$-dimensional identity matrix. Training
performance is measured as in the LS case. Parameters used in Fig. 1 are
$N=400$, $\sigma^{\text{R}}=1$ and $\sigma^{\text{pert}}=0.1$. Parameters used
in Fig. 5 are $N=400$,
$\left(\sigma^{\text{R}}\right)^{2}=L^{\prime}/2\cdot\sigma^{2}$ and
$\sigma^{\text{pert}}=0$; $\sigma^{2}=10^{-7}$ is the regularization parameter
used for regression in the Fourier space (see Appendix 4.10), which was used
to compute the theoretical prediction for outlier eigenvalues.

#### RLS training

Training is performed in the closed-loop setup, from $t=0$ to
$t=T^{\text{tot}}$, with $T=N^{\text{tot}}2\pi/\omega$ and $N^{\text{tot}}=20$
(Fig. 1) or 10 (Fig. 5). At $t=0$, an $N\times N$-dimensional matrix
$\mathbf{P}$ is initialized as: $\bm{P}=\mathbf{I}/\alpha$, where $\mathbf{I}$
indicates the $N$-dimensional identity matrix and $\alpha$ is a free
parameter. Matrix $\mathbf{P}$ represents a running estimate of the inverse of
the activity cross-correlation matrix [7]. Readout vector $\mathbf{n}$ is
further initialized with zero entries. At every learning step, closed-loop
activity is simulated from $t_{0}$ to $t_{0}+\tau$ (Eq. (3)), with
$\tau=(2\pi/\omega)/500$. Activity at $t=t_{0}+\tau$ is stored in an
$N$-dimensional vector $\bm{\Phi}$. Matrix $\mathbf{P}$ is then updated as
[7]:

$\mathbf{P}\leftarrow\mathbf{P}-\frac{\mathbf{P}\bm{\Phi}^{\top}\bm{\Phi}\mathbf{P}}{1+\bm{\Phi}^{\top}\mathbf{P}\bm{\Phi}}$
(32)

The readout vector $\mathbf{n}$ is then updated as:

$\mathbf{n}\leftarrow\mathbf{n}-e\mathbf{P}\bm{\Phi}$ (33)

where the error $e$ is measured as: $e=z(t_{0}+\tau)-f(t_{0}+\tau)$. Once
training is completed, performance is measured as in the LS case. Parameters
used in Fig. 1 are $N=400$, $\alpha=1$ and $\sigma^{\text{pert}}=0.1$.
Parameters used in Fig. 5 are $N=400$, $\alpha=1$ and
$\sigma^{\text{pert}}=0$.

Figure 7:  Emergence of preferred frequency in non-linear feedback networks,
supplementary results obtained for different hyper-parameters. A. LS training.
Parameters are as in Fig. 1, except $\sigma^{\text{pert}}=0.01$ (top) and
$\sigma^{\text{LS}}=0.001$ (bottom). B. Ridge training. Parameters are as in
Fig. 1, except $\sigma^{\text{pert}}=0.01$ (top) and $\sigma^{\text{R}}=0.1$
(bottom). C. RLS training. Parameters are as in Fig. 1, except
$N^{\text{tot}}=10$ (top and bottom), $\sigma^{\text{pert}}=0.01$ (top) and
$\alpha=5$ (bottom). In the bottom row, we have removed from the plot the
preferred frequencies $\bar{\omega}$ at low $g$ values in the cases where
training fails (i.e. it is characterized by very high error) for every value
of $\omega$ tested.  Figure 8:  Emergence of preferred frequency in non-linear
feedback networks, example trials. Example trials as in Fig. 1B where training
is performed via Ridge regression (A) or RLS (B). Parameters are as in Fig.
1B, except $N^{\text{tot}}=3$ in B.

###  4.2 Analysis of linear open-loop reservoirs

For a general input $f(t)$, the system of linear equations Eq. (4) admits the
asymptotic solution ($t\rightarrow\infty$)

$\mathbf{x}(t)=\int_{0}^{t}e^{\mathbf{(J-I)}(t-\tau)}\>\mathbf{m}f(\tau)\>\mathop{}\\!\mathrm{d}\tau,$
(34)

which for a complex exponential $f(t)=e^{st}$ for $s\in\mathbb{C}$, simplifies
to

$\mathbf{x}(t)=[(1+s)\mathbf{I}-\mathbf{J}]^{-1}\mathbf{m}e^{st}.$ (35)

In particular, for $f(t)=\cos(\omega t)=\frac{1}{2}(e^{i\omega t}+e^{-i\omega
t})$, one finds the expression in Eq. (5) in the main text. Note that
$\mathbf{x}_{\pm}$ defined in Eq. (5) correspond to the amplitude of the peaks
of the Fourier transform of $\mathbf{x}(t)$; indeed, we have:

$\begin{split}\mathbf{\hat{x}}(\hat{\omega})&=\int_{-\infty}^{\infty}\mathbf{x}(t)e^{-i\hat{\omega}t}\mathop{}\\!\mathrm{d}t=\frac{1}{2}[\mathbf{x_{+}}\delta(\hat{\omega}-\omega)+\mathbf{x_{-}}\delta(\hat{\omega}+\omega)].\end{split}$
(36)

In deriving Eq. (7), we defined:

$\begin{split}\mathbf{v}_{+}&:=\frac{1}{2}\left(\mathbf{x}_{+}+\mathbf{x}_{-}\right)\\\
\mathbf{v}_{-}&:=\frac{i}{2}\left(\mathbf{x}_{+}-\mathbf{x}_{-}\right).\end{split}$
(37)

###  4.3 Statistics of spanning vectors $\mathbf{v}_{+}$ and $\mathbf{v}_{-}$

In this section, we characterize the geometry of vectors $\mathbf{v}_{+}$ and
$\mathbf{v}_{-}$ in terms of their norms and overlap.

We start by evaluating the dot product:

$\left[\left(\mathbf{I}-\frac{\mathbf{J}}{a}\right)^{-1}\mathbf{m}\right]\cdot\left[\left(\mathbf{I}-\frac{\mathbf{J}}{b}\right)^{-1}\mathbf{m}\right]=\mathbf{m}^{\top}\left[\left(\mathbf{I}-\frac{\mathbf{J}}{a}\right)^{-1}\right]^{\top}\left(\mathbf{I}-\frac{\mathbf{J}}{b}\right)^{-1}\mathbf{m}$
(38)

with $a$, $b\in\mathbb{C}$. If the eigenvalues of $\mathbf{J}/a$ and
$\mathbf{J}/b$ have real part smaller than one, we can use the power series
expansion

$\left(\mathbf{I}-\frac{\mathbf{J}}{a}\right)^{-1}=\sum_{p=0}^{\infty}\frac{\mathbf{J}^{p}}{a^{p}},$
(39)

so that Eq. (38) becomes

$\mathbf{m}^{T}\sum_{p=0}^{\infty}\sum_{q=0}^{\infty}\frac{\left(\mathbf{J}^{p}\right)^{\top}\mathbf{J}^{q}}{a^{p}b^{q}}\mathbf{m}.$
(40)

Since $\mathbf{J}$ is random, the value of this expression randomly fluctuates
across different realizations of recurrent connectivity $\mathbf{J}$. We thus
turn to a statistical characterization, and evaluate Eq. (38) by computing its
mean and variance with respect to different realizations of $\mathbf{J}$.

The mean yields, to the leading order in $N$ [12]:

$\begin{split}&\mathbf{m}^{T}\sum_{p=0}^{\infty}\sum_{q=0}^{\infty}\frac{\left\langle(\mathbf{J}^{p})^{\top}\mathbf{J}^{q}\right\rangle_{\mathbf{J}}}{a^{p}b^{q}}\mathbf{m}\\\
&=\mathbf{m}^{T}\mathbf{m}\sum_{p=0}^{\infty}\left(\frac{g^{2}}{ab}\right)^{p}\\\
&=N\sum_{p=0}^{\infty}\left(\frac{g^{2}}{ab}\right)^{p}\\\
&=N\frac{1}{1-\frac{g^{2}}{ab}}=N\frac{ab}{ab-g^{2}}.\end{split}$ (41)

We have used:

$\left\langle(\mathbf{J}^{p})^{T}\mathbf{J}^{q}\right\rangle_{\mathbf{J}}=\delta_{pq}\>\mathbf{I}\>{g^{2p}}$
(42)

which comes from observing that $\mathbf{J}^{p}$ is a random matrix, which is
uncorrelated to $\mathbf{J}^{q}$ for $q\neq p$ has variance $g^{2p}/N$. This
yields:

$\sum_{k=1}^{N}\langle(J^{p})_{ki}(J^{q})_{kj}\rangle_{\mathbf{J}}=\sum_{k=1}^{N}\delta_{ij}\delta_{pq}\frac{g^{2p}}{N}=\delta_{ij}\delta_{pq}g^{2p},$
(43)

from which we obtain Eq. (42).

The variance can be computed in a similar way. Like the mean, the variance is
characterized by $\mathcal{O}(N)$ scaling [12]; as a consequence, variability
due to different realizations of $\mathbf{J}$ does not enter the dot product
Eq. (38) to the leading order in $N$, and dot products can be replaced with
their mean (Eq. (41)) when $N\rightarrow\infty$.

We can now compute the mean norm of the spanning vectors from combining Eqs.
(37) and (41):

$\begin{split}\frac{\|\mathbf{v}_{+}\|^{2}}{N}&=\frac{1}{4N}\left(\mathbf{x}_{+}\cdot\mathbf{x}_{+}+2\mathbf{x}_{+}\cdot\mathbf{x}_{-}+\mathbf{x}_{-}\cdot\mathbf{x}_{-}\right)\\\
&=\frac{1}{4}\left(\frac{1}{(1+i\omega)^{2}-g^{2}}+\frac{2}{1+\omega^{2}-g^{2}}+\frac{1}{(1-i\omega)^{2}-g^{2}}\right)\\\
&=\frac{(1-g^{2})^{2}+\omega^{2}(2-(1-g^{2}))}{((1-g^{2}-\omega^{2})^{2}+4\omega^{2})(1-g^{2}+\omega^{2})};\end{split}$
(44)

the calculation of the norm of $\mathbf{v}_{-}$ is very similar, and only
differs in the sign of the first summand in the second line above:

$\begin{split}\frac{\|\mathbf{v}_{-}\|^{2}}{N}&=\frac{1}{4}\left(-\frac{2(1-g^{2}-\omega^{2})}{(1-g^{2}-\omega^{2})^{2}+4\omega^{2}}+\frac{2}{1-g^{2}+\omega^{2}}\right)\\\
&=\frac{\omega^{2}(\omega^{2}+2-(1-g^{2}))}{((1-g^{2}-\omega^{2})^{2}+4\omega^{2})(1-g^{2}+\omega^{2})}.\end{split}$
(45)

Finally, when computing the dot product between the two vectors, the cross
terms cancel to yield

$\begin{split}\frac{\mathbf{v}_{+}\cdot\mathbf{v}_{-}}{N}&=\frac{1}{4}\left(-\frac{1}{(1+i\omega)^{2}-g^{2}}+\frac{1}{(1-i\omega)^{2}-g^{2}}\right)\\\
&=\frac{\omega}{(1-g^{2}-\omega^{2})^{2}+4\omega^{2}}.\end{split}$ (46)

With these expressions, the angle $\theta$ between $\mathbf{v}_{+}$ and
$\mathbf{v}_{-}$ can be written as

$\text{cos}(\theta)=\frac{1-g^{2}+\omega^{2}}{\sqrt{(1-g^{2})^{2}+\omega^{2}(2-(1-g^{2}))}\sqrt{\omega^{2}+2-(1-g^{2})}},$
(47)

By denoting $\varepsilon=1-g^{2}$, we now summarize the statistics of the
spanning vectors $\mathbf{v}_{\pm}$:

$\boxed{\begin{split}&\frac{\|\mathbf{v}_{+}\|^{2}}{N}=\frac{\omega^{2}(2-\varepsilon)+\varepsilon^{2}}{(\varepsilon+\omega^{2})((\varepsilon-\omega^{2})^{2}+4\omega^{2})}\\\
&\frac{\|\mathbf{v}_{-}\|^{2}}{N}=\frac{\omega^{2}(2-\varepsilon+\omega^{2})}{(\varepsilon+\omega^{2})((\varepsilon-\omega^{2})^{2}+4\omega^{2})}\\\
&\frac{\mathbf{v}_{+}\cdot\mathbf{v}_{-}}{N}=\frac{\omega}{(\varepsilon-\omega^{2})^{2}+4\omega^{2}}\end{split}}$
(48)

and angle between the two spanning vectors reads:

$\boxed{\text{cos}(\theta)=\frac{\mathbf{v}_{+}\cdot\mathbf{v}_{-}}{\|\mathbf{v}_{+}\|\|\mathbf{v}_{-}\|}=\frac{\varepsilon+\omega^{2}}{\sqrt{2-\varepsilon+\omega^{2}}\sqrt{\varepsilon^{2}+\omega^{2}(2-\varepsilon)}}.}$
(49)

###  4.4 Analysis of geometric properties of driven trajectories

We can use the expressions computed in Appendix 4.3 to evaluate the
participation ratio $d$:

$d=\frac{(\nu_{1}+\nu_{2})^{2}}{\nu_{1}^{2}+\nu_{2}^{2}}$ (50)

where $\nu_{1}$ and $\nu_{2}$ are the eigenvalues of the reduced cross-
correlation matrix $C^{R}$ (Eq. (9)). We express the eigenvalues as

$\nu_{1,2}=\frac{1}{2}\left(\gamma\pm\sqrt{\Delta^{2}}\right),$ (51)

where $\gamma:=\text{Tr}(C^{R})$ and
$\Delta^{2}:=\text{Tr}^{2}(C^{R})-4\cdot\text{det}(C^{R})$. This yields:

$d=2\frac{\gamma^{2}}{\gamma^{2}+\Delta^{2}}$ (52)

from which Eq. (10) follows. Inserting the results from Appendix 4.3, we also
have:

$\boxed{d=\frac{\varepsilon^{2}-2\varepsilon\omega^{2}+4\omega^{2}+\omega^{4}}{\varepsilon^{2}+2\omega^{2}+\omega^{4}}.}$
(53)

We now consider some limiting behaviours of the quantities computed above. In
the high-frequency limit $\omega\to\infty$, the norms and dot product of
$\mathbf{v}_{\pm}$ vanish, whereas the overlap is given by

$\lim_{\omega\to\infty}\text{cos}(\theta)=\frac{1}{\sqrt{2-\varepsilon}}=\frac{1}{\sqrt{1+g^{2}}};$
(54)

the dimensionality $d$ tends to 1 in this limit.

In the limit of low frequency, for any $0<g<1$ the norms and dot product can
be evaluated directly by substituting $\omega=0$ to obtain

$\|\mathbf{v}_{-}\|=\mathbf{v}_{+}\cdot\mathbf{v}_{-}=0$ (55)

and

$\|\mathbf{v}_{+}\|=\frac{1}{\sqrt{2\left(1-g^{2}\right)}},$ (56)

while the overlap assumes the same non-zero value as in the other extreme:

$\text{cos}(\theta)=\frac{1}{\sqrt{1+g^{2}}}$ (57)

and the dimensionality is one.

In between the two extremes, the dot product, the participation ratio, and the
norm $\|\mathbf{v}_{-}\|$ obtain a maximum value at a finite frequency. For
any $g$, we compute the frequency yielding the minimum angle, by searching for
local extrema of Eq. (49):

$\frac{\partial\text{cos}^{2}(\theta)}{\partial\omega}=\frac{\partial}{\partial\omega}\frac{\left(\omega^{2}+\varepsilon\right)^{2}}{\left(2-\varepsilon\right)\left(\omega^{4}-2\omega^{2}\varepsilon+\varepsilon^{2}\right)+4\omega^{2}}=0,$
(58)

the numerator of which, after some algebra, can be written as
$8\omega\left(\omega^{4}-\varepsilon^{2}\right)$. The angle thus has a local
minimum at

$\omega^{*}=\sqrt{\varepsilon}=\sqrt{1-g^{2}}.$ (59)

Likewise, from Eq. (53) we have

$\frac{\partial
d}{\partial\omega}=4\frac{\omega(\varepsilon-1)\left[\omega^{4}-\varepsilon^{2}\right]}{(\varepsilon^{2}+2\omega^{2}+\omega^{4})^{2}},$
(60)

which, again, implies a nontrivial maximum of $d$ at
$\omega^{*}(g)=\sqrt{1-g^{2}}$.

Finally, we observe that $g=1$ is a singular point, as in the limit of low
frequency we have

$\lim_{\omega\to 0}\|\mathbf{v}_{+}\|=\lim_{\omega\to
0}\|\mathbf{v}_{-}\|=\lim_{\omega\to
0}\mathbf{v}_{+}\cdot\mathbf{v}_{-}=\infty.$ (61)

but the overlap vanishes

$\lim_{\omega\to 0}\text{cos}(\theta)=0,$ (62)

and the participation ratio attains its global maximum

$\lim_{\omega\to 0}d=2.$ (63)

Note that the limits $\lim_{g\to 1,\omega\to 0}d$ and $\lim_{g\to 1,\omega\to
0}\text{cos}(\theta)$ do not exist, since they depend on the order of limits
taken. To see this, compare Eqs. (57) and (62).

###  4.5 Distribution of response phases

The phase spread in response of different units (Eq. (15)) was computed as a
numerical integral performed over the probability distribution $p(\phi)$ whose
analytical form is available in [24]. We used:

$p(\phi)=\frac{1}{\pi\|\mathbf{v}_{+}\|\|\mathbf{v}_{-}\|\sqrt{1-\rho^{2}}}\left[\frac{\cos^{2}(\phi)}{\|\mathbf{v}_{+}\|^{2}(1-\rho^{2})}+\frac{\sin^{2}(\phi)}{\|\mathbf{v}_{-}\|^{2}(1-\rho^{2})}-\frac{2\rho\sin(\phi)\cos(\phi)}{\|\mathbf{v}_{+}\|\|\mathbf{v}_{-}\|(1-\rho^{2})}\right]^{-1}$
(64)

where $\rho=\cos(\theta)$. From [24] we also used:

$\bar{\phi}=\frac{1}{2}\arctan\left(\frac{2\rho\|\mathbf{v}_{+}\|\|\mathbf{v}_{-}\|}{\|\mathbf{v}_{+}\|^{2}-\|\mathbf{v}_{-}\|^{2}}\right).$
(65)

###  4.6 _From-k_ regression

In this section, we explain how _from-k_ least-squares regression (Fig. 4) is
performed.

For $2\leq k\leq N$, we define the cropped spanning vectors as

$\left[\mathbf{v}_{\pm}^{k}\right]_{i}=\begin{cases}\left[\mathbf{v}_{\pm}\right]_{i},&\text{if}\
1\leq i\leq k\\\ 0,&\text{otherwise}\end{cases}$ (66)

where $\left[\mathbf{v}_{\pm}\right]_{i}$ indicates the $i$-th element of
vectors $\mathbf{v}_{\pm}$. For every $k$, the from-k LS regressor is given by
the pseudo-inverse, which yields the LS regressor of minimum norm:

$\displaystyle\mathbf{n}^{k}_{LS}=\begin{pmatrix}\mathbf{v}_{+}^{k}&\mathbf{v}_{-}^{k}\end{pmatrix}\begin{pmatrix}\|\mathbf{v}_{+}^{k}\|^{2}&\mathbf{v}_{+}^{k}\cdot\mathbf{v}_{-}^{k}\\\
\mathbf{v}_{+}^{k}\cdot\mathbf{v}_{-}^{k}&\|\mathbf{v}_{-}^{k}\|^{2}\end{pmatrix}^{-1}\begin{pmatrix}1\\\
0\end{pmatrix}.$ (67)

Note that $\mathbf{n}_{LS}^{N}=\mathbf{n}_{LS}$.

In the following, we show that any _from-k_ readout vector
$\mathbf{n}_{LS}^{k}$ with $k<N$ is not fully contained in the plane spanned
by vectors $\mathbf{v}_{\pm}$. The _from-k_ readout vector
$\mathbf{n}_{LS}^{k}$ is contained in the plane spanned by cropped vectors
$\mathbf{v}^{k}_{\pm}$. Consider now any vector $\mathbf{a}$ which is
orthogonal to the reservoir trajectory plane spanned by $\mathbf{v}_{\pm}$. We
have

$0=\mathbf{a}\cdot\mathbf{v}_{\pm}=\sum_{i=1}^{N}(\mathbf{a})_{i}(\mathbf{v}_{\pm})_{i}.$
(68)

We have that cropped vectors $\mathbf{v}_{\pm}^{k}$ do overlap with vector
$\mathbf{a}$, because:

$\mathbf{a}\cdot\mathbf{v}_{\pm}^{k}=\sum_{i=1}^{k}(\mathbf{a})_{i}(\mathbf{v}_{\pm})_{i}=-\sum_{i=k+1}^{N}(\mathbf{a})_{i}(\mathbf{v}_{\pm})_{i}\neq
0.$ (69)

As a result, the readout vector $\mathbf{n}^{k}_{LS}$ also has a nonzero
overlap with $\mathbf{a}$.

###  4.7 Analysis of least-squares regression: norm

We analytically compute the norm of the least-squares readout solution
$\mathbf{n}_{LS}$ ($k=N$, Eqs. (18) and (19)). We start from Eq. (19) to
write:

$\displaystyle\mathbf{n}_{LS}$
$\displaystyle=\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}\left(C^{R}\right)^{-1}\begin{pmatrix}1\\\
0\end{pmatrix}$ (70)
$\displaystyle=\frac{1}{\|\mathbf{v}_{+}\|^{2}\|\mathbf{v}_{-}\|^{2}(1-\cos^{2}(\theta))}\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}\begin{pmatrix}\|\mathbf{v}_{-}\|^{2}&-\mathbf{v}_{+}\cdot\mathbf{v}_{-}\\\
-\mathbf{v}_{+}\cdot\mathbf{v}_{-}&\|\mathbf{v}_{+}\|^{2}\end{pmatrix}\begin{pmatrix}1\\\
0\end{pmatrix}$
$\displaystyle=\frac{1}{\|\mathbf{v}_{+}\|^{2}\left(1-\text{cos}^{2}(\theta)\right)}\left(\mathbf{v}_{+}-\frac{\mathbf{v}_{+}\cdot\mathbf{v}_{-}}{\|\mathbf{v}_{-}\|^{2}}\mathbf{v}_{-}\right)$
$\displaystyle:=\frac{\|\mathbf{v}_{+}\|\|\mathbf{v}_{-}\|^{2}}{\text{det}\left(C^{R}\right)}\left(\hat{\mathbf{v}}_{+}-\cos(\theta)\hat{\mathbf{v}}_{-}\right)$

The vector within the parenthesis above is contained in the activity-spanning
plane and is orthogonal to $\mathbf{v}_{-}$, and has norm

$\|\mathbf{v}_{+}\|\sqrt{1-\text{cos}^{2}(\theta)};$

the norm of $\mathbf{n}_{LS}$ is therefore given by

$\|\mathbf{n}_{LS}\|=\frac{1}{N\|\mathbf{v}_{+}\|\sqrt{1-\text{cos}^{2}(\theta)}}$
(71)

We find that this expression is monotonically increasing in both $g$ and
$\omega$, as shown in Supp. Fig. 9.

Figure 9:  Norm of LS-trained readout vector $\mathbf{n}$. Computed via Eq.
(71) for a range of values $g$ and $\omega$.

###  4.8 Analysis of least-squares regression: outlier eigenvalues

In this section, we compute the outlier eigenvalues of $\bar{\mathbf{J}}$
which result from full LS regression ($k=N$, Eqs. (18) and (19)). As derived
in the main text, outliers obey:

$\begin{split}1&={\mathbf{n}_{LS}}^{\top}[({\lambda}\mathbf{I}-\mathbf{J})^{-1}\mathbf{m}]:={\mathbf{n}_{LS}}^{\top}{\mathbf{x}}_{\lambda}.\end{split}$
(72)

Because of Eq. (17), we know that the equation above admits the solutions
$\lambda=\lambda_{\pm}=1\pm\omega$. In the following, we show that
$\lambda_{\pm}$ is in fact the only solution admitted. To this end, we use Eq.
(19) to rewrite the equation as:

$1=\begin{pmatrix}1&0\end{pmatrix}P^{\top}\begin{pmatrix}\mathbf{v}_{+}^{\top}\\\
\mathbf{v}_{-}^{\top}\end{pmatrix}\mathbf{x}_{\lambda}$ (73)

where we defined the short-hand notation $P:=\left(C^{R}\right)^{-1}$. A
little algebra yields:

$\begin{split}1&=\frac{1}{2}\begin{pmatrix}P_{11}&P_{21}\end{pmatrix}\begin{pmatrix}(\mathbf{x}_{+}+\mathbf{x}_{-})^{\top}\\\
i(\mathbf{x}_{+}-\mathbf{x}_{-})^{\top}\end{pmatrix}\mathbf{x}_{\lambda}\\\
&=\begin{pmatrix}P_{11}(\mathbf{x}_{+}+\mathbf{x}_{-})^{\top}+iP_{21}(\mathbf{x}_{+}-\mathbf{x}_{-})^{\top}\end{pmatrix}\mathbf{x}_{\lambda}\\\
&=\frac{1}{2}\left[\mathbf{x}_{+}^{\top}{\mathbf{x}}_{\lambda}(P_{11}+iP_{21})+\mathbf{x}_{-}^{\top}{\mathbf{x}}_{\lambda}(P_{11}-iP_{21})\right].\end{split}$
(74)

$P_{11}$ and $P_{21}$ are the elements of the first column of
$P=\left(C^{R}\right)^{-1}$, namely:

$\begin{split}&P_{11}=\frac{\|\mathbf{v}_{-}\|^{2}}{\|\mathbf{v}_{+}\|^{2}\|\mathbf{v}_{-}\|^{2}(1-\cos^{2}(\theta))}\\\
&P_{21}=-\frac{\mathbf{v}_{+}\cdot\mathbf{v}_{-}}{\|\mathbf{v}_{+}\|^{2}\|\mathbf{v}_{-}\|^{2}(1-\cos^{2}(\theta))}.\end{split}$
(75)

We can evaluate dot products in the form:

$\mathbf{x}_{\pm}^{\top}{\mathbf{x}}_{\lambda}=\left[({\lambda}_{\pm}\mathbf{I}-\mathbf{J})^{-1}\mathbf{m}\right]^{\top}\left[({\lambda}\mathbf{I}-\mathbf{J})^{-1}\mathbf{m}\right]$
(76)

by following Eq. (41), which was derived in Appendix 4.3 by averaging over the
random connectivity $\mathbf{J}$. This yields:

$\begin{split}&1=\frac{N}{2}\left[\frac{P_{11}+iP_{21}}{\lambda_{+}{\lambda}-g^{2}}+\frac{P_{11}-iP_{21}}{\lambda_{-}{\lambda}-g^{2}}\right]\end{split}$
(77)

which can be re-cast as a quadratic equation in $\lambda$:

$(1+\omega^{2}){\lambda}^{2}-\left[2g^{2}+N(P_{11}+\omega
P_{21})\right]{\lambda}+g^{4}+Ng^{2}P_{11}=0.$ (78)

As we know that the quadratic equation above is satisfied by
$\lambda=\lambda_{\pm}$, we conclude that Eq. (72) cannot admit other
solutions beyond these two.

###  4.9 Condition number of the cross-correlation matrix

The condition number of the reduced cross-correlation matrix $C^{R}$ is
defined as:

$c=\frac{\nu_{1}}{\nu_{2}}$ (79)

where $\nu_{1}$ and $\nu_{2}$ are the two eigenvalues of $C^{R}$. Their value
can be computed as a function of the statistics of the spanning vectors
$\mathbf{v}_{\pm}$, which in turns depend on $\omega$ and $g$ (see Appendix
4.3).

In this section, we show that for fixed $g$, the condition number $c$ is
minimized at the same value of $\omega$ which maximizes the participation
ratio $d$; this frequency coincides with $\omega^{*}$ (see Appendix 4.4).
Using:

$\nu_{1,2}=\frac{1}{2}\left(\gamma\pm\sqrt{\Delta^{2}}\right).$ (80)

we have:

$d=2\frac{\gamma^{2}}{\gamma^{2}+\Delta^{2}}=\frac{2}{1+(\frac{\Delta}{\gamma})^{2}},$
(81)

while

$c=\frac{\gamma+\Delta}{\gamma-\Delta}=1+\frac{2}{\frac{\gamma}{\Delta}-1}.$
(82)

The participation ratio $d$ is maximized when the quantity
$\frac{\Delta}{\gamma}$ is minimized, which is precisely where the condition
number $c$ is minimized.

###  4.10 Analysis of ridge regression

We start by computing the readout vector which performs ridge regression in
the Fourier space. The ridge regressor of Eq. (17) can be written in terms of
the $\mathbf{v}_{\pm}$ spanning vectors as [18]:

$\tilde{\mathbf{n}}_{R}=\begin{pmatrix}\mathbf{v}_{+}&\mathbf{v}_{-}\end{pmatrix}\tilde{P}\begin{pmatrix}1\\\
0\end{pmatrix}$ (83)

where $\tilde{P}=\left(C^{R}+N\sigma^{2}\mathbf{I}\right)^{-1}$ i.e.

$\tilde{P}=\frac{1}{\det\left(C^{R}+N\sigma^{2}\mathbf{I}\right)}\begin{pmatrix}\|\mathbf{v}_{-}\|^{2}+N\sigma^{2}&-\mathbf{v}_{+}\cdot\mathbf{v}_{-}\\\
-\mathbf{v}_{+}\cdot\mathbf{v}_{-}&\|\mathbf{v}_{+}\|^{2}+N\sigma^{2}\end{pmatrix}.$
(84)

This yields the ridge regressor readout:

$\tilde{\mathbf{n}}_{R}=\frac{\|\mathbf{v}_{-}\|^{2}\|\mathbf{v}_{+}\|}{\det\left(C^{R}+N\sigma^{2}\mathbf{I}\right)}\left[\left(1+\frac{N\sigma^{2}}{\|\mathbf{v}_{-}\|^{2}}\right)\mathbf{\hat{v}_{+}}-\cos(\theta)\mathbf{\hat{v}_{-}}\right],$
(85)

where $\mathbf{\hat{v}}$ indicates a normalized vector. By comparison with the
LS regressor, Eq. (70), we observe that $\sigma$ has two effects on the
readout: first, it reduces the norm, as expected from a regularizer. Second,
it biases the readout vector towards $\mathbf{v}_{+}$.

The outlier eigenvalues $\tilde{\lambda}_{\pm}$ imposed by the ridge regressor
can be found by utilizing the same strategy as in Appendix 4.8. We insert
$P_{11}=\tilde{P}_{11}$, $P_{21}=\tilde{P}_{21}$ and
$\lambda_{\pm}=1\pm\omega$ into Eq. (77) to obtain the equation for the
outlier eigenvalues $\tilde{\lambda}$:

$(1+\omega^{2})\tilde{\lambda}^{2}-\left(2g^{2}+N\frac{-\omega\mathbf{v}_{+}\cdot\mathbf{v}_{-}+\|\mathbf{v}_{-}\|^{2}+\sigma^{2}}{\det\left(C^{R}+N\sigma^{2}\mathbf{I}\right)}\right)\tilde{\lambda}+\left(g^{4}+g^{2}N\frac{\|\mathbf{v}_{-}\|^{2}+\sigma^{2}}{\det\left(C^{R}+N\sigma^{2}\mathbf{I}\right)}\right)=0.$
(86)

Depending on the values of $\sigma$, $\omega$ and $g$, the equation above
admits real or complex conjugate eigenvalues (see Supp. Fig. 10). For low
frequencies, the dynamics are characterised by two real eigenvalues
$\tilde{\lambda}_{\pm}$. As $\omega$ increases, a complex conjugate pair of
eigenvalues is formed. Importantly, their real part is always smaller than 1,
yielding stable closed-loop dynamics. To see this, we approximate the real
part of the solution to Eq. (86) by assuming that $\sigma\ll 1$:

$\mathcal{R}(\tilde{\lambda}_{+})\approx\frac{g^{2}}{1+\omega^{2}}+\frac{\omega\mathbf{v}_{+}\cdot\mathbf{v}_{-}+\|\mathbf{v}_{-}\|^{2}+\sigma^{2}}{2(1+\omega^{2})\left[\|\mathbf{v}_{+}\|^{2}\|\mathbf{v}_{-}\|^{2}\sin^{2}(\theta)+\sigma^{2}\left(\|\mathbf{v}_{+}\|^{2}+\|\mathbf{v}_{-}\|^{2}\right)\right]}.$
(87)

Now, by use of Eq. (48) we evaluate

$\|\mathbf{v}_{+}\|^{2}+\|\mathbf{v}_{-}\|^{2}=\frac{1}{1+\omega^{2}-g^{2}},$

implying that the pre-factor for the $\sigma^{2}$ term in the denominator is
always larger than 1; as a consequence, $\sigma>0$ always reduces the real
part.

To conclude, note that the analysis above allows to predict the behaviour of
outlier eigenvalues when ridge regression is performed in the Fourier space
(i.e. from the 2-dimensional system of equations in Eq. (17)). In Figs. 1 and
5, however, regression is performed in the temporal domain, on a higher-
dimensional ($L^{\prime}$-dimensional) set of equations (see Appendix 4.1). In
order to compare the analytical prediction with trained networks, in Fig. 5 we
thus scale the regularization parameter w.r.t. the value of $\sigma$ which is
used to derive analytical predictions, i.e. we set
$\left(\sigma^{\text{R}}\right)^{2}=\sigma^{2}\cdot L^{\prime}/2$ (see
Appendix 4.1).

Figure 10:  Bifurcation diagram for networks trained through ridge regression.
For a range of the regularization parameter $\sigma$ (grey shades), the curve
represents the location on the $g-\omega$ plane at which the outlier
eigenvalues $\tilde{\lambda}_{\pm}$ bifurcate from two real eigenvalues into a
pair or complex conjugate eigenvalues.

###  4.11 Characterization of the edge-of-chaos in non-linear networks

In analysing non-linear networks (Figs. 1 and 6), we varied the strength of
internal connectivity $g$ in such a way that the open-loop dynamics driven by
the target function $f(t)$ remains non-chaotic [15] for every value of the
forcing frequency tested. The critical value of connectivity strength $g^{c}$
at which open-loop dynamics becomes chaotic depends on the target frequency
$\omega$ and amplitude $A$ [25], and was investigated numerically (Fig. 11).

In order to find the critical values $g^{c}$, we start by computing the
Lyapunov dimension $d_{L}$ of driven activity [28, 29], which is defined based
on the Lyapunov spectrum $\Lambda=\\{\mu_{i}\\}_{i=1}^{N}$. Intuitively,
Lyapunov exponents $\mu_{i}$ quantify the rate of exponential convergence or
divergence of nearby trajectories along the different directions in state-
space. If we order the exponents such that
$\mu_{1}\geq\mu_{2}\geq\ldots\geq\mu_{N}$, the Kaplan-Yorke conjecture asserts
that the Lyapunov dimension of the attractor is given by the index $j$ for
which the number of contracting and expanding dimensions of the dynamics are
balanced, namely the index $j$ for which

$\sum_{i=1}^{j}\mu_{i}\geq 0,$

and

$\sum_{i=1}^{j+1}\mu_{i}<0.$

In order to allow for attractors of fractal dimension, the dimension is
defined as

$d_{L}=j+\frac{\sum_{i=1}^{j}\mu_{i}}{|\mu_{j+1}|}.$ (88)

Note that $d_{L}$ is bound from below by $0$. An attractor is considered
chaotic iff the maximal Lyapunov exponent $\mu_{1}$ is positive, implying
$d_{L}>0$.

The Lyapunov spectrum $\Lambda$ is computed numerically by evaluating the mean
logarithmic growth of perturbations in the tangent space of the dynamics. To
do so, we follow a QR-decomposition method, outlined in [29]. The resulting
Lyapunov dimension of activity in our driven reservoirs is show in Fig. 11A
for a range of $g$, $\omega$ and $A$ values. From these maps we extract, for
each value of the target amplitude $A$, the edge-of-chaos connectivity
strength $g^{c}(\omega)$ which corresponds to the minimal value of $g$ for
which $d_{L}>0$ (Fig. 11B).

Figure 11:  Lyapunov dimensionality of driven non-linear networks. A. Lyapunov
dimensionality (color code) of driven activity in non-linear networks as a
function of connectivity strength $g$ and target frequency $\omega$. The
target amplitude is taken to be $A=0.1$ (left), $A=1$ (center) and $A=2$
(right). B. Edge-of-chaos curves extracted from the heat maps in A: we plot
the minimal value of connectivity strength $g^{c}$ for which the Lyapunov
dimensionality is significantly non-zero.