Abstract
We present a framework for compactly summarizing many recent results inefficient and/or biologically plausible online training of recurrent neuralnetworks (RNN). The framework organizes algorithms according to severalcriteria: (a) past vs. future facing, (b) tensor structure, (c) stochastic vs.deterministic, and (d) closed form vs. numerical. These axes reveal latentconceptual connections among several recent advances in online learning.Furthermore, we provide novel mathematical intuitions for their degree ofsuccess. Testing various algorithms on two synthetic tasks shows thatperformances cluster according to our criteria. Although a similar clusteringis also observed for gradient alignment, alignment with exact methods does notalone explain ultimate performance, especially for stochastic algorithms. Thissuggests the need for better comparison metrics.
Quick Read (beta)
A Unified Framework of Online Learning Algorithms for Training Recurrent Neural Networks
Abstract
We present a framework for compactly summarizing many recent results in efficient and/or biologically plausible online training of recurrent neural networks (RNN). The framework organizes algorithms according to several criteria: (a) past vs. future facing, (b) tensor structure, (c) stochastic vs. deterministic, and (d) closed form vs. numerical. These axes reveal latent conceptual connections among several recent advances in online learning. Furthermore, we provide novel mathematical intuitions for their degree of success. Testing various algorithms on two synthetic tasks shows that performances cluster according to our criteria. Although a similar clustering is also observed for gradient alignment, alignment with exact methods does not alone explain ultimate performance, especially for stochastic algorithms. This suggests the need for better comparison metrics.
Center for Neural Science
New York University
New York, NY 10006, USA Kyunghyun Cho [email protected]
New York University
Facebook AI Research
CIFAR Azrieli Global Scholar Cristina Savin [email protected]
Center for Neural Science
Center for Data Science
New York University
Editor:
Keywords: realtime recurrent learning, backpropagation through time, approximation, biologically plausible learning, local, online
1 Introduction
Training recurrent neural networks (RNN) to learn sequence data is traditionally done with stochastic gradient descent (SGD), using the backpropagation through time algorithm (BPTT, Werbos et al., 1990) to calculate the gradient. This requires “unrolling” the network over some range of time steps $T$ and performing backpropagation as though the network were feedforward under the constraint of sharing parameters across time steps (“layers”). BPTT’s success in a wide range of applications (Mikolov et al., 2010; Graves, 2013; Bahdanau et al., 2016, 2014; Cho et al., 2015; Graves et al., 2016) has made it the industry standard; however, there exist alternative online algorithms for training RNNs. These compute gradients in real time as the network runs forward, without explicitly referencing past activity or averaging over batches of data. There are two reasons for considering online alternatives to BPTT. One is practical: computational costs do not scale with $T$. The other is conceptual: human brains are able to learn longterm dependencies without explicitly memorizing all past brain states, and understanding online learning is a key step in the larger project of understanding human learning.
The classic online learning algorithm is realtime recurrent learning (RTRL, Williams and Zipser, 1989), which is equivalent to BPTT in the limit of a small learning rate (Murray, 2019). RTRL recursively updates the total derivative of the hidden state with respect to the parameters, eliminating the need to reference past activity but introducing an order $n\text{(hidden units)}\times {n}^{2}\text{(parameters)}={n}^{3}$ memory requirement. In practice, this is often more computationally demanding than BPTT (order $nT$), hence not frequently used in applications. Nor is RTRL at face value a good model of biological learning, for the same reason: no known biological mechanism exists to store—let alone manipulate—a float for each synapseneuron pair. Thus RTRL and online learning more broadly have remained relatively obscure footnotes to both the deep learning revolution itself and its impact on computational neuroscience.
Recent advances in recurrent network architectures have brought the issue of online learning back into the spotlight. While vanishing gradients used to significantly limit the extent of the temporal dependencies that an RNN could learn, new architectures like LSTMs (Hochreiter and Schmidhuber, 1997) and GRUs (Cho et al., 2014) have dramatically expanded this learnable time horizon. Unfortunately, taking advantage of this capacity requires an equally dramatic expansion in computational resources, if using BPTT. This has led to an explosion of novel online learning algorithms (Tallec and Ollivier, 2017; Mujika et al., 2018; Roth et al., 2019; Murray, 2019; Jaderberg et al., 2017) which aim to improve on the efficiency of RTRL, in many cases using update rules that might be implemented by a biological circuit.
The sheer number and variety of these approaches pose challenges for both theory and practice. It is not always completely clear what makes various algorithms different from one another, how they are conceptually related, or even why they might work in the first place. There is a pressing need in the field for a cohesive framework for describing and comparing online methods. Here we aim to provide a thorough overview of modern online algorithms for training RNNs, in a way that provides a clearer understanding of the mathematical structure underlying all these different approaches. Our framework organizes the existing literature along several axes that encode meaningful conceptual distinctions:

a)
Past facing vs. future facing

b)
The tensor structure of the algorithm

c)
Stochastic vs. deterministic update

d)
Closed form vs. numerical solution for update
These axes will be explained in detail later, but briefly: the past vs. future axis is a root distinction that divides algorithms by the type of gradient they calculate, while the other three describe their representations and update principles. Table 1 contains (to our knowledge) all recently published online learning algorithms for RNNs, categorized according to these criteria. We can already see that many combinations of these characteristics manifest in the literature, suggesting that new algorithms could be developed by mixing and matching properties. (We provide a concrete example of this in §3.4.)
Algorithm  Facing  Tensor  Update  Memory  Time  

