We consider an animal performing a task while both behavior and neural activity are recorded. For example, the task might be to produce a periodic motion, described by the output $z(t)$ of Figure 1A. For simplicity, we assume that the behavioral output can be decoded linearly from the neural activity (Mante et al., 2013; Sadtler et al., 2014; Gallego et al., 2017; Russo et al., 2018; Willett et al., 2021). We can thus write

with readout weights ${𝐰}_{\mathrm{out}}$. The activity of neuron $i\in \{1,\dots ,N\}$ is given by $x_i(t)$, and we refer to the vector $𝐱$ as the state of the network.

Figure 1

Download asset Open asset

Neural activity has to generate the output in some subspace of the state space, where each axis represents the activity of one neuron. In the simplest case (Figure 1B, top), the output is produced along the largest PCs of activity, as shown by the fact that projecting the neural activity $𝐱(t)$ onto the largest PCs returns the target oscillation (Figure 1D, top). We call such dynamics ‘aligned’ because of the alignment between the subspace spanned by the largest PCs and the output vector (red).

There is, however, another possibility. Neural activity may have many other components not directly related to the output and these other components may even dominate the overall activity. In this case (Figure 1B and D, bottom), the two largest PCs are not enough to read out the output, and smaller PCs are needed. We call such dynamics ‘oblique’ because the subspace spanned by the largest PCs and the output vector are poorly aligned.

We consider these two possibilities as distinct dynamical regimes, noting that intermediate situations are also possible. The actual regime of neural dynamics has important consequences for how one interprets neural recordings. For aligned dynamics, analyzing the dynamics within the largest PCs may lead to insights about the computations generating the output (Vyas et al., 2020). For oblique dynamics, such an analysis is hampered by the dissociation between the large PCs and the components generating the output (Russo et al., 2018).

What determines which regime a neural network operates in? Given that behavior is the same, but representations differ, we study this question using trained RNNs. This framework constrains what networks do, but not how they do it (Sussillo, 2014; Barak, 2017). The specific property of representation we are interested in is the alignment, or correlation, between output weights and states:

where the vector norms $\Vert {𝐰}_{\mathrm{out}}\Vert$ and $\Vert 𝐱\Vert$ quantify the magnitude of each vector.

For aligned dynamics, the correlation is large, corresponding to the alignment between the leading PCs of the neural activity and the output weights (Figure 1B, top). In contrast, for oblique dynamics, this correlation is small (Figure 1B, bottom). Note that the concept of correlation can be generalized to accommodate multiple time points and multidimensional output (see section Generalized correlation).

Studying the same task means that the output $z$ is the same, so it is instructive to express it in terms of the correlation $\rho$:

Recent work on feedforward networks showed that $\Vert {𝐰}_{\mathrm{out}}\Vert$ can have a large effect on the resulting representations (Jacot et al., 2018; Chizat et al., 2019). Equation 3 shows that $\rho$ is indeed linked to $\Vert {𝐰}_{\mathrm{out}}\Vert$, but $\Vert 𝐱(t)\Vert$ can also vary. In a detailed analysis in the Methods (sections Analysis of solutions under noiseless conditions to Oblique solutions arise for noisy, nonlinear systems), we show that for recurrent networks, stability considerations preclude $\Vert 𝐱(t)\Vert$ from being small. This implies that if we choose a readout norm $\Vert {𝐰}_{\mathrm{out}}\Vert$ and then train the RNN on a given task, the correlation must compensate.

If we choose small output weights, we expect aligned dynamics, because a large correlation is necessary to generate sufficiently large output (Figure 1B, top). If instead we choose large output weights, we expect oblique dynamics, because only a small correlation keeps the output magnitude from growing too large.

We tested whether output weights can serve as a control knob to select dynamical regimes using an RNN model trained on an abstract version of the cycling task introduced in Russo et al., 2018. The networks were trained to generate a 2D signal that rotated in the plane spanned by two outputs $z_1(t)$ and $z_2(t)$ (Figure 2A). An input pulse at the beginning of each trial indicated the desired direction of rotation. We set up two models with either small or large output weights and trained the recurrent weights of each with gradient descent (section Details on RNN models and training).

Figure 2

Download asset Open asset

After both models learned the task, we projected the network activity into a three-dimensional (3D) space spanned by the two largest PCs of the dynamics $𝐱(t)$. A third direction, ${𝐰}_{\mathrm{out},⟂}$, spanned the remaining part of the first output vector ${𝐰}_{\mathrm{out},1}$. The resulting plots, Figure 2B, corroborate our hypothesis: Small output weights led to aligned dynamics with a large correlation between the largest PCs and the output weights. In contrast, the large output weights of the second network were almost orthogonal, or oblique, to the two leading PCs. Further qualitative differences between the two solutions in terms of the direction of trajectories will be discussed below.

Another way to quantify these regimes is by the ability to reconstruct the output from the large PCs of neural activity, as quantified by the coefficient of determination $R^2$. For the aligned network, the projection on the two largest PCs (Figure 2C, solid) already led to a good reconstruction. For the oblique networks, the two largest PCs were not sufficient. We needed the first eight dimensions (Figure 2C, dashed) to obtain a good reconstruction ($R^2>0.9$). In contrast to the differences in these fits, the neural dynamics themselves were much more similar between the networks. Specifically, 90% of the variance was explained by four and five dimensions for the aligned and oblique networks, respectively.

Can we use the output weights to induce aligned or oblique dynamics in more general settings? We trained RNN models with small or large initial output weights on five different neuroscience tasks (section Task details). All weights (input, recurrent, and output) were trained using the Adam algorithm (section Details on RNN models and training). After training, we measured the three quantities of Equation 3: magnitudes of neural activity and output weights, and the correlation between the two. The results in Figure 3A show that across tasks, initialization with large output weights led to oblique dynamics (small correlation), and with small output weights to aligned dynamics (large correlation). While training could, in principle, change the initially small output weights to large ones (and vice versa), we noticed that this does not happen. Small output weights did increase with training, but the large gap in norms remained. This shows that setting the output weights at initialization can serve to determine their scale after learning under realistic settings. While explaining this observation is beyond the scope of this work, we note that (1) changing the internal weights suffices to solve the task, and (2) the extent to which the output weights change during learning depends on the algorithm and specific parameterization (Jacot et al., 2018; Geiger et al., 2020; Yang and Hu, 2020).

Figure 3

Download asset Open asset

In Figure 3B and C, we adopted the perspective of Figure 2C and quantified how well we can reconstruct the output from a projection of $𝐱$ onto its largest $D$ PCs (section Regression). As expected, both the variance of $𝐱$ explained and the quality of output reconstruction increased, for an increasing number of PCs $D$ (Figure 3B). How both quantities increased, however, differs between the two regimes. While the variance explained increased similarly in both cases, the quality of the reconstruction increased much more slowly for the model with large output weights. We quantified this phenomenon by comparing the dimensions at which either the variance of $𝐱$ explained or $R^2$ reaches 90%, denoted by $D_{x,90}$ and $D_{\mathrm{fit},90}$, respectively.

In Figure 3C, we compare $D_{x,90}$ and $D_{\mathrm{fit},90}$ across multiple networks and tasks. Generally, larger output weights led to larger numbers for both. However, for large output weights, the number of PCs necessary to obtain a good reconstruction increased much more drastically than the dimension of the data. Thus, the output was less well represented by the large PCs of the dynamics for networks with large output weights, in agreement with our notion of oblique dynamics.

Importantly, reaching the aligned and oblique regimes relies on ensuring robust and stable dynamics, which we achieve by adding noise to the dynamics during training. This yields a similar magnitude of neural activity $\Vert 𝐱\Vert$ across networks and tasks (Figure 3A). We show in Methods, section Analysis of solutions under noiseless conditions, that learning in simple, noise-free conditions with large output weights can lead to solutions not captured by either aligned or oblique dynamics; those solutions, however, are unstable. Furthermore, we observed that some of the qualitative differences between aligned and oblique dynamics are less pronounced if we initialized networks with small recurrent weights and initially decaying dynamics (Appendix 1—figure 2).

What are the functional consequences of the two regimes? A hint might be seen in an intriguing qualitative difference between the aligned and oblique solutions for the cycling task in Figure 2. For the aligned network, the two trajectories for the two different contexts (green and purple) are counter-rotating (Figure 2B, top). This agrees with the output, which also counter-rotates as demanded by the task (Figure 2A). In contrast, the neural activity of the oblique network co-rotates in the leading two PCs (Figure 2B, bottom). This is despite the counter-rotating output, since this network also solves the task (not shown). The co-rotation also indicates why reconstructing the output from the leading two PCs is not possible (Figure 2C). Naturally, the dynamics also contain counter-rotating trajectories for producing the correct output, but these are only present in low-variance PCs. (Note also that for aligned networks, one can also observe co-rotation in low-variance PCs, see Appendix 1—figure 3.) Taken together, aligned and oblique dynamics differ in the coupling between leading neural dynamics and output. For aligned dynamics, the two are strongly coupled. For oblique dynamics, the two decouple qualitatively.

Such a decoupling for oblique, but not aligned, dynamics leads to a prediction regarding the universality of solutions (Maheswaranathan et al., 2019; Turner et al., 2021; Pagan et al., 2022). For aligned dynamics, the coupling implies that the internal dynamics are strongly constrained by the task. We thus expect different learners to converge to similar solutions, even if their initial connectivity is random and unstructured. In Figure 4A, we show the dynamics of three randomly initialized aligned networks trained on the cycling task, projected onto the three leading PCs. Apart from global rotations, the dynamics in the three networks are very similar.

Figure 4

Download asset Open asset

For oblique dynamics, the task-defined output exerts weaker constraints on the internal dynamics. Any variability experienced during learning can potentially build up, and eventually create qualitatively different solutions. Three examples of oblique networks solving the cycling tasks indeed show visibly different dynamics (Figure 4B). Further analysis shows that the models also differ in the frequency components in the leading dynamics (Appendix 1—figure 4).

The degree of variability between learners depends on the task. The observable differences in the PC projections were most striking for the cycling task. For the flip-flop task, for example, solutions were generally noisier in the oblique regime than in the aligned but did not have observable qualitative differences in either regime (Appendix 1—figure 5). We quantified the difference between models for the different neuroscience tasks considered before. To compare different neural dynamics, we used a dissimilarity measure invariant under rotation (section Dissimilarity measure) (Williams et al., 2021). The results are shown in Figure 4C. Two observations stand out: First, across tasks, the dissimilarity was higher for networks in the oblique regime than for those in the aligned. Second, both overall dissimilarity and the discrepancy between regimes differed strongly between tasks. The largest dissimilarity (for oblique dynamics) and the largest discrepancy between regimes was found for the cycling. The smallest discrepancy between regimes was found for the flip-flop task. Such a difference between tasks is consistent with the differences in the range of possible solutions for different tasks, as reported in Turner et al., 2021; Maheswaranathan et al., 2019.

What are the underlying mechanisms for the qualitative decoupling in oblique, but not aligned networks? For aligned dynamics, we saw that the small output weights demand large activity to generate the output. In other words, the activity along the largest PCs must be coupled to the output. For oblique dynamics, this constraint is not present, which opens the possibility for small components outside the largest PCs to generate the output. If this is the case, we have a decoupling, such as the observed co-rotation in the cycling task, and the possible variability between solutions. We discuss this point in more detail in the Methods, section Mechanisms behind decoupling of neural dynamics and output.

In the following two sections, we will explore how the decoupling between neural dynamics and output for oblique, but not aligned, dynamics influences the response to perturbations and the effects of noise during learning.

Understanding how networks respond to external perturbations and internal noise requires some insight into how dynamics are generated. Dynamics of trained networks are mostly generated internally, through recurrent interactions. In robust networks, these internally generated dynamics are a prominent part of the largest PCs (among input-driven components; section Analysis of solutions under noiseless conditions and Oblique solutions arise for noisy, nonlinear systems). Internally generated dynamics are sustained by positive feedback loops, through which neurons excite each other. Those loops are low-dimensional, with activity along a few directions of the dynamics being amplified and fed back along the same directions. This results in dynamics being driven by effective feedback loops along the largest PCs (Figure 5A). As shown above, the largest PCs can either be aligned or not aligned, with the output weights. This leads to predictions about how aligned and oblique networks differentiate in their responses to perturbations along different directions.

Figure 5

Download asset Open asset

Our intuition about feedback loops suggests that networks respond strongly to a perturbation that is aligned with the directions contributing to the feedback loop, but weakly to a perturbation that is orthogonal to them. In particular, if a perturbation is applied along the output weights, aligned and oblique dynamics should dissociate, with a strong disruption of dynamics for aligned, but not for oblique dynamics (Figure 5A).

To test this, we compare the response to perturbations along the output direction and the largest PCs. We apply perturbations to the neural activity at a single point in time: $𝐱(t)$ evolves undisturbed until time $t_p$. At that point, it is shifted to $𝐱(t_p)+\mathrm{\Delta }𝐱$. After the perturbation, we let the network evolve freely and compare this evolution to that of an unperturbed copy. Such a perturbation mimics a very short optogenetic perturbation applied to a selected neural population (O’Shea et al., 2022; Finkelstein et al., 2021). In Figure 5B, we show the output after such perturbations for an aligned (top) and an oblique network (bottom) trained on the cycling task. The time point and amplitude are the same for both directions and networks. For each network and type of perturbation, there is an immediate deflection and a long-term response. For both networks, perturbing along the PCs (blue) leads to a long-term phase shift. Only in the aligned network, however, perturbation along the output direction (red) leads to a visible long-term response. In the oblique network, the amplitude of the immediate response is larger, but the long-term response is smaller. Our results for the oblique network, but not for the aligned, agree with simulations of networks generating EMG data from the cycling experiments (Saxena et al., 2022).

To quantify the relative long-term susceptibility of networks to perturbations along output weights or PCs, we sampled from different times $t_p$ and different directions in the 2D subspaces spanned either by the two output vectors or by the two largest PCs. For each perturbation, we measured the loss of the perturbed networks on the original task (excluding the immediate deflection after the perturbation by starting to compute the loss at $t_p+5$). Figure 5C shows that the aligned network is almost equally susceptible to perturbations along the PCs and the output weights. In contrast, the oblique network is much more susceptible to perturbations along the PCs.