RTRL  Past  ${M}_{kij}$  Deterministic  Closedform  ${n}^{3}$  ${n}^{4}$ 
UORO  Past  ${A}_{k}{B}_{ij}$  Stochastic  Closedform  ${n}^{2}$  ${n}^{2}$ 
KFRTRL  Past  ${A}_{j}{B}_{ki}$  Stochastic  Closedform  ${n}^{2}$  ${n}^{3}$ 
RKFRTRL  Past  ${A}_{i}{B}_{jk}$  Stochastic  Closedform  ${n}^{2}$  ${n}^{3}$ 
$r$OK  Past  ${\sum}_{l=1}^{r}{A}_{lj}{B}_{lki}$  Stochastic  Numerical  $r{n}^{2}$  $r{n}^{3}$ 
KeRNL  Past  ${A}_{ki}{B}_{ij}$  Deterministic  Numerical  ${n}^{2}$  ${n}^{2}$ 
RFLO  Past  ${\delta}_{ki}{B}_{ij}$  Deterministic  Closedform  ${n}^{2}$  ${n}^{2}$ 
EBPTT  –  –  Deterministic  Closedform  $nT$  ${n}^{2}$ 
FBPTT  Future  –  Deterministic  Closedform  $nT$  ${n}^{2}T$ 
DNI  Future  ${A}_{li}$  Deterministic  Numerical  ${n}^{2}$  ${n}^{2}$ 
Here we describe each algorithm in unified notation that makes clear their classification by these criteria. In the process, we generate novel intuitions about why different approximations can be successful and discuss some of the finer points of their biological plausibility. Finally, we simulate each algorithm on a common set of synthetic tasks with vanilla RNN architecture for simplicity. We compare performance and analyze gradient alignments to see to what extent their empirical similarity is predicted by their similarity according to our framework. Algorithm performance roughly clusters according to criteria (a)(d) across tasks, lending credence to our approach. Curiously, gradient alignment with exact methods (RTRL and BPTT) does not predict performance, despite its ubiquity as a tool for analyzing approximate learning algorithms.
2 Past and futurefacing perspectives of online learning
Before we dive into the details of these algorithms, we first articulate what we mean by past and futurefacing, related to the “reverse/forward accumulation” distinction concurrently described by Cooijmans and Martens (2019). Consider a recurrent neural network that contains, at each time step $t$, a state ${\mathbf{a}}^{(t)}\in {\mathbb{R}}^{n}$. This state is updated via a function ${F}_{\mathbf{w}}:{\mathbb{R}}^{m}\to {\mathbb{R}}^{n}$, which is parameterized by a flattened vector of parameters $\mathbf{w}\in {\mathbb{R}}^{P}$. Here $m=n+{n}_{in}+1$ counts the total number of input dimensions, including the recurrent inputs ${\mathbf{a}}^{(t1)}\in {\mathbb{R}}^{n}$, task inputs ${\mathbf{x}}^{(t)}\in {\mathbb{R}}^{{n}_{\text{in}}}$, and an additional input clamped to $1$ (to represent bias). For some initial state ${\mathbf{a}}^{(0)}$, ${F}_{\mathbf{w}}$ defines the network dynamics by
$${\mathbf{a}}^{(t)}={F}_{\mathbf{w}}({\mathbf{a}}^{(t1)},{\mathbf{x}}^{(t)}).$$ 
At each time step an output ${\mathbf{y}}^{(t)}\in {\mathbb{R}}^{{n}_{\text{out}}}$ is computed by another function ${F}_{{\mathbf{w}}_{\text{o}}}^{\text{out}}:{\mathbb{R}}^{n}\to {\mathbb{R}}^{{n}_{out}}$, parameterized by ${\mathbf{w}}_{\text{o}}\in {\mathbb{R}}^{{P}_{\text{o}}}$. We will typically choose an affinesoftmax readout for ${F}_{{\mathbf{w}}_{\text{o}}}^{\text{out}}$, with output weights/bias ${\mathbf{W}}^{\text{out}}\in {\mathbb{R}}^{{n}_{\text{out}}\times (n+1)}$. A loss function $L({\mathbf{y}}^{(t)},{\mathbf{y}}^{*(t)})$ calculates an instantaneous loss ${L}^{(t)}$, quantifying to what degree the predicted output ${\mathbf{y}}^{(t)}$ matches the target output ${\mathbf{y}}^{*(t)}$.
The goal is to train the network by gradient descent (or other gradientbased optimizers such as ADAM from Kingma and Ba, 2014) on the total loss $\mathcal{L}=\sum {L}^{(t)}$ w.r.t. the parameters $\mathbf{w}$ and ${\mathbf{w}}_{\text{o}}$. It is natural to learn ${\mathbf{w}}_{\text{o}}$ online, because only information at present time $t$ is required to calculate the gradient $\partial {L}^{(t)}/\partial {\mathbf{w}}_{\text{o}}$. So the heart of the problem is to calculate $\partial \mathcal{L}/\partial \mathbf{w}$.
The parameter $\mathbf{w}$ is applied via ${F}_{\mathbf{w}}$ at every time step, and we denote a particular application of $\mathbf{w}$ at time $s$ as ${\mathbf{w}}^{(s)}$. Of course, a recurrent system is constrained to share parameters across time steps, so a perturbation $\delta \mathbf{w}$ is effectively a perturbation across all applications $\delta {\mathbf{w}}^{(s)}$, i.e., $\partial {\mathbf{w}}^{(s)}/\partial \mathbf{w}={\mathbf{I}}_{P}$. In principle, each application of the parameters affects all future losses ${L}^{(t)}$, $t\ge s$. The core of any recurrent learning algorithm is to estimate the influence $\partial {L}^{(t)}/\partial {\mathbf{w}}^{(s)}$ of one parameter application ${\mathbf{w}}^{(s)}$ on one loss ${L}^{(t)}$, since these individual terms are necessary and sufficient to define the global gradient
$$\frac{\partial \mathcal{L}}{\partial \mathbf{w}}=\sum _{t}\frac{\partial {L}^{(t)}}{\partial \mathbf{w}}=\sum _{t}\sum _{s\le t}\frac{\partial {L}^{(t)}}{\partial {\mathbf{w}}^{(s)}}\frac{\partial {\mathbf{w}}^{(s)}}{\partial \mathbf{w}}=\sum _{t}\sum _{s\le t}\frac{\partial {L}^{(t)}}{\partial {\mathbf{w}}^{(s)}}.$$  (1) 
This raises the question of how to sum these components to produce individual gradients to pass to the optimizer. In truncated BPTT, one unrolls the graph over some range of time steps and sums $\partial {L}^{(t)}/\partial {\mathbf{w}}^{(s)}$ for all $t,s$ in that range with $t\ge s$ (see §4.1.1). This does not qualify as an “online” learning rule, because it requires two independent time indices—at most one can represent “real time” leaving the other to represent the future or the past. If we can account for one of the summations via dynamic updates, then the algorithm is online or temporally local, i.e. not requiring explicit reference to the past or future. As depicted in Fig. 1, there are two possibilities. If $t$ from Eq. (1) corresponds to real time, then the gradient passed to the optimizer is
$${\nabla}_{\mathbf{w}}\mathcal{L}(t)=\sum _{s=0}^{t}\frac{\partial {L}^{(t)}}{\partial {\mathbf{w}}^{(s)}}=\frac{\partial {L}^{(t)}}{\partial \mathbf{w}}.$$  (2) 
In this case, we say learning is past facing, because the gradient is a sum of the influences of past applications of $\mathbf{w}$ on the current loss. On the other hand, if $s$ from Eq. (1) represents real time, then the gradient passed to the optimizer is
$${\nabla}_{\mathbf{w}}\mathcal{L}(s)=\sum _{t=s}^{\mathrm{\infty}}\frac{\partial {L}^{(t)}}{\partial {\mathbf{w}}^{(s)}}=\frac{\partial \mathcal{L}}{\partial {\mathbf{w}}^{(s)}}.$$  (3) 
Here we say learning is future facing, because the gradient is a sum of influences by the current application of $\mathbf{w}$ on future losses.
2.1 Pastfacing online learning algorithms
Here we derive a fundamental relation leveraged by pastfacing (PF) online algorithms. Let $t$ index real time, and define the influence matrix ${\mathbf{M}}^{(t)}\in {\mathbb{R}}^{n\times P}$, where $n$ and $P$ are respectively the number of hidden units and the number of parameters defining ${F}_{\mathbf{w}}$. ${\mathbf{M}}^{(t)}$ tracks the derivatives of the current state ${\mathbf{a}}^{(t)}$ with respect to each parameter ${w}_{p}$:
$${M}_{kp}^{(t)}=\frac{\partial {a}_{k}^{(t)}}{\partial {w}_{p}}.$$  (4) 
Let’s rewrite Eq. (4) with matrix notation and unpack it by one time step:
${\mathbf{M}}^{(t)}$  $={\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial \mathbf{w}}}={\displaystyle \sum _{s\le t}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(s)}}}={\displaystyle \sum _{s\le t1}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(s)}}}+{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(t)}}}$  
$={\displaystyle \sum _{s\le t1}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{a}}^{(t1)}}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t1)}}{\partial {\mathbf{w}}^{(s)}}}+{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(t)}}}$  
$={\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{a}}^{(t1)}}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t1)}}{\partial \mathbf{w}}}+{\displaystyle \frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(t)}}}$  
$\equiv {\mathbf{J}}^{(t)}{\mathbf{M}}^{(t1)}+{\overline{\mathbf{M}}}^{(t)}.$  (5) 
A simple recursive formula emerges, wherein the influence matrix is updated by multiplying its current value by the Jacobian ${\mathbf{J}}^{(t)}=\partial {\mathbf{a}}^{(t)}/\partial {\mathbf{a}}^{(t1)}\in {\mathbb{R}}^{n\times n}$ of the network and then adding the immediate influence ${\overline{\mathbf{M}}}^{(t)}=\partial {\mathbf{a}}^{(t)}/\partial {\mathbf{w}}^{(t)}\in {\mathbb{R}}^{n\times P}$. To compute the gradient that ultimately gets passed to the optimizer, we simply use the chain rule over the current hidden state ${\mathbf{a}}^{(t)}$:
$$\frac{\partial {L}^{(t)}}{\partial \mathbf{w}}=\frac{\partial {L}^{(t)}}{\partial {\mathbf{a}}^{(t)}}\frac{\partial {\mathbf{a}}^{(t)}}{\partial \mathbf{w}}\equiv {\overline{\mathbf{c}}}^{(t)}{\mathbf{M}}^{(t)},$$  (6) 
where the immediate credit assignment vector ${\overline{\mathbf{c}}}^{(t)}\in {\mathbb{R}}^{n}$ is defined to be $\partial {L}^{(t)}/\partial {\mathbf{a}}^{(t)}$ and is calculated by backpropagating the error ${\bm{\delta}}^{(t)}$ through the derivative of the output function ${F}_{{\mathbf{w}}_{\text{o}}}^{\text{out}}$ (or approximated by Feedback Alignment, see Lillicrap et al., 2016). In the end, we compute a derivative in Eq. (6) that is implicitly a sum over the many terms of Eq. (2), using formulae that depend explicitly only on times $t$ and $t1$. For this reason, such a learning algorithm is online, and it is past facing because the gradient computation is of the form in Eq. (2).
2.2 Futurefacing online learning algorithms
Here we show a symmetric relation for futurefacing (FF) online algorithms. The credit assignment vector ${\mathbf{c}}^{(t)}\in {\mathbb{R}}^{n}$ is a row vector defined as the gradient of the loss $\mathcal{L}$ with respect to the hidden state ${\mathbf{a}}^{(t)}$. It plays a role analogous to ${\mathbf{M}}^{(t)}$ and has a recursive update similar to Eq. (5):
${\mathbf{c}}^{(t)}$  $={\displaystyle \frac{\partial \mathcal{L}}{\partial {\mathbf{a}}^{(t)}}}={\displaystyle \sum _{s\ge t}}{\displaystyle \frac{\partial {L}^{(s)}}{\partial {\mathbf{a}}^{(t)}}}={\displaystyle \frac{\partial {L}^{(t)}}{\partial {\mathbf{a}}^{(t)}}}+{\displaystyle \sum _{s\ge t+1}}{\displaystyle \frac{\partial {L}^{(s)}}{\partial {\mathbf{a}}^{(t)}}}$  
$={\displaystyle \frac{\partial {L}^{(t)}}{\partial {\mathbf{a}}^{(t)}}}+{\displaystyle \sum _{s\ge t+1}}{\displaystyle \frac{\partial {L}^{(s)}}{\partial {\mathbf{a}}^{(t+1)}}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t+1)}}{\partial {\mathbf{a}}^{(t)}}}$  
$={\displaystyle \frac{\partial {L}^{(t)}}{\partial {\mathbf{a}}^{(t)}}}+{\displaystyle \frac{\partial \mathcal{L}}{\partial {\mathbf{a}}^{(t+1)}}}{\displaystyle \frac{\partial {\mathbf{a}}^{(t+1)}}{\partial {\mathbf{a}}^{(t)}}}$  
$={\overline{\mathbf{c}}}^{(t)}+{\mathbf{c}}^{(t+1)}{\mathbf{J}}^{(t+1)}.$  (7) 
As in the PF case, the gradient is ultimately calculated using the chain rule over ${\mathbf{a}}^{(t)}$:
$$\frac{\partial \mathcal{L}}{\partial {\mathbf{w}}^{(t)}}=\frac{\partial \mathcal{L}}{\partial {\mathbf{a}}^{(t)}}\frac{\partial {\mathbf{a}}^{(t)}}{\partial {\mathbf{w}}^{(t)}}\equiv {\mathbf{c}}^{(t)}{\overline{\mathbf{M}}}^{(t)}.$$  (8) 
The recursive relations for PF and FF algorithms are of identical form given the following changes: (1) swap the roles of $\mathcal{L}$ and $\mathbf{w}$, (2) swap the roles of $t1$ and $t+1$, and (3) flip the direction of all derivatives. This clarifies the fundamental tradeoff between the PF and FF approaches to online learning. On the one hand, memory requirements favor FF because $\mathcal{L}$ is a scalar while $\mathbf{w}$ is a matrix. On the other, only PF can truly be run online, because the time direction of the update in FF is opposite the forward pass. Thus, efficient PF algorithms must compress ${\mathbf{M}}^{(t)}$, while efficient FF algorithms must predict ${\mathbf{c}}^{(t+1)}$.
3 Pastfacing algorithms
3.1 RealTime Recurrent Learning
The RealTime Recurrent Learning (RTRL, Williams and Zipser, 1989) algorithm directly applies Eqs. (5) and (6) as written. We call the application of Eq. (5) the “update” to the learning algorithm, which is deterministic and in closed form. Implementing Eq. (5) requires storing $nP\approx \mathcal{O}({n}^{3})$ floats in ${\mathbf{M}}^{(t)}$ and performing $\mathcal{O}({n}^{4})$ multiplications in ${\mathbf{J}}^{(t)}{\mathbf{M}}^{(t)}$, which is neither especially efficient nor biologically plausible. However, several efficient (and in some cases, biologically plausible) online learning algorithms have recently been developed, including Unbiased Online Recurrent Optimization (UORO; Tallec and Ollivier, 2017; §3.2), KroneckerFactored RTRL (KFRTRL; Mujika et al., 2018; §3.3), Kernel RNN Learning (KeRNL; Roth et al., 2019; §3.5), and RandomFeedback Online Learning (RFLO; Murray, 2019; §3.6). We claim that these learning algorithms, whether explicitly derived as such or not, are all implicitly approximations to RTRL, each a special case of a general class of techniques for compressing ${\mathbf{M}}^{(t)}$. In the following section, we clarify how each of these learning algorithms fits into this broad structure.
3.1.1 Approximations to RTRL
To concretely illuminate these ideas, we will work with a special case of ${F}_{\mathbf{w}}$, a timecontinuous vanilla RNN:
$${\mathbf{a}}^{(t)}={F}_{\mathbf{w}}({\mathbf{a}}^{(t1)},{\mathbf{x}}^{(t)})=(1\alpha ){\mathbf{a}}^{(t1)}+\alpha \varphi (\mathbf{W}{\widehat{\mathbf{a}}}^{(t1)}),$$  (9) 
where ${\widehat{\mathbf{a}}}^{(t1)}=\text{concat}({\mathbf{a}}^{(t1)},{\mathbf{x}}^{(t)},1)\in {\mathbb{R}}^{m}$, $\mathbf{W}\in {\mathbb{R}}^{n\times m}$, $\varphi :{\mathbb{R}}^{n}\to {\mathbb{R}}^{n}$ is some pointwise nonlinearity (e.g. $\mathrm{tanh}$), and $\alpha \in (0,1]$ is the network’s inverse time constant. The trainable parameters ${w}_{p}$ are folded via the indexing $p=i\times n+j$ into the weight matrix ${W}_{ij}$, whose columns hold the recurrent weights, the input weights, and a bias. By reshaping ${w}_{p}$ into its natural matrix form ${W}_{ij}$, we can write the influence matrix as an order3 influence tensor
$${M}_{kij}^{(t)}=\partial {a}_{k}^{(t)}/\partial {W}_{ij}.$$ 
Thus ${M}_{kij}^{(t)}$ specifies the effect on the $k$th unit of perturbing the direct connection from the $j$th unit to the $i$th unit. The immediate influence can also be written as a tensor. By differentiating Eq. (9), we see it takes the sparse form
$${\overline{M}}_{kij}^{(t)}=\partial {a}_{k}^{(t)}/\partial {W}_{ij}^{(t)}=\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)},$$ 
because ${W}_{ij}$ can affect the $k$th unit directly only if $k=i$. Many approximations of RTRL involve a decomposition of ${M}_{kij}^{(t)}$ into a product of lowerorder tensors. For example, UORO represents ${M}_{kij}^{(t)}$ by an outer product ${A}_{k}^{(t)}{B}_{ij}^{(t)}$, which has a memory requirement of only $\mathcal{O}({n}^{2})$. Similarly, KFRTRL uses a Kroneckerproduct decomposition ${A}_{j}^{(t)}{B}_{ki}^{(t)}$. We can generalize these cases into a set of six possible decompositions of ${M}_{kij}^{(t)}$ into products of lowerorder tensors ${A}^{(t)}$ and ${B}^{(t)}$:
$${M}_{kij}^{(t)}\approx \{\begin{array}{cc}{A}_{k}^{(t)}{B}_{ij}^{(t)}\hfill & \text{UORO},\mathrm{\S}\mathit{\text{3.2}}\hfill \\ {A}_{j}^{(t)}{B}_{ki}^{(t)}\hfill & \text{KFRTRL},\mathrm{\S}\mathit{\text{3.3}}\hfill \\ {A}_{i}^{(t)}{B}_{kj}^{(t)}\hfill & \text{\u201cReverse\u201d KFRTRL},\mathrm{\S}\mathit{\text{3.4}}\hfill \\ {A}_{ki}^{(t)}{B}_{ij}^{(t)}\hfill & \text{KeRNL/RFLO},\mathrm{\S}\mathit{\text{3.5}}/\mathrm{\S}\mathit{\text{3.6}}\hfill \\ {A}_{kj}^{(t)}{B}_{ij}^{(t)}\hfill & \text{Unexplored}\hfill \\ {A}_{ki}^{(t)}{B}_{kj}^{(t)}\hfill & \text{Unexplored}\hfill \end{array}.$$ 
Each such decomposition has a memory requirement of $\mathcal{O}({n}^{2})$. Of course, it is not sufficient to write down an idealized decomposition for a particular time point; there must exist some efficient way to update the decomposition as the network runs forwards. We now go through each algorithm and show the mathematical techniques used to derive update equations and categorize them by the criteria outlined in Table 1.
3.2 Unbiased Online Recurrent Optimization (UORO)
Tallec and Ollivier (2017) discovered a technique for approximating ${\mathbf{M}}^{(t)}\in {\mathbb{R}}^{n\times P}$ as an outer product ${\mathbf{A}}^{(t)}{\mathbf{B}}^{(t)}$, where ${\mathbf{A}}^{(t)}\in {\mathbb{R}}^{n\times 1}$ and ${\mathbf{B}}^{(t)}\in {\mathbb{R}}^{1\times P}$. The authors proved a crucial lemma (see Appendix A or Tallec and Ollivier, 2017) that gives, in closed form, an unbiased rank1 estimate of a given matrix over the choice of a random vector $\bm{\nu}\in {\mathbb{R}}^{n}$ with $\mathbb{E}[{\nu}_{i}{\nu}_{j}]\propto {\delta}_{ij}$ and $\mathbb{E}[{\nu}_{i}]=0$. They leverage this result to derive a closedform update rule for ${\mathbf{A}}^{(t)}$ and ${\mathbf{B}}^{(t)}$ at each time step, without ever having to explicitly (and expensively) calculate ${\mathbf{M}}^{(t)}$. We present an equivalent formulation in terms of tensor components, i.e.,
$${M}_{kij}^{(t)}\approx {A}_{k}^{(t)}{B}_{ij}^{(t)},$$ 
where ${B}_{ij}^{(t)}$ represents the “rolledup” components of ${\mathbf{B}}^{(t)}$, as in ${W}_{ij}$ w.r.t. $\mathbf{w}$. Intuitively, the $kij$th component of the influence matrix is constrained to be the product of the $k$th unit’s “sensitivity” ${A}_{k}^{(t)}$ and the $ij$th parameter’s “efficacy” ${B}_{ij}^{(t)}$. Eqs. (10) and (11) show the form of the update and why it is unbiased over $\bm{\nu}$, respectively:
${A}_{k}^{(t)}{B}_{ij}^{(t)}$  $=\left({\rho}_{0}{\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{A}_{{k}^{\prime}}^{(t1)}+{\rho}_{1}{\nu}_{k}\right)\left({\rho}_{0}^{1}{B}_{ij}^{(t1)}+{\rho}_{1}^{1}{\displaystyle \sum _{{k}^{\prime}}}{\nu}_{{k}^{\prime}}{\overline{M}}_{{k}^{\prime}ij}^{(t)}\right)$  
$={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{A}_{{k}^{\prime}}^{(t1)}{B}_{ij}^{(t1)}+{\displaystyle \sum _{{k}^{\prime}}}{\nu}_{k}{\nu}_{{k}^{\prime}}{\overline{M}}_{{k}^{\prime}ij}^{(t)}$  
$+{\displaystyle \sum _{{k}^{\prime}}}{\nu}_{{k}^{\prime}}\left[{\rho}_{1}{\rho}_{0}^{1}{\delta}_{k{k}^{\prime}}{B}_{ij}^{(t1)}+{\rho}_{0}{\rho}_{1}^{1}{\overline{M}}_{{k}^{\prime}ij}^{(t)}{\displaystyle \sum _{{k}^{\prime \prime}}}{J}_{{k}^{\prime}{k}^{\prime \prime}}^{(t)}{A}_{{k}^{\prime \prime}}^{(t1)}\right]$  (10)  
$\u27f9\mathbb{E}\left[{A}_{k}^{(t)}{B}_{ij}^{(t)}\right]$  $={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}\mathbb{E}\left[{A}_{{k}^{\prime}}^{(t1)}{B}_{ij}^{(t1)}\right]+{\displaystyle \sum _{{k}^{\prime}}}\mathbb{E}[{\nu}_{k}{\nu}_{{k}^{\prime}}]{\overline{M}}_{{k}^{\prime}ij}^{(t)}$  
$+{\displaystyle \sum _{{k}^{\prime}}}\mathbb{E}[{\nu}_{{k}^{\prime}}]\left(\text{cross terms}\right)$  
$={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\displaystyle \sum _{{k}^{\prime}}}{\delta}_{k{k}^{\prime}}{\overline{M}}_{{k}^{\prime}ij}^{(t)}+{\displaystyle \sum _{{k}^{\prime}}}0\times (\text{cross terms})$  
$={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$={M}_{kij}^{(t)}.$  (11) 
The cross terms vanish in expectation because $\mathbb{E}[{\nu}_{k}]=0$. Thus, by induction over $t$, the estimate of ${M}_{kij}^{(t)}$ remains unbiased at every time step. The constants ${\rho}_{0},{\rho}_{1}\in {\mathbb{R}}^{>0}$ are chosen at each time step to minimize total variance of the estimate by balancing the norms of the cross terms. This algorithm’s update is stochastic due to its reliance on the random vector $\bm{\nu}$, but it is in closed form because it has an explicit update formula (Eq. 10). Both its memory and computational complexity are $\mathcal{O}({n}^{2})$.
3.3 KroneckerFactored RTRL (KFRTRL)
Mujika et al. (2018) leverage the same lemma as in UORO, but using a decomposition of ${\mathbf{M}}^{(t)}$ in terms of a Kronecker product ${\mathbf{A}}^{(t)}\otimes {\mathbf{B}}^{(t)}$, where now ${\mathbf{A}}^{(t)}\in {\mathbb{R}}^{1\times m}$ and ${\mathbf{B}}^{(t)}\in {\mathbb{R}}^{n\times n}$. This decomposition is more natural, because the immediate influence ${\overline{\mathbf{M}}}^{(t)}$ factors exactly as a Kronecker product ${\widehat{\mathbf{a}}}^{(t)}\otimes {\mathbf{D}}^{(t)}$ for vanilla RNNs, where ${D}_{ki}^{(t)}=\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)})$. To derive the update rule for UORO, one must first generate a rank1 estimate of ${\overline{\mathbf{M}}}^{(t)}$ as an intermediate step, introducing more variance, but in KFRTRL, this step is unnecessary. In terms of components, the compression takes the form
$${M}_{kij}^{(t)}\approx {A}_{j}^{(t)}{B}_{ki}^{(t)},$$ 
which is similar to UORO, modulo a cyclic permutation of the indices. Given a sample $\bm{\nu}\in {\mathbb{R}}^{2}$ of only 2 i.i.d. random variables, again with $\mathbb{E}[{\nu}_{i}{\nu}_{j}]\propto {\delta}_{ij}$ and $\mathbb{E}[{\nu}_{i}]=0$, the update takes the form shown in Eqs. (12) and (13):
${A}_{j}^{(t)}$  $=\left({\nu}_{0}{\rho}_{0}{A}_{j}^{(t1)}+{\nu}_{1}{\rho}_{1}{\widehat{a}}_{j}^{(t1)}\right)$  (12)  
${B}_{ki}^{(t)}$  $=\left({\nu}_{0}{\rho}_{0}^{1}{\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{B}_{{k}^{\prime}i}^{(t1)}+{\nu}_{1}{\rho}_{1}^{1}\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)})\right)$  (13)  
$\u27f9{A}_{j}^{(t)}{B}_{ki}^{(t)}$  $={\nu}_{0}^{2}{\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{A}_{j}^{(t1)}{B}_{{k}^{\prime}i}^{(t1)}+{\nu}_{1}^{2}\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}+\text{crossterms}$  
$\u27f9\mathbb{E}\left[{A}_{j}^{(t)}{B}_{ki}^{(t)}\right]$  $={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}\mathbb{E}\left[{A}_{j}^{(t1)}{B}_{{k}^{\prime}i}^{(t1)}\right]+\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$  
$={\displaystyle \sum _{k{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$={M}_{kij}^{(t)}.$ 
As in UORO, the cross terms vanish in expectation, and the estimate is unbiased by induction over $t$. This algorithm’s updates are also stochastic and in closed form. Its memory complexity is $\mathcal{O}({n}^{2})$, but its computation time is $\mathcal{O}({n}^{3})$ because of the matrixmatrix product in Eq. (13).
3.4 Reverse KFRTRL (RKFRTRL)
Our exploration of the space of different approximations naturally raises a question: is an approximation of the form
$${M}_{kij}^{(t)}\approx {A}_{i}^{(t)}{B}_{kj}^{(t)}$$  (14) 
also possible? We refer to this method as “Reverse” KFRTRL (RKFRTRL) because, in matrix notation, this would be formulated as ${\mathbf{M}}^{(t)}\approx {\mathbf{B}}^{(t)}\otimes {\mathbf{A}}^{(t)}$, where ${\mathbf{A}}^{(t)}\in {\mathbb{R}}^{1\times n}$ and ${\mathbf{B}}^{(t)}\in {\mathbb{R}}^{n\times m}$. We propose the following update for ${A}_{i}^{(t)}$ and ${B}_{kj}^{(t)}$ in terms of a random vector $\bm{\nu}\in {\mathbb{R}}^{n}$:
${A}_{i}^{(t)}{B}_{kj}^{(t)}$  $=\left({\rho}_{0}{A}_{i}^{(t1)}+{\rho}_{1}{\nu}_{i}\right)\left({\rho}_{0}^{1}{\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{B}_{{k}^{\prime}j}^{(t1)}+{\rho}_{1}^{1}{\displaystyle \sum _{{i}^{\prime}}}{\nu}_{{i}^{\prime}}{\overline{M}}_{k{i}^{\prime}j}^{(t)}\right)$  (15)  
$={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{A}_{i}^{(t1)}{B}_{{k}^{\prime}j}^{(t1)}+{\displaystyle \sum _{{i}^{\prime}}}{\nu}_{i}{\nu}_{{i}^{\prime}}{\overline{M}}_{k{i}^{\prime}j}^{(t)}+\text{crossterms}$  
$\u27f9\mathbb{E}\left[{A}_{i}^{(t)}{B}_{kj}^{(t)}\right]$  $={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}\mathbb{E}\left[{A}_{i}^{(t1)}{B}_{{k}^{\prime}j}^{(t1)}\right]+{\overline{M}}_{kij}^{(t)}$  
$={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$={M}_{kij}^{(t)}.$  (16) 
Eq. (16) shows that this estimate is unbiased, using updates that are stochastic and in closed form, like its sibling algorithms. Its memory and computational complexity are $\mathcal{O}({n}^{2})$ and $\mathcal{O}({n}^{3})$, respectively. RKFRTRL is actually more similar to UORO than KFRTRL, because ${\overline{M}}_{kij}^{(t)}$ does not naturally factor like Eq. (14), introducing more variance. Worse, it has the computational complexity of KFRTRL due to the matrixmatrix multiplication in Eq. (15). KFRTRL stands out as the most effective of these 3 algorithms, because it estimates ${\mathbf{M}}^{(t)}$ with the lowest variance due to its natural decomposition structure. (See Mujika et al., 2018 for variance calculations.)
3.4.1 Optimal KroneckerSum Approximation (OK)
We briefly mention an extension of KFRTRL by Benzing et al. (2019), where the influence matrix is approximated not by 1 but rather a sum of $r$ Kronecker products, or, in components
$${M}_{kij}^{(t)}\approx \sum _{l=1}^{r}{A}_{lj}^{(t)}{B}_{lki}^{(t)}.$$ 
On the RTRL update, the $k$ index of ${B}_{lki}^{(t)}$ is propagated forward by the Jacobian, and then the immediate influence—itself a Kronecker product—is added. Now ${M}_{kij}^{(t)}$ is approximated by $r+1$ Kronecker products
$${M}_{kij}^{(t)}\approx \sum _{l=1}^{r}{A}_{lj}^{(t1)}{J}_{k{k}^{\prime}}^{(t)}{B}_{l{k}^{\prime}i}^{(t1)}+\alpha {\widehat{a}}_{j}^{(t1)}{\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}),$$ 
but the authors developed a technique to optimally reduce this sum back to $r$ Kronecker products, keeping the memory complexity $\mathcal{O}(r{n}^{2})$ and computational complexity $\mathcal{O}(r{n}^{3})$ constant. This update is stochastic because it requires explicit randomness in the flavor of the above algorithms, and it is numerical because there is no closed form solution to the update. We leave the details to the original paper.
3.5 Kernel RNN Learning (KeRNL)
Roth et al. (2019) developed a learning algorithm for RNNs that is essentially a compression of the influence matrix of the form ${M}_{kij}^{(t)}\approx {A}_{ki}{B}_{ij}^{(t)}$. We will show that this algorithm is also an implicit approximation of RTRL, although the update rules are fundamentally different than those for UORO, KFRTRL and RKFRTRL. The eligibility trace ${\mathbf{B}}^{(t)}\in {\mathbb{R}}^{n\times m}$ updates by temporally filtering the immediate influences $\alpha {\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$ with unitspecific, learnable timescales ${\alpha}_{i}$:
$${B}_{ij}^{(t)}=(1{\alpha}_{i}){B}_{ij}^{(t1)}+\alpha {\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}.$$  (17) 
The sensitivity matrix $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is chosen to approximate the multistep Jacobian $\partial {a}_{k}^{(t)}/\partial {a}_{i}^{({t}^{\prime})}$ with help from the learned timescales:
$$\frac{\partial {a}_{k}^{(t)}}{\partial {a}_{i}^{({t}^{\prime})}}\approx {A}_{ki}{(1{\alpha}_{i})}^{(t{t}^{\prime})}.$$  (18) 
We will describe how $\mathbf{A}$ is learned later, but for now we assume this approximation holds and use it to show how the KeRNL update is equivalent to that of RTRL. We have dropped the explicit timedependence from $\mathbf{A}$, because it updates too slowly for Eq. (18) to be specific to any one time point. If we unpack this approximation by one time step, we uncover the consistency relation
$${A}_{ki}(1{\alpha}_{i})\approx \sum _{{k}^{\prime}}{J}_{k{k}^{\prime}}^{(t)}{A}_{{k}^{\prime}i}.$$  (19) 
By taking $t={t}^{\prime}$ in Eq. (18) and rearranging Eq. (19), we see this approximation implicitly assumes both
$${A}_{ki}\approx \{\begin{array}{cc}{\delta}_{ki}\hfill & \\ {(1{\alpha}_{i})}^{1}{\sum}_{{k}^{\prime}}{J}_{k{k}^{\prime}}^{(t)}{A}_{{k}^{\prime}i}\hfill & \end{array}.$$  (20) 
Then the eligibility trace update effectively implements the RTRL update, assuming inductively that ${M}_{kij}^{(t1)}$ is well approximated by ${A}_{ki}{B}_{ij}^{(t1)}$:
${A}_{ki}{B}_{ij}^{(t)}$  $={A}_{ki}\left[(1{\alpha}_{i}){B}_{ij}^{(t1)}+\alpha {\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}\right]$  
$={A}_{ki}(1{\alpha}_{i}){B}_{ij}^{(t1)}+\alpha {A}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$  
$\approx {\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{A}_{{k}^{\prime}i}{B}_{ij}^{(t1)}+\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$  (21)  
$={\displaystyle \sum _{k{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$={M}_{kij}^{(t)}.$ 
In Eq. (21), we use each of the special cases from Eq. (20). Of course, the ${A}_{ki}$ and ${\alpha}_{i}$ have to be learned, and Roth et al. (2019) use gradient descent to do so. We leave details to the original paper; briefly, they run in parallel a perturbed forward trajectory to estimate the LHS of Eq. (18) and then perform SGD on the squared difference between the LHS and RHS, giving gradients for ${A}_{ki}$ and ${\alpha}_{i}$.
KeRNL uses deterministic updates because it does not need explicit random variables. While the ${B}_{ij}^{(t)}$ update is in closed form via Eq. (17), the updates for ${A}_{ki}$ and ${\alpha}_{i}$ are numerical because of the need for SGD to train them to obey Eq. (18). Both its memory and computational complexities are $\mathcal{O}({n}^{2})$.
3.6 RandomFeedback Online Learning (RFLO)
Coming from a computational neuroscience perspective, Murray (2019) developed a beautifully simple and biologically plausible learning rule for RNNs, which he calls RandomFeedback Online Learning (RFLO). He formulates the rule in terms of an eligibility trace ${B}_{ij}^{(t)}$ that filters the nonzero immediate influence elements ${\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$ by the network inverse time constant $\alpha $:
$${B}_{ij}^{(t)}=(1\alpha ){B}_{ij}^{(t1)}+\alpha {\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}.$$ 
Then the approximate gradient is ultimately calculated^{1}^{1} 1 As the “random feedback” part of the name suggests, Murray goes a step further in approximating ${\overline{c}}_{k}^{(t)}$ by random feedback weights á la Lillicrap et al., 2016, but we assume exact feedback in this paper for easier comparisons with other algorithms. as
$$\frac{\partial {L}^{(t)}}{\partial {W}_{ij}}\approx {\overline{c}}_{i}^{(t)}{B}_{ij}^{(t)}.$$ 
By observing that
$${\overline{c}}_{i}^{(t)}{B}_{ij}^{(t)}=\sum _{k}{\overline{c}}_{k}^{(t)}{\delta}_{ki}{B}_{ij}^{(t)},$$ 
we see that RFLO is a special case of KeRNL, in which we fix ${A}_{ki}={\delta}_{ki}$, ${\alpha}_{i}=\alpha $. Alternatively, and as hinted in the original paper, we can view RFLO as a special case of RTRL under the approximation ${J}_{k{k}^{\prime}}^{(t)}\approx (1\alpha ){\delta}_{k{k}^{\prime}}$, because the RTRL update reduces to RFLO with ${M}_{kij}^{(t)}={\delta}_{ki}{B}_{ij}^{(t)}$ containing ${B}_{ij}^{(t)}$ along the diagonals:
${M}_{kij}^{(t)}$  $={\displaystyle \sum _{{k}^{\prime}}}{J}_{k{k}^{\prime}}^{(t)}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$=(1\alpha ){\displaystyle \sum _{{k}^{\prime}}}{\delta}_{k{k}^{\prime}}{M}_{{k}^{\prime}ij}^{(t1)}+{\overline{M}}_{kij}^{(t)}$  
$=(1\alpha ){M}_{kij}^{(t1)}+\alpha {\delta}_{ki}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}.$  (22) 
Fig. 2 illustrates how ${\mathbf{B}}^{(t)}$ is contained in the influence matrix ${\mathbf{M}}^{(t)}$. This algorithm’s update is deterministic and in closed form, with memory and computational complexity $\mathcal{O}({n}^{2})$.
4 Futurefacing algorithms
4.1 Backpropagation Through Time (BPTT)
For many applications, a recurrent network is unrolled only for some finite number of time steps, and backpropagation through time (BPTT) manifests as the computation of the sum $\partial {L}^{(t)}/\partial {\mathbf{w}}^{(s)}$ over every $s\le t$ in the graph. This can be efficiently accomplished using
$${\mathbf{c}}^{(t)}={\overline{\mathbf{c}}}^{(t)}+{\mathbf{J}}^{(t+1)}{\mathbf{c}}^{(t+1)}$$  (23) 
(see Eq. 7) to propagate credit assignment backwards. However, in our framework, where a network is run on an infinitetime horizon, there are two qualitatively different ways of unrolling the network. We call them “efficient” and “futurefacing” BPTT.
4.1.1 Efficient backpropagation through time (EBPTT)
For this method, we simply divide the graph into nonoverlapping segments of truncation length $T$ and perform BPTT between $tT$ and $t$ as described above, using Eq. (23). It takes $\mathcal{O}({n}^{2}T)$ computation time to compute one gradient, but since this computation is only performed once every $T$ time steps, the computation time is effectively $\mathcal{O}({n}^{2})$, with memory requirement $\mathcal{O}(nT)$. A problem with this approach is that it does not treat all time points the same: an application of $\mathbf{w}$ occurring near the end of the graph segment has less of its future influence accounted for than applications of $\mathbf{w}$ occurring before it, as can be visualized in Fig. 3. And since any one gradient passed to the optimizer is a sum across both $t$ and $s$, it is not an online algorithm by the framework we presented in §2. Therefore, for the purpose of comparing with online algorithms, we also show an alternative version of BPTT that calculates a futurefacing gradient (up to truncation) $\partial \mathcal{L}/\partial {\mathbf{w}}^{(t)}$ for every $t$.
4.1.2 Futurefacing backpropagation through time (FBPTT)
In this version of BPTT, we keep a dynamic list of truncated credit assignment estimates ${\widehat{\mathbf{c}}}^{(s)}$ for times $s=tT,\mathrm{\cdots},t1$:
$$[{\widehat{\mathbf{c}}}^{(tT)},\mathrm{\cdots},{\widehat{\mathbf{c}}}^{(t1)}],$$ 
where each truncated credit assignment estimate includes the influences of ${\mathbf{a}}^{(s)}$ only up to time $t1$:
$${\widehat{\mathbf{c}}}^{(s)}=\sum _{{t}^{\prime}=s}^{t1}\frac{\partial {L}^{({t}^{\prime})}}{\partial {\mathbf{a}}^{(s)}}.$$ 
At current time $t$, every element ${\widehat{\mathbf{c}}}^{(s)}$ is extended by adding $\partial {L}^{(t)}/\partial {\mathbf{a}}^{(s)}$, calculated by backpropagating from the current loss ${L}^{(t)}$, while the explicit credit assignment ${\overline{\mathbf{c}}}^{(t)}$ is appended to the front of the list. To compensate, the oldest credit assignment estimate ${\widehat{\mathbf{c}}}^{(tT)}$ is removed and combined with the immediate influence to form a (truncated) gradient
$${\widehat{\mathbf{c}}}^{(tT)}{\overline{\mathbf{M}}}^{(tT)}=\sum _{{t}^{\prime}=tT}^{t}\frac{\partial {L}^{({t}^{\prime})}}{\partial {\mathbf{a}}^{(tT)}}\frac{\partial {\mathbf{a}}^{(tT)}}{\partial {\mathbf{w}}^{(tT)}}=\sum _{{t}^{\prime}=tT}^{t}\frac{\partial {L}^{({t}^{\prime})}}{\partial {\mathbf{w}}^{(tT)}}\approx \frac{\partial \mathcal{L}}{\partial {\mathbf{w}}^{(tT)}},$$ 
which is passed to the optimizer to update the network. This algorithm is “online” in that it produces strictly futurefacing gradients at each time step, albeit delayed by the truncation time $T$ and requiring memory of the network states from $tT$. Each update step requires $\mathcal{O}({n}^{2}T)$ computation, but since the update is performed at every time step, computation remains a factor of $T$ more expensive than EBPTT. Memory requirement is still $\mathcal{O}(nT)$. Fig. 3 illustrates the differences among these methods and RTRL, using a triangular lattice as a visualization tool. Each point in the lattice is one derivative $\partial {L}^{(t)}/\partial {\mathbf{w}}^{(s)}$ with $t\ge s$, and the points are grouped together into discrete gradients passed to the optimizer.
4.2 Decoupled Neural Interfaces (DNI)
Jaderberg et al. (2017) developed a framework for online learning by predicting credit assignment. Whereas PF algorithms face the problem of a large influence tensor ${M}_{kij}^{(t)}$ that needs a compressed representation, FF algorithms face the problem of incomplete information: at time $t$, it is impossible to calculate ${\mathbf{c}}^{(t)}$ without access to future network variables. The approach of Decoupled Neural Interfaces (DNI) is to simply make a linear prediction of ${\mathbf{c}}^{(t)}$ (Czarnecki et al., 2017) based on the current hidden state ${\mathbf{a}}^{(t)}$ and the current labels ${\mathbf{y}}^{*(t)}$:
$${c}_{i}^{(t)}\approx \sum _{l}{\stackrel{~}{a}}_{l}^{(t)}{A}_{li},$$ 
where ${\stackrel{~}{\mathbf{a}}}^{(t)}=\text{concat}({\mathbf{a}}^{(t)},{\mathbf{y}}^{*(t)},1)\in {\mathbb{R}}^{{m}^{\prime}}$, ${m}^{\prime}=n+{n}_{\text{out}}+1$, and ${A}_{li}$ are the components of a matrix $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, which parameterizes what the authors call the synthetic gradient function. The parameters ${A}_{li}$ are trained to minimize the loss
$${L}_{\text{SG}}^{(t)}=\frac{1}{2}{\parallel \sum _{l}{\stackrel{~}{a}}_{l}^{(t)}{A}_{li}{c}_{i}^{(t)}\parallel}^{2}$$  (24) 
via gradient descent, similar to KeRNL’s treatment of ${A}_{ki}$ and ${\alpha}_{i}$ (and we drop the time dependence of ${A}_{li}$ for the same reason). Of course, this begs the question—the whole point is to avoid calculating ${\mathbf{c}}^{(t)}$ explicitly, but calculating the error in Eq. (24) requires access to ${\mathbf{c}}^{(t)}$. So the authors propose a “bootstrapping” technique analogous to the Bellman equation in Reinforcement Learnin (Sutton and Barto, 2018). If we take the FF relation we derived in Eq. (7)
$${\mathbf{c}}^{(t)}={\overline{\mathbf{c}}}^{(t)}+{\mathbf{c}}^{(t+1)}{\mathbf{J}}^{(t+1)}$$  (25) 
and approximate the appearance of ${\mathbf{c}}^{(t+1)}$ with the synthetic gradient estimate ${\stackrel{~}{\mathbf{a}}}^{(t+1)}\mathbf{A}$, then Eq. (25) provides an estimate of ${c}_{i}^{(t)}$ to use in Eq. (24). Then the update for $\mathbf{A}$ can be written as
$$\mathrm{\Delta}{A}_{li}\propto {\stackrel{~}{a}}_{l}^{(t)}\left[\sum _{{l}^{\prime}}{\stackrel{~}{a}}_{{l}^{\prime}}^{(t)}{A}_{{l}^{\prime}i}\left({\overline{c}}_{i}^{(t)}+\sum _{m}\sum _{{l}^{\prime}}{\stackrel{~}{a}}_{{l}^{\prime}}^{(t+1)}{A}_{{l}^{\prime}m}{J}_{mi}^{(t+1)}\right)\right]$$  (26) 
with learning rate chosen as a hyperparameter. As in Eq. (8), the gradient is calculated by combining the estimated credit assignment for the $i$th unit with the explicit influence by the $ij$th parameter:
$$\frac{\partial \mathcal{L}}{\partial {W}_{ij}^{(t)}}=\frac{\partial \mathcal{L}}{\partial {a}_{i}^{(t)}}\frac{\partial {a}_{i}^{(t)}}{\partial {W}_{ij}^{(t)}}={\overline{c}}_{i}^{(t)}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}\approx \left(\sum _{l}{\stackrel{~}{a}}_{l}^{(t)}{A}_{li}\right){\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$$ 
This algorithm is future facing because it ultimately estimates the effect of applying $\mathbf{w}$ at current time $t$ on future losses. Its updates are deterministic and numerical, because no explicit randomness is required, but the minimization problem over ${A}_{li}$ implied by Eq. (24) is approximated by gradient descent rather than solved in closed form. It requires $\mathcal{O}({n}^{2})$ memory for $\mathbf{A}$ and $\mathcal{O}({n}^{2})$ computation for the matrixvector multiplications in Eq. (26).
4.2.1 Biological approximation to DNI
While many of the algorithms we have presented are biologically plausible in the abstract, i.e. temporally/spatially local and requiring no more than $\mathcal{O}({n}^{2})$ memory, we have not yet discussed any explicit biological implementations. There are a handful of additional considerations for evaluating an algorithm with respect to biological plausibility:

i
Any equation describing synaptic strength changes (weight updates) must be local, i.e. any variables needed to update a synapse connecting the $i$th and $j$th units must be physically available to those units.

ii
Matrixvector multiplication can be implemented by networkwide neural transmission, but input vectors must represent firing rates (e.g. postactivations $\mathbf{a}$) and not membrane potentials (e.g. preactivations $\mathbf{h}$), since interneuron communication is mediated by spiking.

iii
Feedback weights used to calculate $\overline{\mathbf{c}}$ cannot be perfectly symmetric with ${\mathbf{W}}^{\text{out}}$, since there is no evidence for biological weight symmetry (see Lillicrap et al., 2016).

iv
Matrices (e.g. $\mathbf{J}$ or $\mathbf{A}$) must represent a set of synapses, whose strengths are determined by some local update.
With a few modifications, many of the presented algorithms can satisfy these requirements. We briefly illustrate one particular case with DNI, as shown in Marschall et al. (2019). To address (i), the result of the synthetic gradient operation ${\sum}_{l}{\stackrel{~}{a}}_{l}^{(t)}{A}_{li}$ can be stored in an electrically isolated neural compartment, in a manner similar to biological implementations of feedforward backpropagation (Guerguiev et al., 2017; Sacramento et al., 2018), to allow for local updates to ${W}_{ij}$. For (ii), simply pass the bootstrapped estimate of ${\overline{\mathbf{c}}}^{(t+1)}$ from Eq. (26) through the activation function $\varphi $ so that it represents a neural firing rate. For (iii), one can use fixed, random feedback weights ${\mathbf{W}}^{\text{fb}}$ instead of the output weights to calculate ${\overline{\mathbf{c}}}^{(t)}$, as in Lillicrap et al. (2016). And for (iv), one can train a set of weights ${\mathcal{J}}_{ij}$ online to approximate the Jacobian by performing SGD on ${L}_{J}^{(t)}={\parallel {a}_{i}^{(t)}{\sum}_{j}{\mathcal{J}}_{ij}{a}_{j}^{(t1)}\parallel}^{2}$, which encodes the error of the linear prediction of the next network state by ${\mathcal{J}}_{ij}$. The update rule manifests as
$$\mathrm{\Delta}{\mathcal{J}}_{ij}\propto \left({a}_{i}^{(t)}\sum _{{j}^{\prime}}{\mathcal{J}}_{i{j}^{\prime}}{a}_{{j}^{\prime}}^{(t1)}\right){a}_{j}^{(t1)},$$ 
essentially a “perceptron” learning rule, which is local and biologically realistic. Although this approximation brings no traditional computation speed benefits, it offers a plausible mechanism by which a neural circuit can access its own Jacobian for learning purposes. This technique could be applied to any other algorithm discussed in this paper. We refer to this altered version of DNI as DNI(b) in the experiments section.
5 Experiments
We run a number of experiments to empirically validate our categorizations and compare performance of the algorithms reviewed here, using two different synthetic tasks.
5.1 Setup
We implemented every algorithm presented here in a custom NumPybased Python module.^{2}^{2} 2 Link to public code repository to be included upon acceptance. In every simulation, we use gradient descent with a learning rate of ${10}^{4}$, the fastest learning rate for which all algorithms are able to converge to a stable steadystate performance. We restrict ourselves to using a batch size of 1, because, in an online setting, the network must learn as data arrive in real time. Most algorithms demand additional configuration decisions and hyperparameter choices: the truncation horizon $T$ (FBPTT), the initial values of the tensors $\mathbf{A}$ and $\mathbf{B}$ (all approximations), the initial values of the learned timescales ${\alpha}_{i}$ (KeRNL), the distribution from which $\bm{\nu}$ is sampled for stochastic updates (UORO, KFRTRL, RKFRTRL), and the learning rate for the numerical updates (KeRNL, DNI). For each algorithm, we independently optimize these choices by hand (see Appendix B for details).
We evaluate each algorithm’s ability to learn two different synthetic tasks: an additive dependencies task (Add) inspired by Pitis (2016) and a mimic target RNN task (Mimic). In both tasks, a stream of i.i.d. Bernoulli inputs ${\mathbf{x}}^{(t)}\in {\{0,1\}}^{{n}_{\text{in}}}$ is provided to the RNN. For Add, ${n}_{\text{in}}={n}_{\text{out}}=1$. The label ${y}^{*(t)}$ has a baseline value of 0.5 that increases (or decreases) by 0.5 (or 0.25) if $x=1$ at $t{t}_{1}$ (or $t{t}_{2}$), for specified lags ${t}_{1}$ and ${t}_{2}$. The longer the lags of the dependencies, the more difficult the task. We choose ${t}_{1}=6$ and ${t}_{2}=10$. In the Mimic task, the labels ${\mathbf{y}}^{*(t)}\in {\mathbb{R}}^{{n}_{\text{out}}}$ are determined by the outputs of a randomly generated, untrained target RNN that is fed the same input stream $\{{\mathbf{x}}^{({t}^{\prime})}:{t}^{\prime}\le t\}$. We use ${n}_{\text{in}}={n}_{\text{out}}=32$ in Mimic, chosen so that learning $\mathbf{W}$ is necessary for strong performance.
For each task, we consider a version on two different time scales: when the network is perfectly discrete ($\alpha =1$, see Eq. 9) and when the network update has some time continuity ($\alpha =0.5$). For the $\alpha =0.5$ case, the tasks are stretched over time by a factor of 2 to compensate, and the dependencies are reduced to ${t}_{1}=3$ and ${t}_{2}=5$ in the Add task to keep the difficulty roughly the same.
5.2 Add task: result and analysis
Fig. 4 shows the performance of each learning algorithm on the Add task, in both $\alpha =1$ and $\alpha =0.5$ conditions, including a fixed $\mathbf{W}$ algorithm that does not train $\mathbf{W}$ at all but does train ${\mathbf{W}}^{\text{out}}$, as a baseline. As expected, since they compute exact gradients up to truncation, RTRL and FBPTT perform best, although KFRTRL is a sufficiently robust approximation to virtually tie with them. RKFRTRL and UORO perform similarly and worse than KFRTRL does, as expected, since these approximations carry significantly more variance than KFRTRL. However, in the $\alpha =0.5$ condition, their performance is similar to that of KFRTRL.
KeRNL is theoretically a stronger approximation of RTRL than RFLO, because of its ability to learn optimal ${A}_{ki}$ and ${\alpha}_{i}$ whereas RFLO has fixed ${A}_{ki}={\delta}_{ki}$ and ${\alpha}_{i}=\alpha $. However, the numerical procedure for updating ${A}_{ki}$ and ${\alpha}_{i}$ depends on several configuation/hyperparameter choices. Despite significant effort in exploring this space, we are not able to get KeRNL to perform better than RFLO, suggesting that the procedure for training ${A}_{ki}$ and ${\alpha}_{i}$ does more harm than good. In the original paper, Roth et al. (2019) show promising results using RMSProp (Tieleman and Hinton, 2012) and batched training, so we suspect that the perturbationbased method for training ${A}_{ki}$ is simply too noisy for online learning.
DNI sits somewhere between the RFLOKeRNL cluster and the rest, with its biologically realistic sibling DNI(b) performing slightly worse than DNI, as to be expected, since it is an approximation on top of an approximation. As with KeRNL, DNI’s numerical update of ${A}_{li}$ introduces more hyperparameters and implementation choices, but there is a larger space of configurations in which the updates improve rather than hinder the algorithm’s performance.
5.3 Mimic task: result and analysis
For the Mimic task (Fig. 5), we see similar clustering of the algorithms but not in the same order. RTRL, FBPTT, and KFRTRL perform the best. UORO and RKFRTRL perform similarly to each other, but relatively worse than they did on Add. Conversely, DNI, DNI(b), KeRNL, and RFLO perform relatively better on this task than they did on Add. The information content that must be memorized for Add is relatively small (${n}_{\text{in}}=1$), but the time horizon is quite long. In Mimic, it is difficult to quantify the exact time horizon, but it is clear that the network must memorize much more information (${n}_{\text{in}}=32$). Thus perhaps UORO and RKFRTRL are effective at maintaining information over time, but the stochasticity in the updates places a limit on how much information can be retained.
5.4 Gradient similarity analysis
We conduct an indepth investigation looking beyond task accuracy by directly comparing the gradients (or approximate gradients) produced by each algorithm. Fig. 6 shows how a given pair of algorithms align on a timestepbytimestep basis for the Add and Mimic tasks. Each subplot is a histogram giving the distribution of normalized dot products
$$\mathrm{cos}\left({\theta}_{XY}^{(t)}\right)=\frac{{\mathrm{\Delta}}_{X}^{(t)}\mathbf{W}\cdot {\mathrm{\Delta}}_{Y}^{(t)}\mathbf{W}}{\parallel {\mathrm{\Delta}}_{X}^{(t)}\mathbf{W}\parallel \parallel {\mathrm{\Delta}}_{Y}^{(t)}\mathbf{W}\parallel}$$  (27) 
for (flattened) weight updates ${\mathrm{\Delta}}_{X}^{(t)}\mathbf{W}$ and ${\mathrm{\Delta}}_{Y}^{(t)}\mathbf{W}$ prescribed by algorithms $X$ and $Y$, respectively, at time $t$.
Figs. 6 a,c show qualitatively similar trends:

i
As shown directly in Figs. 6b,d, PF algorithms (UORO, KFRTRL, RKFRTRL, RFLO, KeRNL) align better with RTRL than with FBPTT, and vice versa for FF algorithms (DNI).

ii
The deterministic PF algorithms (RFLO and KeRNL) align better with RTRL than the stochastic algorithms (UORO, KFRTRL, and RKFRTRL) align with RTRL.

iii
RFLO and KeRNL align more strongly with each other than any other pair.

iv
UORO and RKFRTRL do not align strongly with any other algorithms, despite their ability to train the RNN effectively. UORO’s mean alignments with RTRL are 0.043 on Add and 0.084 on Mimic, while RKFRTRL’s mean alignments with RTRL are 0.050 on Add and 0.092 on Mimic (Fig. 6b,d), which are much lower than all other approximate algorithms, even those that perform worse on the task in some cases.
Observations (i)(iii) validate our categorizations, as similarity according to the normalized alignment corresponds to similarity by the pastfacing/futurefacing, tensor structure, and stochastic/deterministic criteria. Observation (iv) is puzzling, as it shows that angular alignment with exact algorithms is not predictive of learning performance.
How are UORO and RKFRTRL able to learn at all if their gradients are almost orthogonal with RTRL on average? We address this question for both UORO and RKFRTRL by examining the joint distribution of the gradient’s alignment with RTRL and the gradient’s norm (Fig. 7). All 4 cases show a statistically significant positive linear correlation between the normalized alignment and the common log of the gradient norm. This observation may partially explain (iv), because larger weight updates occur when UORO happens to strongly align with RTRL. However, these correlations are fairly weak even if statistically significant, and we argue that better algorithm similarity metrics are needed to account for observed differences in performance.
5.5 RFLO analysis
Among all approximate algorithms, RFLO stands out as having the simplest tensor structure and update rule, and it has been empirically demonstrated to be able to train RNNs on longterm dependencies. This is such a severe approximation of RTRL, yet it works so well in practice—and there is no clear understanding why. Although Murray (2019) goes into detail showing how loss decreases despite the random feedback used to approximately calculate ${\overline{\mathbf{c}}}^{(t)}$, he does not address the more basic mystery of how RFLO is able to learn despite the significant approximation ${\mathbf{J}}^{(t)}\approx (1\alpha )\mathbf{I}$. In this section, we provide some intuition for how this simple learning rule is so successful and empirically validate our claims.
We hypothesize that, rather than learning dynamics that actively retain useful bits of the past like RTRL and BPTT, RFLO works by training what is essentially a highcapacity feedforward model to predict labels from the natural memory traces of previous inputs contained in the hidden state. This is reminiscient of reservoir computing (Lukoševičius and Jaeger, 2009). We illustrate this idea in the special case of a perfectly discrete network ($\alpha =1$), where the learning rule still performs remarkably well despite ${B}_{ij}^{(t)}={\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$ containing no network history. As Fig. 8a depicts, ${\mathbf{a}}^{(t1)}$ ultimately maps to ${\mathbf{y}}^{(t)}$ via a singlehiddenlayer feedforward network parameterized by $\mathbf{W}$ and ${\mathbf{W}}^{\text{out}}$.
The RFLO learning rule in the discrete case corresponds exactly to training $\mathbf{W}$ by backpropagation:
$$\frac{\partial {L}^{(t)}}{\partial {W}_{ij}}=\frac{\partial {L}^{(t)}}{\partial {a}_{i}^{(t)}}\frac{\partial {a}_{i}^{(t)}}{\partial {W}_{ij}}={\overline{c}}_{i}^{(t)}{\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}={\overline{c}}_{i}^{(t)}{B}_{ij}^{(t)}.$$  (28) 
While every learning algorithm additionally trains ${\mathbf{W}}^{\text{out}}$ online to best map ${\mathbf{a}}^{(t)}$ to ${\mathbf{y}}^{*(t)}$, this purely linear model cannot perfectly capture the complex ways that information about past inputs ${\mathbf{x}}^{({t}^{\prime})}$, ${t}^{\prime}\le t$ implicit in ${\mathbf{a}}^{(t)}$ relates to labels ${\mathbf{y}}^{*(t)}$. Adding a hidden layer improves the ability of the network to predict ${\mathbf{y}}^{*(t)}$ from whatever evidence of ${\mathbf{x}}^{({t}^{\prime})}$ is naturally retained in ${\mathbf{a}}^{(t1)}$, analogous to how a singlehiddenlayer feedforward network outperforms a simple softmax regression on MNIST (Deng, 2012).
To empirically validate our explanation, we first show that the strength of natural memory traces in the RNN depends on its recurrent weights. We measure this “memory” by running an untrained RNN with fixed weights forwards for $20$k time steps of the Add task and calculating the ${r}^{2}$ value of a linear regression of ${\mathbf{a}}^{(t+\mathrm{\Delta}t)}$ onto ${\mathbf{y}}^{*(t)}$ for different values of the time shift $\mathrm{\Delta}t$. The sudden jumps in information occur as $\mathrm{\Delta}t$ passes the time lags of the inputoutput dependencies explicitly included in the task (${t}_{1}=6,{t}_{2}=10$), followed by a slow decay as the information relevant for predicting ${\mathbf{y}}^{*(t)}$ gets corrupted by running the network forwards. The speed of this decay differs by the choice of $\mathbf{W}$. As Fig. 8b shows, orthogonal and Gaussian $\mathbf{W}$ seem to best preserve information over time, while symmetric and diagonal $\mathbf{W}$ lose information quite rapidly, likely due to their real spectra. In separate simulations (Fig. 8c), we trained networks initialized in these ways using the RFLO rule, and only the networks initialized with orthogonal or Gaussian $\mathbf{W}$ are able to learn at all with RFLO. This validates our hypothesis, that RFLO works by a static prediction of ${\mathbf{y}}^{*(t)}$ based on evidence in ${\mathbf{a}}^{(t1)}$, because in cases where this evidence is absent (or at least weak) due to rapid decay, RFLO fails.
6 Discussion
We have presented a framework for conceptually unifying many recent advances in efficient online training of recurrent neural networks. These advances come from multiple perspectives: pure machine learning, biologically plausible machine learning, and computational neuroscience. We started by articulating a crucial distinction, (a) past facing vs. future facing, that divides algorithms by what exactly they calculate. We then presented a few other properties that characterize each new technique: (b) the tensor structure of the algorithm, (c) whether its update requires explicit randomness, and (d) whether its update can be derived in closed form from the relations (5)–(7) or must be approximated with SGD or some other numerical approach. Along the way, we clarified the relationship between some of the modern online approximations and exact methods. Further, we showed it’s possible to create new algorithms by identifying unexplored combinations of properties (a)(d).
We empirically validated these ideas with synthetic tasks on which all algorithms could perform reasonably well in an online setting, despite using vanilla RNN architecture; gradient descent optimization with fixed learning rate; and no standard machine learning techniques such as gradient clipping, batch/layer normalization, ${L}^{2}$ regularization, etc. We saw that training errors roughly cluster according to their categorical distinctions, across tasks and hyperparameter choices. But performance of these clusters differed across tasks: for example, UORO and RKFRTRL performed relatively well (poorly) on Add (Mimic), while for DNI and DNI(b) this effect was flipped. We speculate that this can be explained in terms of the task demands: the Add task has a long time horizon (10 time steps) of explicit dependencies for small inputs (${n}_{\text{in}}=1$), while the Mimic task requires more information to be memorized (${n}_{\text{in}}=32$) that decays with the target RNN’s forward pass (harder to measure, but likely shorter than 10 time steps). Thus perhaps UORO and RKFRTRL produce gradients with a limited amount of information that survives many updates, while DNI and DNI(b) have a larger information capacity but a limited time horizon. Future work should more systematically explore how these different task features interact with the different algorithmic approximations.
Following common practice, we used the pairwise vector alignments of the (approximate) gradients calculated by each algorithm as a way to analyze the precision of different approximations. This similarity turned out to reflect the natural clustering of the algorithms along the axes proposed here. In particular, pastfacing approximations had stronger alignment with the exact pastfacing gradients calculated by RTRL compared to the exact (up to truncation) futurefacing gradients calculated by FBPTT, and vice versa for the futurefacing approximation, DNI. Reassuringly, KeRNL and RFLO, which have the same tensor structure, featured strong alignment.
Importantly, the angular alignment of the gradients did not account for task performance: UORO and RKFRTRL performed quite well despite their weak alignment with RTRL and BPTT, while KeRNL performed relatively poorly despite its strong alignment with RTRL and BPTT. Analyzing the magnitudes of the gradients partially explained this observation, as UORO and RKFRTRL aligned with RTRL more strongly when their gradients were larger in norm. However, this effect was subtle, so probably not enough to account for the performance differences. This exposes a limitation of current ways of analyzing this class of algorithms. There clearly is a need for better similarity metrics that go beyond timepointwise gradient alignment and instead compare longterm learning trajectories.
Notably, it is the stochastic algorithms that have particularly weak alignment with RTRL—even KFRTRL, yet its performance is indistinguishable from the exact algorithms. Averaged over many time steps of learning, corresponding to many samples of $\bm{\nu}$, these stochastic methods do contain complete information about ${M}_{kij}^{(t)}$ (i.e. are unbiased), but at any one time point the alignment is heavily corrupted by the explicit randomness. In contrast, deterministic approximations, such as KeRNL, RFLO and DNI, may partially align with exact methods by construction, but their errors have no reason to average out, hence their inability to find the same minima as exact methods (Fig. 9). This may also explain why stochastic approximations do not align with each other despite their conceptual similarity.
We chose the simplest possible setup for comparison, using basic RNN architecture and gradient descent, aiming to be as inclusive as possible of various algorithms. Of course, there are alternative online methods that do not naturally fall into this framework. Any online algorithm that exploits a particular RNN architecture necessarily cannot fit into our classification. For example, Ororbia et al. (2017, 2018) propose a specialized neural architecture (Temporal Neural Coding Network), whose learning algorithm (Discrepancy Reduction) is naturally local, due to the network structure. In contrast, the algorithms reviewed here were explicitly derived as approximations to nonlocal algorithms (RTRL and BPTT), and their locality ends up being more of a bug than a feature, as evidenced by their impaired performance. From a machine learning perspective, there are limits to generalpurpose approximations of this kind, and future progress in online methods will likely come about by further exploring specialized algorithmarchitecture pairs.
The same is true from a biological perspective. In the brain, cortical architecture and synaptic plasticity rules evolved together, while physical constraints dictate that plasticity is necessarily local. Still, the details of how local plasticity rules interact with neural circuits remain mysterious, and this is a current focus of research (Guerguiev et al., 2017; Sacramento et al., 2018; Lillicrap and Santoro, 2019). Exploring which architectures allow locality to manifest as a consequence, rather than a constraint, of learning is a potentially fruitful point of interaction between artificial intelligence and computational neuroscience.
Acknowledgments
We thank Guangyu Yang and James Murray for helpful conversations.
References
 Bahdanau et al. (2014) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. Neural machine translation by jointly learning to align and translate. arXiv preprint arXiv:1409.0473, 2014.
 Bahdanau et al. (2016) Dzmitry Bahdanau, Jan Chorowski, Dmitriy Serdyuk, Philemon Brakel, and Yoshua Bengio. Endtoend attentionbased large vocabulary speech recognition. In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4945–4949. IEEE, 2016.
 Benzing et al. (2019) Frederik Benzing, Marcelo Matheus Gauy, Asier Mujika, Anders Martinsson, and Angelika Steger. Optimal kroneckersum approximation of real time recurrent learning. arXiv preprint arXiv:1902.03993, 2019.
 Cho et al. (2014) Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoderdecoder for statistical machine translation. arXiv preprint arXiv:1406.1078, 2014.
 Cho et al. (2015) Kyunghyun Cho, Aaron Courville, and Yoshua Bengio. Describing multimedia content using attentionbased encoderdecoder networks. IEEE Transactions on Multimedia, 17(11):1875–1886, 2015.
 Cooijmans and Martens (2019) Tim Cooijmans and James Martens. On the variance of unbiased online recurrent optimization. arXiv preprint arXiv:1902.02405, 2019.
 Czarnecki et al. (2017) Wojciech Marian Czarnecki, Grzegorz Swirszcz, Max Jaderberg, Simon Osindero, Oriol Vinyals, and Koray Kavukcuoglu. Understanding synthetic gradients and decoupled neural interfaces. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 904–912. JMLR. org, 2017.
 Deng (2012) Li Deng. The mnist database of handwritten digit images for machine learning research [best of the web]. IEEE Signal Processing Magazine, 29(6):141–142, 2012.
 Graves (2013) Alex Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:1308.0850, 2013.
 Graves et al. (2016) Alex Graves, Greg Wayne, Malcolm Reynolds, Tim Harley, Ivo Danihelka, Agnieszka GrabskaBarwińska, Sergio Gómez Colmenarejo, Edward Grefenstette, Tiago Ramalho, John Agapiou, et al. Hybrid computing using a neural network with dynamic external memory. Nature, 538(7626):471, 2016.
 Guerguiev et al. (2017) Jordan Guerguiev, Timothy P Lillicrap, and Blake A Richards. Towards deep learning with segregated dendrites. ELife, 6:e22901, 2017.
 Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. Long shortterm memory. Neural computation, 9(8):1735–1780, 1997.
 Jaderberg et al. (2017) Max Jaderberg, Wojciech Marian Czarnecki, Simon Osindero, Oriol Vinyals, Alex Graves, David Silver, and Koray Kavukcuoglu. Decoupled neural interfaces using synthetic gradients. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 1627–1635. JMLR. org, 2017.
 Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
 Lillicrap and Santoro (2019) Timothy P Lillicrap and Adam Santoro. Backpropagation through time and the brain. Current Opinion in Neurobiology, 55:82 – 89, 2019. ISSN 09594388. doi: https://doi.org/10.1016/j.conb.2019.01.011. URL http://www.sciencedirect.com/science/article/pii/S0959438818302009. Machine Learning, Big Data, and Neuroscience.
 Lillicrap et al. (2016) Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman. Random synaptic feedback weights support error backpropagation for deep learning. Nature communications, 7:13276, 2016.
 Lukoševičius and Jaeger (2009) Mantas Lukoševičius and Herbert Jaeger. Reservoir computing approaches to recurrent neural network training. Computer Science Review, 3(3):127–149, 2009.
 Marschall et al. (2019) Owen Marschall, Kyunghyun Cho, and Cristina Savin. Using local plasticity rules to train recurrent neural networks. In Computational and Systems Neuroscience (CoSyNe), number II90, Lisbon, Portugal, Feb 2019.
 Mikolov et al. (2010) Tomáš Mikolov, Martin Karafiát, Lukáš Burget, Jan Černockỳ, and Sanjeev Khudanpur. Recurrent neural network based language model. In Eleventh annual conference of the international speech communication association, 2010.
 Mnih et al. (2015) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. Humanlevel control through deep reinforcement learning. Nature, 518(7540):529, 2015.
 Mujika et al. (2018) Asier Mujika, Florian Meier, and Angelika Steger. Approximating realtime recurrent learning with random kronecker factors. In Advances in Neural Information Processing Systems, pages 6594–6603, 2018.
 Murray (2019) James M Murray. Local online learning in recurrent networks with random feedback. eLife, 8:e43299, 2019.
 Ororbia et al. (2018) Alexander Ororbia, Ankur Mali, C Lee Giles, and Daniel Kifer. Online learning of recurrent neural architectures by locally aligning distributed representations. arXiv preprint arXiv:1810.07411, 2018.
 Ororbia et al. (2017) II Ororbia, G Alexander, Patrick Haffner, David Reitter, and C Lee Giles. Learning to adapt by minimizing discrepancy. arXiv preprint arXiv:1711.11542, 2017.
 Pitis (2016) Silviu Pitis. Recurrent neural networks in tensorflow 1. r2rt.com/recurrentneuralnetworksintensorflowi, 2016. Accessed: 20181113.
 Roth et al. (2019) Christopher Roth, Ingmar Kanitscheider, and Ila Fiete. Kernel RNN learning (keRNL). In International Conference on Learning Representations, 2019. URL https://openreview.net/forum?id=ryGfnoC5KQ.
 Sacramento et al. (2018) João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn. Dendritic cortical microcircuits approximate the backpropagation algorithm. In Advances in Neural Information Processing Systems, pages 8721–8732, 2018.
 Sutton and Barto (2018) Richard S Sutton and Andrew G Barto. Reinforcement learning: An introduction. MIT press, 2018.
 Tallec and Ollivier (2017) Corentin Tallec and Yann Ollivier. Unbiased online recurrent optimization. arXiv preprint arXiv:1702.05043, 2017.
 Tieleman and Hinton (2012) T. Tieleman and G. Hinton. Lecture 6.5—RmsProp: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural Networks for Machine Learning, 2012.
 Werbos et al. (1990) Paul J Werbos et al. Backpropagation through time: what it does and how to do it. Proceedings of the IEEE, 78(10):1550–1560, 1990.
 Williams and Zipser (1989) Ronald J Williams and David Zipser. A learning algorithm for continually running fully recurrent neural networks. Neural computation, 1(2):270–280, 1989.
Appendix A.
Lemma for generating rank1 unbiased estimates
For completeness, we state the Lemma from Tallec and Ollivier (2017) in components notation. Given a decomposition of a matrix $\mathbf{M}\in {\mathbb{R}}^{n\times m}$ into $r$ rank1 components
$${M}_{ij}=\sum _{k=1}^{r}{A}_{ik}{B}_{kj},$$  (29) 
a vector of i.i.d. random variables $\bm{\nu}\in {\mathbb{R}}^{r}$ with $\mathbb{E}[{\nu}_{k}]=1$, $\mathbb{E}[{\nu}_{k}{\nu}_{{k}^{\prime}}]={\delta}_{k{k}^{\prime}}$, and a list of $r$ positive constants ${\rho}_{k}>0$, then
$${\stackrel{~}{M}}_{ij}=\left(\sum _{k=1}^{r}{\rho}_{k}{\nu}_{k}{A}_{ik}\right)\left(\sum _{k=1}^{r}{\rho}_{k}^{1}{\nu}_{k}{B}_{kj}\right)$$  (30) 
is a rank1, unbiased estimate of ${M}_{ij}$ over the choice of $\bm{\nu}$.
Appendix B.
Implementation details
For reproducibility, we describe in fuller detail our implementation configurations for each simulation. Table 2 shows hyperparameter/configuration choices that apply across all algorithms. Table 3 shows the algorithmspecific hyperparameter choices we made for each task. In Table 2, we reference submatrices of $\mathbf{W}=[{\mathbf{w}}^{\text{rec}},{\mathbf{w}}^{\text{in}},{\mathbf{b}}^{\text{rec}}]$ and ${\mathbf{W}}^{\text{out}}=[{\mathbf{w}}^{\text{out}},{\mathbf{b}}^{\text{rec}}]$, since they are initialized differently.
hyperparameter  value  explanation 

learning rate  ${10}^{4}$  learning rate for SGD w.r.t. $\mathbf{W}$ and ${\mathbf{W}}^{\text{out}}$ 
$n$  32  number of hidden units in the network 
$\varphi $  tanh  nonlinearity used in RNN forward dynamics 
init. ${\mathbf{w}}^{\text{in}}$  $\sim \mathcal{N}(0,1/\sqrt{{n}_{\text{in}}})$  initial value for input weights 
init. ${\mathbf{w}}^{\text{rec}}$  rand. orth.  initial value for recurrent weights 
init. ${\mathbf{b}}^{\text{rec}}$  0  initial value for recurrent bias 
init. ${\mathbf{w}}^{\text{out}}$  $\sim \mathcal{N}(0,1/\sqrt{n})$  initial value for output weights 
init. ${\mathbf{b}}^{\text{out}}$  0  initial value for output bias 
init. ${\mathbf{W}}^{\text{FB}}$  $\sim \mathcal{N}(0,1/\sqrt{{n}_{\text{out}}})$  value for fixed feedback weights used in DNI(b) 
init. ${\mathbf{b}}_{\text{targ.}}^{\text{rec}}$  $\sim \mathcal{N}(0,0.1)$  initial value for target recurrent bias in Mimic 
init. ${\mathbf{b}}_{\text{targ.}}^{\text{out}}$  $\sim \mathcal{N}(0,0.1)$  initial value for target output bias in Mimic 
algorithm  initial values  $\bm{\nu}$ dist.  LR  pert. $\sigma $  $T$ 
UORO  ${A}_{k}^{(0)}\sim \mathcal{N}(0,1),{B}_{ij}^{(0)}\sim \mathcal{N}(0,1)$  unif. $\{1,1\}$  
KFRTRL  ${A}_{j}^{(0)}\sim \mathcal{N}(0,1),{B}_{ki}^{(0)}\sim \mathcal{N}(0,1/\sqrt{n})$  unif. $\{1,1\}$  
RKFRTRL  ${A}_{i}^{(0)}\sim \mathcal{N}(0,1),{B}_{kj}^{(0)}\sim \mathcal{N}(0,1/\sqrt{n})$  unif. $\{1,1\}$  
KeRNL  ${A}_{ki}^{(0)}={\delta}_{ki},{B}_{ij}^{(0)}=0,{\alpha}_{i}^{(0)}=0.8$  $5$  ${10}^{3}$  
DNI  ${A}_{li}^{(0)}\sim \mathcal{N}(0,1/\sqrt{n})$  ${10}^{3}$  
DNI(b)  ${A}_{li}^{(0)}\sim \mathcal{N}(0,1/\sqrt{n}),{\mathcal{J}}_{ij}^{(0)}={W}_{ij}^{\text{rec}}$  ${10}^{3}$  
FBPTT  10 
Some miscellaneous implementation details below:

•
For the Add task in the $\alpha =1$ condition, we changed the DNI/DNI(b) learning rate to $5\times {10}^{2}$ for ${A}_{li}$ and ${10}^{2}$ for ${\mathcal{J}}_{ij}$ (DNI(b)). In other cases, the learning rates for ${A}_{li}$ and ${\mathcal{J}}_{ij}$ are identical.

•
There are two appearances of the synthetic gradient weights ${A}_{li}$ in Eq. (26). Although we wrote them as one matrix $\mathbf{A}$ for brevity, in implementation we actually keep two separate values, $\mathbf{A}$ and ${\mathbf{A}}^{*}$, the latter of which we use for for the righthand appearance ${A}_{{l}^{\prime}m}$ (specifically to calculate the bootstrapped estimate of the SG training label). We update $\mathbf{A}$ every time step but keep ${\mathbf{A}}^{*}$ constant, replacing it with the latest value of $\mathbf{A}$ only once per $\tau \in \mathbb{N}$ time steps. This integer $\tau $ introduces another hyperparameter, which we choose to be 5. (Inspired by an analogous technique used in deep Qlearning from Mnih et al., 2015.)

•
In the original paper, Roth et al. (2019) use $(1\mathrm{exp}({\gamma}_{i}))$ rather than ${\alpha}_{i}$ as a temporal filter for ${B}_{ij}^{(t)}$. We made this change so that ${\alpha}_{i}$ makes sense in terms of the $\alpha $ in the forward dynamics of the network and RFLO. Of course, these are equivalent via ${\gamma}_{i}=\mathrm{log}(1{\alpha}_{i})$, but the gradient w.r.t. ${\alpha}_{i}$ must be rescaled by a factor of $1/(1{\alpha}_{i})$ to compensate.

•
For KeRNL, there is a choice for how to update the eligibility trace (Eq. 17): one can scale the righthand term ${\varphi}^{\prime}({h}_{i}^{(t)}){\widehat{a}}_{j}^{(t1)}$ by either the learned timescale ${\alpha}_{i}$ or the RNN timescale $\alpha $. We chose the latter because it has stronger empirical performance and it theoretically recovers the RTRL equation under the approximating assumptions about $\mathbf{A}$.

•
Perturbations for calculating gradients for ${A}_{ki}$ and ${\alpha}_{i}$ in KeRNL are sampled i.i.d. ${\zeta}_{i}\sim \mathcal{N}(0,\sigma )$.

•
In our implementation of the Add task, we use ${n}_{\text{in}}={n}_{\text{out}}=2$ for a “onehot” representation of the input ${x}^{(t)}\in \{0,1\}$ and label ${y}^{*(t)}\in \{0.25,0.5,0.75,1\}$, such that ${\mathbf{x}}^{(t)}=[{x}^{(t)},1{x}^{(t)}]$ and ${\mathbf{y}}^{*(t)}=[{y}^{(t)},1{y}^{(t)}]$.

•
In our implementation of Mimic, the target RNN was initialized in the same way as the RNNs we train, with the exception of the recurrent and output biases (see Table 2).