We repeated this analysis for oblique and aligned networks trained on the five different tasks. We computed the area under the curve (AUC) for both loss profiles in Figure 5C. We then defined the ‘relative susceptibility’ as the ratio ${\mathrm{AUC}}_{{𝐰}_{\mathrm{out}}}/{\mathrm{AUC}}_{\mathrm{PC}}$, Figure 5D. For aligned networks (red), the relative susceptibility was close to 1 indicating similarly strong responses to both types of perturbations. For oblique networks (yellow), it was much smaller than 1, indicating that long-term responses to perturbations along the output direction were weaker than those to perturbations along the PCs.

In the oblique regime, the output weights are large. To produce the correct output (and not a too large one), the large PCs of the dynamics are almost orthogonal to the output weights. The large output weights, however, pose a robustness problem: Small noise in the direction of the output weights is also amplified at the level of the readout. We show that learning leads to a slow process of sculpting noise statistics to avoid this effect (Figure 11). Specifically, a negative feedback loop is generated that suppresses fluctuations along the output direction (Figure 6A, Figure 10; Kadmon et al., 2020). Because the positive feedback loop that gives rise to the large PCs is mostly orthogonal to the output direction, it remains unaffected by this additional negative feedback loop. A detailed analysis of how learning is affected by noise shows that, for large output weights, the network first learns a solution that is not robust to noise. This solution is then transformed to increasingly stable and oblique dynamics over longer time scales (section Learning with noise for linear RNNs and Oblique solutions arise for noisy, nonlinear systems).

Figure 6

Download asset Open asset

To illustrate the effect of the negative feedback loop, we consider the fluctuations around trial averages. We take a collection of states $𝐱(t)$ and then subtract the task-conditioned averages $\overline{𝐱}(t)$ to compute $\delta 𝐱(t)=𝐱(t)-\overline{𝐱}(t)$. We then project $\delta 𝐱(t)$ onto three different direction categories: the largest PCs of the averaged data $\overline{𝐱}(t)$, the output directions, or randomly drawn directions.

How strongly the activity fluctuates along each direction is quantified by the variance of the projections (Figure 6B). For both aligned and oblique dynamics, the variance is much larger along the PCs than along random directions. This is not necessarily expected, because the PCA was performed on the averaged activity, without the fluctuations. Instead, it is a dynamical effect: the same positive feedback that generates the autonomous dynamics also amplifies the noise (section Oblique solutions arise for noisy, nonlinear systems).

The two network regimes, however, dissociate when considering the variance along the output direction. For aligned dynamics, there is no negative feedback loop, and ${𝐰}_{\mathrm{out}}$ is correlated with the PCs. The variance along the output direction is hence similar to that along the PCs, and larger than along random directions. For oblique dynamics, the negative feedback loop suppresses the fluctuations along the output direction, so that they become weaker than along random directions.

In Figure 6C, we quantify this dissociation across different tasks. We measured the ratio between variance along output and random directions. Aligned networks have a ratio much larger than one, indicating that the fluctuations along the output direction are increased due to the autonomous dynamics along the PCs. In contrast, oblique networks have a ratio smaller than 1 for all tasks, which indicates noise compression along the output.

For the cycling task, we observed that dynamics were qualitatively different for the two regimes, with trajectories either counter- or co-rotating (Figure 2B). Interestingly, the experimental results of Russo et al., 2018, matched the oblique, but not the aligned dynamics. The authors observed co-rotating dynamics in the leading PCs of motor cortex activity despite counter-rotating activity of simultaneously recorded muscle activity. Here, we test whether our theory can help to more clearly and quantitatively distinguish between the two regimes in experimental settings.

In typical experimental settings, we do not have direct access to output weights. We can, however, approximate these by fitting neural data to simultaneously recorded behavioral output, such as hand velocity in a motor control experiment (Figure 7A, top). Following the model above, where the output is a weighted average of the states, we reconstruct the output from the neural activity with linear regression. To quantify the dynamical regime, we then compute the correlation $\rho$ between the weights from fitting and the neural data. Additionally, we can also quantify the alignment by the ‘relative fitting dimension’ $D_{\mathrm{fit},90}/D_{x,90}$, where $D_{\mathrm{fit},90}$ is the number of PCs necessary to recover the output and $D_{x,90}$ number to the number of PCs necessary to represent 90% of the variance of the neural data. We computed both the correlation and the relative fitting dimension for different publicly available data sets (Figure 7B). For details, see section Experimental data.

Figure 7

Download asset Open asset

We started with data sets from two monkeys performing the cycling task (Russo et al., 2018). The data contained motor cortex activity, hand movement, and EMG from the arms, all averaged over multiple trials of the same condition. In Figure 7B, we show results for reconstructing the hand velocity. The correlation was small, $\rho \in [0.04,0.07]$. To obtain a good reconstruction, we needed a substantial fraction of the dimension of the neural data: The relative fitting dimension was $D_{\mathrm{fit},90}/D_{x,90}\in [0.7,0.8]$. Our results agree with previous studies, showing that the best decoding directions are only weakly correlated with the leading PCs of motor cortex activity (Schroeder et al., 2022).

We also analyzed data made available through the Neural Latents Benchmark (NLB) (Pei et al., 2021). In two different tasks, monkeys needed to perform movements along a screen. In a random target task, the monkeys had to point at a randomly generated target position on a screen, with a successive target point generated once the previous one was reached (Makin et al., 2018). In a maze task, the monkeys were trained to follow a trajectory through a maze with their hand (Churchland et al., 2010). In both cases, we reconstructed the finger or hand velocity from neural activity on single trials. The correlation was higher than in the cycling task, $\rho =0.13$. The relative fitting dimension was lower than in the trial-averaged cycling data, albeit still on the same order: $D_{\mathrm{fit},90}/D_{x,90}\in [0.4,0.5]$.

Finally, we considered brain-computer interface (BCI) experiments (Sadtler et al., 2014). In these experiments, monkeys were trained to control a cursor on a screen via activity read out from their motor cortex (Figure 7A, bottom). The output weights generating the cursor velocity were set by the experimenter (so we don’t need to fit). Importantly, the output weights were typically chosen to be spanned by the largest PCs (the ‘neural manifold’), suggesting aligned dynamics. For three example data sets (Golub et al., 2018; Hennig et al., 2018; Degenhart et al., 2020), we obtained higher correlation values, $\rho \in [0.17,0.23]$ (Figure 7B). The relative fitting dimension was much smaller than for the non-BCI data sets, especially for the two largest data sets, where ${\displaystyle D_{\mathrm{f}\mathrm{i}\mathrm{t},90}/D_{x,90}\in [0.03,0.06]}$.

The higher correlation and much smaller relative fitting dimension suggest that, indeed, the neural dynamics arising in BCI experiments are more aligned, and those in non-BCI settings are more oblique. These trends also hold when decoding other behavioral outputs for the cycling task and the NLB tasks (position, acceleration, or EMG), even if the ability to decode and the numerical values for correlation and fitting dimension may fluctuate considerably (Appendix 1—figure 8). Thus, while we do not observe strongly different regimes as in the simulations, we do see an ordering between different data sets according to the alignment between outputs and neural dynamics. It would be interesting to test the differences between BCI and non-BCI data on larger data sets, and different experiments with different dimensions of neural data (Gao et al., 2017; Stringer et al., 2019a).

We consider rate-based RNNs with $N$ neurons. The states $𝐱(t)\in {ℝ}^N$ are governed by

where $W$ is a recurrent weight matrix and $\varphi =\tanh$ a nonlinearity applied element-wise. The network receives a low-dimensional input $𝐬(t)\in {ℝ}^{N_{\mathrm{in}}}$ via input weights $W_{\mathrm{in}}$. It is also driven by white, isotropic noise with zero mean and covariance $𝔼[{\xi }_i(t){\xi }_j(t^{\prime })]=2{\sigma }_{\mathrm{noise}}^2{\delta }_{ij}\delta (t-t^{\prime })$. The initial states $𝐱(0)$ are drawn from a centered normal distribution with variance ${\sigma }_{\mathrm{init}}^2$ at each trial. This serves as additional noise. The output is a low-dimensional, linear projection of the states: $𝐳(t)=W_{\mathrm{out}}𝐱(t)$ with $W_{\mathrm{out}}=[{𝐰}_{\mathrm{out},1},\dots ,{𝐰}_{\mathrm{out},N_{\mathrm{out}}}{]}^T$.

The initial output weights are drawn from centered normal distributions with variance ${\sigma }_{\mathrm{out}}^2/N$. ‘Small’ output weights refer to ${\sigma }_{\mathrm{out}}=1/\sqrt{N}$, and ‘large’ ones to ${\sigma }_{\mathrm{out}}=1$. We have $\Vert {𝐰}_{\mathrm{out},i}\Vert ={\sigma }_{\mathrm{out}}[1+O(1/\sqrt{N})]$ at initialization. Note that large initial output weights are the current default in standard learning environments (Paszke et al., 2017; Yang and Hu, 2020). The recurrent weights were initialized from centered normal distributions with variance $g^2/N$. We chose $g=1.5$ so that dynamics were chaotic before learning (Sompolinsky et al., 1988).

To simulate the noisy RNN dynamics numerically, we used the Euler-Maruyama method (Kloeden and Platen, 1992) with a time step of $\mathrm{\Delta }t$. We used the Adam algorithm (Kingma and Ba, 2014) implemented in PyTorch (Paszke et al., 2017). Apart from the learning rate, we kept the parameters for Adam at the default (some filtering, no weight decay). We selected learning rates and the number of training steps such that learning was relatively smooth and converged sufficiently within the given number of trials. Learning rates were set to $\eta ={\eta }_0/N$. Details for all simulation parameters can be found in Table 1.

For the comparisons over different tasks (Figures 3—6), we trained five networks for each task. All weights ($W_{\mathrm{out}},W,W_{\mathrm{in}}$) were adapted. For the example networks trained on the cycling task (Figures 2, 5, and 6), we used networks with $N=256$ neurons and only changed the recurrent weights $W$. We also trained for longer (5000 training steps) and with a higher learning rate ($\eta =0.1/N$).

The tasks the networks were trained on are taken from the neuroscience literature: a cycling task (Russo et al., 2018), a 3-bit flip-flop task, and a ‘complex sine’ task (with input-dependent frequencies) (Sussillo and Barak, 2013), a context-dependent decision-making task (‘Mante’) (Mante et al., 2013; Schuessler et al., 2020a), and a working memory task comparing the amplitudes of two pulse stimuli (‘Romo’) (Romo et al., 1999; Schuessler et al., 2020a). All tasks have similar structure (Schuessler et al., 2020b): A trial of length $T$ starts with a fixation period (length $t_{\mathrm{fix}}$). This is followed by an input for $t_{\mathrm{stim}}$. For the cycling and flip-flop task, the inputs are pulses of amplitude 1; else see below. After a delay $t_{\mathrm{delay}}$, the output of the network is required to reach an input-dependent value during a time period $t_{\mathrm{dec}}$. During this decision period, we set target points $t_i$ every $\mathrm{\Delta }{t}_{\mathrm{target}}$ time steps. The loss was defined as the mean squared error between network output and target at these time points. Below, we provide further details for each task.

For Figure 3A, we used a generalized correlation measure which allows for multiple output dimensions, multiple time points, and noisy data. Consider neural activity of $N$ neurons at $P$ time points stacked into the matrix $X=[𝐱(1),\dots ,𝐱(P)]\in {ℝ}^{N\times P}$. We assume the states to be centered in time, $\frac{1}{P}\sum _{t=1}^P{x}_i(t)=0$ for $i\in [N]$. The corresponding $D$-dimensional output is summarized in the $D\times P$ matrix

with weights $W_{\mathrm{out}}\in {ℝ}^{N\times D}$. We define the generalized correlation as

The norm is the Frobenius norm, $\Vert X\Vert =\sqrt{\sum _{ij}{X}_{ij}^2}$. In particular, we have

The case of 1D output and a single time step discussed in the main text, Equation 3, is recovered up to the sign, which we discard. Note that in that case, the vectors ${𝐰}_{\mathrm{out}}$ and $𝐱(t)$ should be centered along coordinates to receive a valid correlation. Our numerical results did not change qualitatively when centering across coordinates only or both coordinates and time.

For trajectories with multiple conditions, we stack these instances in a matrix $X\in {ℝ}^{N\times N_c{N}_t}$, with $N_c$ the number of conditions, and $N_t$ the number of time points per trajectory. For noisy trajectories, we first average over multiple instances per condition and time point to obtain a similar matrix $\bar{X}$.

In Figure 3B and C we computed the number of PCs necessary to either represent the dynamics or fit the output. We simulated the trained networks again on their corresponding tasks. We did not apply noise during these simulations, since keeping the same noise as during training would reduce the quality of the output for large output weights; trial averaging yielded similar results to the ones obtained without noise (not shown).

The simulations yielded neural states $X\in {ℝ}^{N\times P}$ and outputs $Z\in {ℝ}^{N_{\mathrm{out}}}\times P$, where $P$ is the number of data points (batch size times number of time points $T$). We applied PCA to the states $X$. The cumulative explained variance ratio obtained from PCA is plotted in Figure 3B. We then projected $X$ onto the first $k$ PCs and fitted these projections to the output with ridge regression (cross-validated, using scikit-learn’s RidgeCV Pedregosa et al., 2011).

For measuring the dissimilarity between learners in Figure 4, we apply a measure following Williams et al., 2021. We define the distance between two sets with $P$ data points $X,Y\in {ℝ}^{N\times P}$ as

where the hat corresponds to centering along the rows, and $Q$ is an orthogonal matrix. The solution to this so-called orthogonal Procrustes’ problem is found via the singular value decomposition $\hat{X}{\hat{Y}}^T=U\mathrm{\Sigma }{V}^T$. The optimal transformation is $Q^{*}=V{U}^T$, and the numerator in Equation 8 is then $\mathrm{Tr}(\hat{X}{\hat{Y}}^T{Q}^{*})=\mathrm{Tr}(\mathrm{\Sigma })$.

Note that this is more restricted than canonical correlation analysis (CCA), which is also commonly used (Gallego et al., 2018; Gallego et al., 2020). In particular, CCA involves whitening the matrices $X$ and $Y$ before applying a rotation (Williams et al., 2021). This sets all singular values to 1. For originally low-D data, this mostly means amplifying the noise, unless the data was previously projected onto a small number of PCs. In the latter case, the procedure still removes the information about how much each PC contributes.

We detail the analyses of neural data in section Different degrees of alignment in experimental settings. We made use of publicly available data sets: data from the cycling task of Russo et al., 2018, two data sets available through the NLB (Pei et al., 2021), and data from monkeys trained on a center-out reaching task with a BCI (Golub et al., 2018; Hennig et al., 2018; Degenhart et al., 2020). For all data sets, we first obtain firing rates $X\in {ℝ}^{N\times T}$, where $N$ is the number of measured neurons, and $T$ the number of data points, (see Figure 7B for these numbers). We also collect the simultaneously measured behavior in the matrix $Z\in {ℝ}^{D\times T}$. In Figure 7, we only analyzed cursor or hand velocity for behavior, so that the output dimension is $D=2$. See Appendix 1—figure 8 for similar results for hand position and acceleration or the largest two PCs of the EMG data recorded for the cycling task.

For the cycling task, the firing rates were binned in 1 ms bins and convolved with a 25 ms Gaussian filter. The mean firing rate was 22 and 18 Hz for the two monkeys, respectively. For the NLB data, spikes came in 1 ms bins. We binned data to 45 ms bins and applied a Gaussian filter with 45 ms width. This increased the quality of the fit, as firing rates were much lower (mean of 5 Hz for both) than in the cycling data set. For the BCI experiments, firing rates came as spike counts in 45 ms bins. The mean firing rate was (45, 45, 55) Hz for the data of Golub et al., 2018; Hennig et al., 2018; Degenhart et al., 2020, respectively. In agreement with the original BCI experiments, we did not apply a filter to the neural data.

For fitting, we centered both firing rates $X$ and output $Z$ across time (but not coordinates). We also added a delay of 100 ms between firing rates and output for the cycling and NLB data sets, which increased the quality of the fits. We then fitted the output $Z$ to the firing rates $X$ with ridge regression, with regularization obtained from cross-validation. We treated the coefficients as output weights $W_{\mathrm{out}}$. The trial average data of the cycling tasks was very well fitted for both monkeys, $R^2=[0.97,0.98]$. For the NLB tasks with single-trial data, the fits were not as good, $R^2=[0.73,0.69]$. For two of the BCI data sets (Golub et al., 2018; Hennig et al., 2018), the output weights were also given, and we checked that the fit recovers these. For the third BCI data set (Degenhart et al., 2020), we did not have access to the output weights, and only access to the cursor velocity after Kalman filtering. Here, fitting yielded $R^2=0.83$.

For the fitting dimension $D_{\mathrm{f}\mathrm{i}\mathrm{t},90}/D_{x,90}$ in Figure 7B, bottom, we used an adapted definition of $D_{\mathrm{f}\mathrm{i}\mathrm{t},90}$: Because $R^2=90$% is not reached for all data sets, we asked for the number of PCs necessary to obtain 90% of the $R^2$ value obtained for the full data set.

We also considered whether the correlation $\rho$ scales with the number of neurons $N$. In our model, oblique and aligned dynamics can be defined in terms of such a scaling: aligned dynamics have highly correlated output weights and low-dimensional dynamics, so that $\rho \sim N^0=1$, i.e., independent of the network size. For oblique dynamics, large output weights with norm ${𝐰}_{\mathrm{out}}\sim 1$ lead to vanishing correlation, $\rho \sim 1/\sqrt{N}$. This is indeed similar to the relation between two random vectors, for which the correlation is precisely $1/\sqrt{N}$ (in the limit of large $N$). In Figure 8, we show the scaling of $R^2$ and $\rho$ with the number of subsampled neurons. For the cycling task and NLB data, the correlation scaled slightly weaker than $\rho \sim N^{-1/2}$. For the BCI data, the scaling was closer to $\rho \sim N^{-1/4}$ which is in between the aligned and oblique regimes of the model. These insights, however, are limited due to the trial averaging for the cycling task and the limited number of time points for the NLB tasks (not enough to reach $R^2=1$). Applying these measures to larger data sets could yield more definitive insights.

Figure 8

Download asset Open asset

In the sections below, we explore in detail under which conditions aligned and oblique solutions arise, and which other solutions arise if these conditions are not met.

We first consider small output weights and show that these lead to aligned solutions. Then, for large output weights, we show that without noise, two different, unstable solutions arise. Finally, we consider how adding noise affects learning dynamics. For a linear model, we can solve the dynamics of learning analytically and show how a negative feedback loop arises, that suppresses noise along the output direction. However, the linear model does not yield an oblique solution, so we also consider a nonlinear model for which we show in detail why oblique solutions arise.

We start by analyzing a simplified version of the network dynamics (Equation 4): autonomous dynamics without noise,

with fixed initial condition $𝐱(0)$. We assume a 1D output $z(t)$ and a target $\hat{z}(t_i)$ defined on a finite set of time points $t_i$.

We illustrate the theory with an example of a simple sine wave task (Figure 9). We demand the network to autonomously produce a sine wave with fixed frequency $f=0.1$. At the beginning of the task, the network receives an input pulse that sets the starting point of the trajectory. We set the noise ${\sigma }_{\mathrm{init}}$ on the initial state $𝐱(0)$ to zero. We define the target as 20 target points in the interval $t\in [1,21]$ (two cycles; purple dots in Figure 9).

Figure 9

Download asset Open asset

In the next two sections, we aim to understand how adding noise affects dynamics and training. We start with a simple setting of a linear RNN which allows us to track the learning dynamics analytically. Despite its simplicity, this setting already captures a range of observations: different time scales for learning the bias and variance part, and the rise of a negative feedback loop for noise suppression. Oblique dynamics, however, do not arise, showing that these need autonomously generated, nonlinear dynamics, covered in section Oblique solutions arise for noisy, nonlinear systems.

We consider a linear network driven by a constant input and additional white noise. The dynamics read

with noise $𝝃$ as in Equation 4. We focus on a simplified task, which is to produce a constant nonzero output $\hat{z}$ once the average dynamics converged, i.e., for large trial times. The output is $z={𝐰}_{\mathrm{out}}^T𝐱=\bar{z}+\delta z$, where the bar denotes average over the noise, and the delta fluctuations around the average. The average is given by $\bar{z}={𝐰}_{\mathrm{out}}^T(I-W{)}^{-1}{𝐰}_{\mathrm{in}}$. We train the network by changing only the recurrent weights $W$ via gradient descent. For small output weights, the fluctuations are too small to affect training: $\delta z=O(1/\sqrt{N})$. Apart from a small correction, learning dynamics are then the same as for small output weights and no noise, a setting analyzed in Schuessler et al., 2020b. Here, we only consider the case of large output weights.

The loss separates into two parts, $L=L_{\mathrm{bias}}+L_{\mathrm{var}}$ with $L_{\mathrm{bias}}=(\bar{z}-\hat{z}{)}^2$ and $L_{\mathrm{var}}=\overline{\delta z^2}=\mathrm{var}(\delta z)$. Learning aims to minimize the sum. We first consider learning based on each part alone and then join both to describe the full learning dynamics.

Learning based on the bias part alone converges to a lazy solution (see section Details linear model: Bias only: lazy learning). For no initial weights, $W_0=0$, we have to leading order

with

Thus, $\mathrm{\Delta }W$ is rank one with norm $\Vert \mathrm{\Delta }W{\Vert }_2\sim \frac{1}{\sqrt{N}}$. Furthermore, for a learning rate $\eta$, it converges in $O(\frac{1}{\eta N})$ time steps. We will see that this is very fast compared to learning the variance part.

Optimizing mean and fluctuations occurs on different time scales

We now consider learning both the mean and the variance part. For zero initial recurrent weights, $W_0=0$, the input and output vectors make up the only relevant directions in space. We thus express the recurrent weights as a rank-two matrix, $W=\hat{U}M{\hat{U}}^T$, with orthonormal basis $\hat{U}=[{\widehat{𝐰}}_{\mathrm{out}},{\widehat{𝐰}}_{\mathrm{in},⟂}]$. The hats indicate normalized vectors. For large networks and large output weights, the first vector is already normalized, ${\widehat{𝐰}}_{\mathrm{out}}={𝐰}_{\mathrm{out}}$. The second vector is the input weights after Gram-Schmidt. Assuming that ${𝐰}_{\mathrm{out}}$ and ${𝐰}_{\mathrm{in}}$ are drawn independently, we have a small, random correlation ${\rho }_{io}=O(1/\sqrt{N})$, and can write ${\widehat{𝐰}}_{\mathrm{in},⟂}={\widehat{𝐰}}_{\mathrm{in}}-{𝐰}_{\mathrm{out}}{\rho }_{io}+O(\frac{1}{N}).$

We computed the learning dynamics in terms of the coefficient matrix $M$, using the same tools introduced above, and the insight that learning the mean is much faster than learning to reduce the variance. The details are relegated to section Details linear model: Bias and variance combined, here, we summarize the results. We obtained

(30) $\begin{array}{r}M(\tau )=\left[\begin{array}{cc}{\lambda }_{-}(\tau )& \frac{b(\tau )}{\sqrt{N}}\\ 0& 0\end{array}\right].\end{array}$

Before discussing the temporal evolution of the two components, we analyze the structure of the matrix. The eigenvalue $M_{11}={\lambda }_{-}$ is a negative feedback loop along the eigenvector ${𝐰}_{\mathrm{out}}$, and $M_{21}=b/\sqrt{N}$ is a small feedforward component $b$ which maps the input to the output. Along the input direction ${\widehat{𝐰}}_{\mathrm{in}}$, learning does not change the dynamics: The second eigenvalue, corresponding to this direction, is zero.

The dynamics unfold on two time scales. First, there is a very fast learning of the bias via the feedforward coefficient

(31) $b(\tau )=(1-e^{-2N\eta \tau })(\hat{z}-\sqrt{N}{\rho }_{io}).$

During this phase, the eigenvalue $M_{11}={\lambda }_{-}$ remains at zero, so that overall weight changes remain small, $\Vert W\Vert \sim 1/\sqrt{N}$. The average fixed point also does not change much, $\Vert \mathrm{\Delta }\overline{𝐱}(\tau )\Vert =b_1(\tau )=O(1),$ in comparison to the fixed point before learning, $\Vert {\overline{𝐱}}_0\Vert =\Vert {𝐰}_{\mathrm{in}}\Vert =\sqrt{N}.$ The loss evolves as

(32) $L_{\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}}=e^{-4N\eta \tau }(\hat{z}-\sqrt{N}{\rho }_{io}{)}^2.$

The variance part of the loss does not change during this phase.

In a second, slower learning phase, the eigenvalue ${\lambda }_{-}$ evolves like above, Equation 28, the case where only the variance part is learned. The second coefficient compensates for the resulting change in the output. This compensation happens at a much faster time scale ($N^3$ times faster than ${\lambda }_{-}$), so we consider its steady state:

(33) $b(\tau )=\hat{z}[1-{\lambda }_{-}(\tau )]-\sqrt{N}{\rho }_{io}.$

Because of this compensation, the bias part of the loss always remains at zero, and the full loss $L(\tau )=L_{\mathrm{var}}(\tau )$ evolves as before, Equation 29. Meanwhile, the average fixed point does not change anymore; we have $\Vert \mathrm{\Delta }\overline{𝐱}\Vert =\hat{z}-{\rho }_{io}\sqrt{N}$.

We compare our theoretical predictions against numerical simulations in Figure 10. For the task, we let the linear network dynamics converge from $𝐱(0)=\mathrm{𝟎}$ until $t=15$, and demand the output $z(t)$ to be at the target $\hat{z}$ during the interval $t\in [15,20]$. Because the first learning phase converges $N$ times faster than the second one (with $N=256$), using a single learning rate $\eta$ is problematic. One can either observe the first phase only (for small $\eta$) or risk unstable learning during the first phase (for large $\eta$). We thus split learning into two parts with adapted learning rates. For the initial phase, we set a learning rate to $\eta ={\eta }_0/N$, with ${\eta }_0=0.002$ (Figure 10, left column). For the second phase, we set $\eta ={\eta }_0$ (Figure 10, right column). Theory and simulation agree well for both phases. Small deviations can be observed for the second phase and long learning times: nonzero coefficients $M_{21}$ and $M_{22}$ and a corresponding increase in the norm $\Vert \mathrm{\Delta }\overline{𝐱}\Vert$. Testing with larger network sizes showed that these errors decreased as $O(1/N)$, which is consistent with our theory above (not shown).

Figure 10

Download asset Open asset

Our results show that learning in this linear system is a hybrid between oblique and aligned during the second phase. We have a large term that compresses the output noise, but a very small, ‘lazy’ correction in the feedforward component that corrects the output, and the states are only marginally changed. Although this system is a somewhat degenerate limiting case, we can still derive important insights. The two most striking features – the different time scales of the two learning processes, and the slow emergence of negative feedback, ${\lambda }_{-}(\tau )\sim -{\tau }^{1/3}$, along the output – are found robustly also for nonlinear networks and other tasks.

We now examine the origin of oblique solutions. The linear system did not yield oblique solutions, so we turn to a nonlinear model. We consider a 1D flip-flop task, where the network has to yield a constant, nonzero output $\hat{z}$ depending on the sign of the last input pulse. We further simplify the analysis by only considering the steady state of a network, not how the input pulse mediates the transition. At the output, we thus only consider the average $\bar{z}$ and the fluctuations $\delta z$. As for the linear network, the loss splits into a bias part $L_{\mathrm{bias}}=(\bar{z}-\hat{z}{)}^2$ and variance part $L_{\mathrm{var}}=\overline{\delta z^2}$. As before, we assume that learning only acts on a low-dimensional parameter matrix $M$.

Because following the learning dynamics in nonlinear networks is difficult, we take a different approach. We develop a mean field theory to show how noise affects the dynamics of a nonlinear network with autonomous fixed points. This allows us to compute the loss components $L_{\mathrm{bias}}$ and $L_{\mathrm{var}}$ in terms of $M$. We then show that the minimum of the loss function corresponds to oblique solutions. This leads to a clear interpretation of the mechanisms pushing for oblique solutions. Finally, we show that the theory quantitatively predicts the outcome of learning with gradient descent.

Rank-two connectivity model with fixed point

We first introduce the connectivity model and compute the latent dynamics using mean field theory. We constrain the recurrent connectivity to a rank-two model of the form $W=\frac{1}{N}UM{U}^T$, where $M$ is a 2×2 coefficient matrix to be learned. The randomly drawn projection matrix $U\in {ℝ}^{N\times 2}$ has entries $U_{ia}$ drawn independently from a standard normal distribution. This implies orthogonality to leading order, $\frac{1}{N}{U}^{T}U=I_2+O(1/\sqrt{N})$. We will discard the correction term, as it does not change our results apart from a constant bias. We further assume that the input and output vectors are spanned by $U$, although not necessarily identified with the components of $U$ as in the previous section. Note also that for clarity we do not normalize $U$ as $\hat{U}$ before. The assumption of Gaussian connectivity greatly simplifies the math, while its restrictions are irrelevant to the task considered here (Schuessler et al., 2020a; Beiran et al., 2021; Dubreuil et al., 2021).

To understand the dynamics Equation 4 analytically, we make use of the low-rank connectivity. Following previous work (Rivkind and Barak, 2017; Kadmon et al., 2020), we split the dynamics into two parts: a parallel part ${𝐱}_{\parallel }$ in the subspace spanned by $U$, and an orthogonal part ${𝐱}_{⟂}$. This yields

(34) ${\dot{\mathbf{x}}}_{\parallel }=-{\mathbf{x}}_{\parallel }+\frac{1}{N}UM{U}^T\varphi ({\mathbf{x}}_{\parallel }+{\mathbf{x}}_{\perp })+\frac{1}{N}U{U}^T\mathit{\xi },$

(35) ${\dot{\mathbf{x}}}_{\perp }=-{\mathbf{x}}_{\perp }+(I-\frac{1}{N}U{U}^T)\mathit{\xi }.$

Notice that the parallel part is partially driven by the orthogonal one, but not vice versa. The parallel part can be expressed in terms of the latent variable

(36) $\mathit{\kappa }=\frac{1}{N}{U}^T\mathbf{x}=\frac{1}{N}{U}^T{\mathbf{x}}_{\parallel }.$

The scaling here ensures that $𝜿$ is order one if ${𝐱}_{\parallel }$ has order-one states. Because the readout is assumed to be spanned by $U$, the output is fully determined by the parallel part. We can write

(37) $z={\mathbf{w}}_{\mathrm{o}\mathrm{u}\mathrm{t}}^T{\mathbf{x}}_{\parallel }=\sqrt{N}{\mathbf{v}}_{\mathrm{o}\mathrm{u}\mathrm{t}}^T\mathit{\kappa },$

with projected output weights ${𝐯}_{\mathrm{out}}=\frac{1}{\sqrt{N}}{U}^T{𝐰}_{\mathrm{out}}$. Note that we assume large output weights, $\Vert {𝐰}_{\mathrm{out}}\Vert =1$, so that ${𝐯}_{\mathrm{out}}$ is also normalized.

We split the latent state into its average over the noise $𝝃$ and fluctuations, $𝜿=\overline{𝜿}+\delta 𝜿$. Similarly, the output splits into $z=\bar{z}+\delta z.$ The loss then has two components, $L=L_{\mathrm{bias}}+L_{\mathrm{var}}$, with

(38) $L_{\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}}=(\overline{z}-\hat{z}{)}^2={\left(\sqrt{N}{\mathbf{v}}_{\mathrm{o}\mathrm{u}\mathrm{t}}^T\overline{\mathit{\kappa }}-\hat{z}\right)}^2,$

(39) $L_{\mathrm{v}\mathrm{a}\mathrm{r}}=\overline{\delta z^2}=N{\mathbf{v}}_{\mathrm{o}\mathrm{u}\mathrm{t}}^T\overline{\delta \mathit{\kappa }\delta {\mathit{\kappa }}^T}{\mathbf{v}}_{\mathrm{o}\mathrm{u}\mathrm{t}}.$

We want to understand situations with small loss. For the bias term, the average $\overline{𝜿}$ must be either small or oblique to the readout weights. For the variance term, the covariance of the fluctuations, $\mathrm{cov}(\delta 𝜿)=\overline{\delta 𝜿\delta {𝜿}^T}$, must be compressed along the output direction: Even though the covariance is $O(1/N)$, it still projects to the output at $O(1)$, as reflected by the factor $N$ in Equation 39. Figure 11 illustrates the situation in a cartoon from this high-level perspective.

Figure 11

Download asset Open asset

To understand the underlying mechanisms in detail, we next explore how the relevant variables $\overline{𝜿}$ and $\delta 𝜿$ are determined by the coupling matrix $M$. We do so by applying mean field theory, following previous works (Mastrogiuseppe and Ostojic, 2018; Schuessler et al., 2020a; Kadmon et al., 2020; Schuecker et al., 2018). Detailed derivations can be found in the section Details nonlinear autonomous system with noise. Here, we present the high-level results. The average dynamics converge to a fixed point determined by the equation for the latent variable,

(40) $\overline{\mathit{\kappa }}=⟨{\varphi }^{\prime }⟩M\overline{\mathit{\kappa }}+O(\frac{1}{\sqrt{N}}).$

The $O(1/\sqrt{N})$ term is a constant offset for any given network that we ignore without loss of generality. The average slope is

(41) $⟨{\varphi }^{\prime }⟩=⟨{\varphi }^{\prime }(\overline{{\sigma }_x}u)⟩=\int \mathcal{D}u{\varphi }^{\prime }\left(\overline{{\sigma }_x}u\right),$

with variance

(42) $\overline{{\sigma }_x^2}=\Vert \overline{\mathit{\kappa }}{\Vert }^2+\overline{{\sigma }_{\perp }^2}.$

(We have $\overline{{\sigma }_x}=\sqrt{\overline{{\sigma }_x^2}}$ because the fluctuations are small.) The orthogonal variance is simply the variance of the noise, ${\displaystyle \overline{{\sigma }_{\perp }^2}={\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^2}$. The fixed point (Equation 40) implies that for a nonzero fixed point, the matrix $M$ must have an eigenvalue

(43) ${\lambda }_{+}=\frac{1}{⟨{\varphi }^{\prime }⟩}.$

This is very similar to the noiseless situation discussed briefly above, Equation 14. However, here the additional variance $\overline{{\sigma }_{⟂}^2}$ decreases the average slope due to saturation of the nonlinearity. This in turn increases the minimal eigenvalue for a nonzero fixed point, which can be found by setting $\overline{𝜿}=\mathrm{𝟎}$, and hence $\overline{{\sigma }_x}=\overline{{\sigma }_{⟂}}$:

(44) ${\lambda }_{+,\mathrm{m}\mathrm{i}\mathrm{n}}=\frac{1}{⟨{\varphi }^{\prime }(\overline{{\sigma }_{\perp }}u)⟩}.$

In other words, the noise decreases the effective gain $⟨{\varphi }^{\prime }⟩$, and thus the connectivity eigenvalue ${\lambda }_{+}$ needs to compensate. From the point of view of the spectrum, we are thus pushed away from the margin $1+ϵ$. These considerations, however, do not exclude the possibility that dynamics converge to an average fixed point that is small and correlated, which is at odds with oblique dynamics. To understand why learning leads to oblique dynamics, we need to move beyond the average $\overline{𝜿}$ and take into account the fluctuations $\delta 𝜿$.

Minimizing the loss by balancing saturation and negative feedback loop

To gain an understanding of how the output fluctuations responsible for $L_{\mathrm{bias}}$ can be reduced, we first consider the case of a symmetric coefficient matrix $M$. Simulations of networks trained with gradient descent below show that this approximation is reasonable. For symmetric $M$, the orthogonal eigenvectors ${𝐯}_{\pm }$ with eigenvalues ${\lambda }_{\pm }$ of the recurrent weights $M$ are also eigenvectors of $A$, in that case corresponding to the eigenvalues ${\gamma }_{\pm }$. This allows to diagonalize the Lyapunov Equation 48 and yields the solution

(49) ${\mathrm{\Sigma }}_A=\frac{{\mathbf{v}}_{+}{\mathbf{v}}_{+}^T}{-{\gamma }_{+}}+\frac{{\mathbf{v}}_{-}{\mathbf{v}}_{-}^T}{-{\gamma }_{-}}.$

For the loss (Equation 39), we further need the relation between the eigenvectors ${𝐯}_{\pm }$ and the output weights. Because the fixed point $\overline{𝜿}$ is parallel to the eigenvector, we have ${𝐯}_{\mathrm{out}}^T{𝐯}_{+}=\rho$. For the other eigenvector, orthogonality yields ${𝐯}_{\mathrm{out}}^T{𝐯}_{-}=\sqrt{1-{\rho }^2}$. Inserting this into Equations 39 and 47 yields an expression in terms of the Jacobian eigenvalues ${\gamma }_{\pm }$, the correlation $\rho$, and the norm of the fixed point $\Vert \overline{𝜿}\Vert$:

(50) $L_{\mathrm{v}\mathrm{a}\mathrm{r}}={\rho }^2\left[\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^2}{-{\gamma }_{+}}+\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^4}{\Vert \overline{\mathit{\kappa }}{\Vert }^2}\frac{{\gamma }_{+}}{2({\gamma }_{+}-2)}\right]+(1-{\rho }^2)\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^2}{-{\gamma }_{-}}.$

For more explicit insight, we choose the nonlinearity $\varphi (x)=\mathrm{erf}(\alpha x)$ with $\alpha =\sqrt{\pi }/2$. Similar to $\tanh$, this function is bounded between ±1 and has slope ${\varphi }^{\prime }(0)=1$ at the origin. For this function, we can explicitly compute the relevant Gaussian integrals. The minimal eigenvalue of $M$ to produce a fixed point, Equation 44, is then given by ${\lambda }_{+,\min }=1+2{\alpha }^2{\sigma }_{\mathrm{noise}}^2$. For ${\lambda }_{+}>{\lambda }_{+,\min }$, the resulting fixed point has norm

(51) $\Vert \overline{\mathit{\kappa }}{\Vert }^2=\frac{{\lambda }_{+}^2-{\lambda }_{+,\mathrm{m}\mathrm{i}\mathrm{n}}^2}{2{\alpha }^2},$

as shown in Figure 12A. The larger eigenvalue of the Jacobian is ${\gamma }_{+}=-({\lambda }_{+}^2-{\lambda }_{+,\min }^2)/{\lambda }_{+}^2$. We next obtain the correlation between fixed point and output weights by assuming that the bias part of the loss (Equation 38) is kept at zero. This is reasonable because it requires only a small adaptation to the weights that leaves the variance part mostly untouched. The resulting correlation is

(52) ${\rho }^2=\frac{1}{N}\frac{{\hat{z}}^2}{\Vert \overline{\mathit{\kappa }}{\Vert }^2},$

Figure 12

Download asset Open asset

so that increasing ${\lambda }_{+}$ also decreases the correlation (Figure 12B). Finally, we obtain an expression for the variance part of the loss only in terms of the eigenvalues of $M$:

(53) $L_{\mathrm{v}\mathrm{a}\mathrm{r}}={\rho }^2\left(\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^2{\lambda }_{+}^2}{{\lambda }_{+}^2-{\lambda }_{+,\mathrm{m}\mathrm{i}\mathrm{n}}^2}+\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^4{\alpha }^2}{3{\lambda }_{+}^2-{\lambda }_{+,\mathrm{m}\mathrm{i}\mathrm{n}}^2}\right)+(1-{\rho }^2)\frac{{\sigma }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^2}{1-\frac{{\lambda }_{-}}{{\lambda }_{+}}}.$

We show $L_{\mathrm{var}}$ over ${\lambda }_{+}$ for different negative feedback loop sizes ${\lambda }_{-}$ (Figure 12C) and different network sizes $N$ (Figure 12D). The first term diverges at the phase transition where the fixed point appears, ${\lambda }_{+}\searrow {\lambda }_{+,\min }$. Learning will thus push the weights away from the phase transition toward larger ${\lambda }_{+}$. With such increasing ${\lambda }_{+}$, the fixed point norm increases, and the fixed point rotates away from the output, decreasing the correlation. This in term emphasizes the last term, scaled by $1-{\rho }^2$. The last term is reduced with increasingly negative ${\lambda }_{-}$, corresponding to the negative feedback loop that suppresses noise.

Learning can in principle strengthen this feedback further and further, ${\lambda }_{-}\to -\infty$ (apart from possible stability issues; Kadmon et al., 2020). However, as for the linear network, section Learning with noise for linear RNNs, this process takes time. We thus assume ${\lambda }_{-}$ to be fixed and search for a minimum across ${\lambda }_{+}$. For ${\lambda }_{-}<0$, the last term in Equation 53 increases with increasing ${\lambda }_{+}$: the effective feedback loop in the full, nonlinear system is weakened by saturation. The loss $L_{\mathrm{var}}$ thus has a minimum at some moderate ${\lambda }_{+}$.

To illustrate the mechanisms described above, we simulated networks at different ${\lambda }_{+}$ and with ${\lambda }_{-}=-5$. For each ${\lambda }_{+}$, we compute $\rho$ according to Equation 52. Setting ${𝐯}_{\mathrm{out}}=[1,0{]}^T$, we then set the resulting symmetric $M=V\mathrm{\Lambda }{V}^T$. For each network sample, we then draw independent random projections $U\in {ℝ}^{N\times 2}$. We started simulations at either one of the two nonzero fixed points. For ${\lambda }_{+}$ just above ${\lambda }_{+,\min }$, the noise pushes activity from one basin of attraction to the next (Figure 12C). The resulting output becomes centered around zero and independent of the initial condition for long simulation times (Figure 12F). For the optimal ${\lambda }_{+}$, the trajectories remain close to either one fixed point (Figure 12F). The output forms two overlapping distributions, each closely matching the target on average (Figure 12I). For larger ${\lambda }_{+}$, the fixed points become increasingly larger (Figure 12G). While this decreases the probability of leaving the basin of attraction even further, the variance along the output weights becomes larger (slightly wider histograms in Figure 12J). Note that the mean starts to deviate from the prediction. This is not covered by our theory and is potentially due to the linearization of the fluctuations.

All-in-all, this section revealed a potential path to oblique solutions, initiated by the large fluctuations close to a phase transition, and the interplay between the negative feedback loop and saturation. In the following section, we show that learning via gradient descent actually follows this path and that the parameters predicted by the minimum of $L_{\mathrm{var}}$ quantitatively predict solutions from learning.