Ensemble perspective for understanding temporal credit assignment

Recurrent neural networks are widely used for modeling spatio-temporal sequences in both nature language processing and neural population dynamics. However, understanding the temporal credit assignment is hard. Here, we propose that each individual connection in the recurrent computation is modeled by a spike and slab distribution, rather than a precise weight value. We then derive the mean-field algorithm to train the network at the ensemble level. The method is then applied to classify handwritten digits when pixels are read in sequence, and to the multisensory integration task that is a fundamental cognitive function of animals. Our model reveals important connections that determine the overall performance of the network. The model also shows how spatio-temporal information is processed through the hyperparameters of the distribution, and moreover reveals distinct types of emergent neural selectivity. It is thus promising to study the temporal credit assignment in recurrent neural networks from the ensemble perspective.

FIG. 1: Illustration of a RNN learning the temporal credit assignment. Each connection is described by the spike and slab distribution; the spike indicates the corresponding synaptic weight absent for a task, while the slab represents a Gaussian distribution responsible for the weight uncertainty. The Gaussian slab is displayed in the plot, and the arrow in the reservoir shows different directions of information flow. Different sizes of connection show how probable these connections should be present during recurrent computation. moreover the selectivity of each recurrent unit, we are able to provide mechanistic understanding about the recurrent dynamics in both engineering applications and computational neuroscience modeling. Our method learns the statistics of the RNN ensemble, producing a dynamic architecture that adapts to temporally varying inputs, which thereby goes beyond traditional training protocols focusing on a stationary network topology. The hyperparameters governing the weight distribution reveal which credit assignments are critical to the network behavior, which further explains distinct functions of each computational layer and the emergent neural selectivity. Therefore, our ensemble theory can be used as a promising tool to explore internal dynamics of widely used RNNs.

II. MODEL AND ENSEMBLE TRAINING
In this study, we consider a recurrent neural network processing a time-dependent input x(t) of time length T . The input signal is sent to the recurrent reservoir via the input weight matrix W in . The number of input units are determined by design details of the task, and w in ij denotes the weight value for the connection from input unit j to recurrent node i. The neural responses of recurrent units are represented by an activity vector r(t) at time step t. We define w ij as the connection from node j to node i in the reservoir. In general, w ij = w ji , indicating different weights for different directions of information flow. In addition, we do not preclude the self-interaction w ii , which can be used to maintain the representation encoded in the history of the internal dynamics [19]. The specific statistics of this self-interaction can be determined by the learning shown below. The internal activity r(t) is read out through the output matrix W out in the form of the time-dependent output y(t).
We first define h i (t) as the time-dependent synaptic current of neuron i. The dynamics of the recurrent network can thus be summarized as follows, where u(·) is the pre-activation function, φ(·) denotes the nonlinear transfer function, for which we use the rectified linear unit (ReLU) function for all tasks. α = ∆t τ , where ∆t is a small time interval (e.g., emerging from discretization of a continuous dynamics [20]), and τ denotes the time constant of the dynamics. n i ∼ N (0, 1) indicates a normally distributed random number with zero mean and unit variance, sampled independently at every time step, and σ controls the strength of the recurrent noise intrinsic to the network. This intrinsic noise is present in modeling cognitive tasks, but absent for engineering applications. Considering a Gaussian white noise and a baseline (x 0 ), we write the input vector x(t) as where the input has been written in a discrete time form, x task t denotes the sequence of the task, and ξ i ∼ N (0, 1), which is independently sampled at each time step. The noise term is also called the external sensory noise of strength σ in , commonly observed in brain circuits [21,22], but is absent in engineering applications. For a MNIST classification task, f (·) is chosen to be the softmax function, specifying the probability over all the classes at the last time step T , i.e., y k (T ) = e z k (T ) j e z j (T ) . We defineŷ k as the target label (one-hot form), and use the cross entropy L = − iŷ i ln y i (T ) as the loss function to be minimized. For a multisensory integration task, f (·) is an identity function, and the mean squared-error is chosen to be the objective function.
To search for the optimal random network ensemble for time-dependent tasks in recurrent neural networks, we model the statistics of the weight matrices by the SaS probability distribution [13] as follows, The spike mass at δ(·) is related to the network compression, indicating the necessary weight resources required for a specific task. The continuous slab, N (w ij |m ij , Ξ ij ), denotes the Gaussian distribution with mean m ij and variance Ξ ij , characterizing the weight uncertainty when the corresponding connection can not be absent (see Fig. 1). Next, we derive the learning equations about how the SaS parameters are updated, based on mean-field approximation. More precisely, we consider the average over the statistics of the network ensemble during training. Notice that, the first and second moments of the weight w ij for three sets of weights can be written in a common form as . Given a large fan-in, the pre-activation u i (t) and the output z i (t) can be re-parametrized by using standard Gaussian random variables , i.e., it is reasonable to assume that they are subject to N (u i (t)|G in , respectively, according to the central-limit-theorem. Note that when the number of fan-in to a recurrent unit is small, the centrallimit-theorem may break. In this situation, we take the deterministic limit. Therefore, the mean-field dynamics of the model becomes where Note that { u (t)} and { out (t)} are both independent random variables sampled from the standard Gaussian distribution with zero mean and unit variance, which are quenched for every single training mini-epoch and also time-step dependent, maintaining the same sequence of values in both feedforward and backward computations. Updating the network parameters (θ θ θ in ik , θ θ θ ik , θ θ θ out ki ) can be achieved by the gradient descent on the objective function L. First of all, we update θ θ θ out ki ≡ (m out ki , π out ki , Ξ out ki ).
where ∂L ∂z k (t) is related to the form of the loss function. For categorization tasks, f (·) is chosen to be the softmax function, and we use the cross entropy as our objective function, for which ∂L ∂z k (t) = (y k (T ) −ŷ k )δ t,T . It is worth noticing that for multi-sensory integration tasks, this derivative does not vanish at intermediate time steps and become thus time-dependent. The other term ∂z k (t) ∂θ out ki can be directly computed, showing how sensitive the network activity is read out under the change of the SaS parameters in the output layer: (7c) Next, we derive the learning equation for the hyper-parameters in the reservoir and input layer, i.e., θ θ θ ij ≡ (m ij , π ij , Ξ ij ) and θ θ θ in ij ≡ (m in ij , π in ij , Ξ in ij ). To get a general form, we setm ij ≡ (m ij , m in ij ),π ij ≡ (π ij , π in ij ), and Ξ ij ≡ (Ξ ij , Ξ in ij ). Based on the chain rule, we then arrive at the following equations, where we have defined δ i (t) ≡ ∂L ∂hi(t) . The auxiliary variable δ i (t) can be computed by the chain rule once again, resulting in a BPTT equation of the error signal starting from the last time-step T : where t = 0, 1, 2, ..., T − 1, and φ (·) denotes the derivative of the transfer function. The first summation in δ i (t) can be directly expanded as j ∂L ∂h j (t + 1) .
The second summation in δ i (t) is given by It is worth noting that for the last time-step T , the error signal δ i (T ) is written as To compute Eq. (8), we need to work out the following derivatives, which characterize the sensitivity of the preactivation under the change of hyper-parameters θ θ θ ij and θ θ θ in ij . We summarize the results as follows, In this learning process, our model learns a RNN ensemble to realize the time-dependent computation, in contrast to the standard BPTT algorithm which gives only a point-estimate of RNN weights. In particular, if we set π = 0 and Ξ = 0, our learning equation reduces to the standard BPTT. Hence, our model can be thought of as a generalized version of BPTT (i.e., gBPTT), as shown in Fig. 2, where each weight parameter should be understood as the SaS hyper-parameters, corresponding to our ensemble setting.
The size of the candidate network space can be captured by the network entropy S = − R D P (w) ln P (w)dw, where D is the number of weight parameters in the network. If we assume the joint distribution of weights to be factorized across individual connections, the overall energy S can be obtained by summing up the entropy of individual weights as S = S . The entropy of each directed connection is derived as follows [13]: where Γ( s ) = ln [π δ(m + √ Ξ s ) + 1−π √ Ξ N ( s |0, 1)] and B denotes the number of standard Gaussian variables s . This model entropy is just an approximate estimate of the true value whose exact computation is impossible.
If π = 0, the probability distribution of w can be written as P (w ) = N (w |m , Ξ ), and the entropy S can be analytically computed as 1 2 ln (2πeΞ ). If Ξ = 0, the Gaussian distribution reduces to a Dirac delta function, and the entropy becomes an entropy of discrete random variables, S = π ln π + (1 − π ) ln (1 − π ). However, there may exist a mixture of discrete and continuous contributions to the entropy. Therefore, the entropy value can take negative values. In practice, we use δ( 2a to approximate the delta-peak with a small value of a. The stochasticity of the SaS distribution can be also decomposed into two levels: one at the choice of connection that is absent and the other at the Gaussian slab itself. The former one is a discrete entropy characterized by the spike probability, and the latter is a continuous entropy just characterized by the variance of the Gaussian slab.

III. RESULTS
In this section, we show applications of our ensemble theory of temporal credit assignment to both engineering and computational cognitive tasks, i.e., pixel-by-pixel MNIST digit classification, and the multisensory integration demonstrating the benefit of multiple sources of information for decision making. We will explore in detail rich properties of trained RNN model accomplishing the above computational tasks of different nature.
Illustration of the generalized back-propagation through time. Xt is a time-dependent input into the network, resulting in a sequence of hidden representations r(t) (i.e., φ(ht)) and a time-dependent loss Lt. The error signal ∂L t ∂h t propagates back from the last time-step T to the first time step, yielding the temporally accumulated gradients for hyper-parameters to be updated. The contribution of ∂L t ∂z t is not shown in this plot. Xt, zt and ht indicate x(t), z(t) and h(t) in the main text, respectively.
FIG. 3: Properties of trained RNN models for the pixel-by-pixel MNIST classification. (a) Test accuracy versus training epoch. The network of N = 100 neurons is trained on the full MNIST dataset (60 000 images), and tested on the unseen data of 10 000 images. Five independent runs of the algorithm are considered. There are no baseline inputs. Other parameters for training are as follows: α = 0.1, 2 regularization strength is 10 −4 , and the initial learning rate lr0 = 0.001. (b) Evolution of the network sparsity per connection. Recurrent I denotes the case of wij (i < j) along which the recurrent feedback passes from neuron j to neuron i; while recurrent II denotes the case of wji (j > i) along which the recurrent feedback passes from neuron i to neuron j. The fluctuation is estimated from five independent runs [the same training conditions as in (a)].

A. Pixel-by-Pixel MNIST digit classification
Training RNNs is hard because of long-term dependency in the sequence of inputs. A challenging task of long-term dependency is to train a RNN to classify the MNIST images when the 784 pixels are fed into the network one by one (note that x 0 = 0), where the network is required to predict the category of the image after seeing all the 784 pixels. This task displays a long range of dependency, because the network reads one pixel at a single time-step in a scan-line order from the top left pixel to the bottom right pixel (Fig. 5), and the information of as long as 784 time steps must be maintained before the final decision. We apply the vanilla RNN with N = 100 recurrent units to achieve this challenging goal. Because the input size is one in this task, and thus the set of hyper-parameters (π in , Ξ in ) are set to zero (i.e., the deterministic limit). The entire MNIST dataset is divided into mini-batches for stochastic gradient descent, and we apply cross entropy as the objective function to be minimized. In this task, the  Fig. 3. Five independent runs are considered. The input weight is deterministic without any stochasticity. In addition, to avoid divergence, we relax all values of Ξ smaller than a to be a in the computation of the entropy. Our results are not sensitive to the choice of a, once the value of a is small. (a,b) We decompose the variability of the SaS distribution into two levels: the variability of selecting zero synapses (π-entropy) and the variability inside the Gaussian slab (Ξ-entropy). (c,d) Evolution of entropy per connection (B = 100). We use the Gaussian distribution of small variance (indicated by a) to approximate the delta-peak. The value of a does not affect the qualitative behavior of the entropy profile.
error signal appears only after the network reads all the 784 pixels. In the current setting, we ignore the noise terms in the dynamics equations [Eq. (1) and Eq. (2)]. Surprisingly, although working at the ensemble level, our model can achieve a comparable or even better performance than the traditional BPTT method, as shown in Fig. 3 (a).
Our simulation reveals that the sparsity per connection, l π l M (a total of M directed connections in the network), achieves a larger value in the recurrent layer compared with the output layer; this behavior does not depend on the specific directions of the coupling [ Fig. 3 (b)]. During training, the sparsity level grows until saturation, suggesting that the training may remove or minimize the impacts of irrelevant information across time steps. Interestingly, the entropy per connection in the recurrent layer also increases at the early stage of the training, but decreases at the late stage [ Fig. 4], showing that the training is able to reorganize the information landscape shaped by the recurrent feedbacks. The network at the end of training becomes more deterministic (e.g., Ξ gets close to zero, or the spike probability goes away from one half). The discrete π-entropy always increases until saturation at a certain level. The entropy profile in the output layer shows a similar behavior yet with a lower entropy value. This is consistent with distinct roles of recurrent and readout layers. The recurrent layer (or computing reservoir) searches for an optimal way of processing temporally complex sensory inputs, while the output layer implements the decoding of the signals hidden in the reservoir dynamics. Note that two types of deterministic weights yield vanishing individual entropy values, i.e., (i) π= 1 (unimportant (UIP) weight); (ii) π = 0, Ξ = 0 and m = 0 (very important (VIP) weight).
To study the network behavior from the perspective of hyper-parameter distributions, we plot the distribution of three sets of parameters (m, Ξ, π) for the output layer, two different directions of connections in the recurrent layer, and m of the input layer ((π in , Ξ in ) are set to zero as explained before), as shown in Fig. 5. The distribution of spike probability π has the shape of L for all layers. The extreme at π = 0 indicates that the corresponding synaptic weight carries important feature-selective information, and are thus critical to the network performance. The shape of π-distribution profile reflects the task difficulty for the considered network (given the same initialization condition), since we observe that in a simpler task for the network to read 28 pixels at each time step (rather than one pixel by one pixel), the profile of π-distribution develops a U shape, with the other peak at π = 1, implying the emergence of a significant fraction of UIP weights. In addition, Ξ has an L-shaped distribution which peaks at zero, suggesting that the corresponding weight distribution takes a deterministic value of m, which becomes the weight value of that connection.
As for the readout layer, the mean of the continuous slab becomes much more dispersed, compared with that of the recurrent layer, while the statistics profile for π and Ξ becomes more converged, which is in an excellent agreement with the entropy profile shown in Fig. 4. We thus conclude that the recurrent layer has a greater variability than the output layer. This great variability allows the recurrent layer to transform the complex sensory inputs with hierarchical spatiotemporal structures in a flexible way. Nevertheless, the minimal variability makes the readout behavior more robust. It is therefore interesting in future works to address the precise mechanism underlying how the statistics of the connections in a RNN facilitates the formation of latent dynamics for computation.
We next look at the specific profile of individual hyper-parameters in a matrix form (Fig. 6). We find that the diagonal (self-interaction) of m emerges from gBPTT (note that the diagonals of π and Ξ nearly vanish), demonstrating the significant role of self-interactions in maintaining long-term memory and thereby the learning performance, in accord with the heuristic strategy used in a recent study [19]. There also appears heterogeneity in the hyper-parameter matrices, i.e., some directed connections play a more important role than the others. In particular, some connections can be eliminated for saving computation resources.
Our method can identify the exact nature of each directed connection in the network. Targeted weight perturbation can thus be performed in the recurrent layer (Fig. 7) . In the course of training, the fraction of VIP connections slowly decreases, although the fraction of VIP connections is not significant (around 1%) [ Fig. 7(b)]. Surprisingly, pruning these VIP weights could strongly deteriorate the test performance of the network (up to the chance level), whereas pruning the same number of randomly selected connections from the entire weight population yields a negligible drop of the test accuracy. We thus conclude that in a RNN, there exists a minority of VIP weights carrying key spatiotemporal information essential for decision making of the network.
We then explore whether the selectivity of neurons for the sensory inputs could emerge from our training. The degree of selectivity for an individual unit, say neuron i, can be characterized by an index, namely SLI [23] as follows where r i,s indicates the response of the neuron to the input stimulus s with the total number of stimulus-class being N s (N s = 10 in the MNIST experiment). Note that r i,s for the RNN depends on the time step as well. The value of SLI ranges from 0 (when the neuron responds identically to all stimuli) to 1 (when the neuron responds only to a single stimulus), and a higher SLI indicates a higher degree of selectivity, and vice versa. Interestingly, we find that most neurons have mixed selectivity shown in Fig. 8 (the right panel), whose SLI lies between 0 and 1, and such neurons respond strongly to a few types of images. The mixed selectivity was also discovered in complex cognitive tasks [24]. We also observe that some neurons are particularly selective to only one type of images with a high value of SLI [ Fig. 8 (middle)], while the remaining part of neurons keeps silent with an SLI near to zero [ Fig. 8 (left)]. As expected, the selectivity property emerges slightly before the decision making of the network. We finally remark that gBPTT produces a network ensemble, characterized by the hyper-parameter set of the SaS distribution. A concrete network can thus be sampled from this ensemble. During training, we use independently time-dependent Gaussian noise u (t) to approximate the statistics of the ensemble. Accordingly, we find that a timedependent concrete network yields the identical test performance with the mean-field propagation (i.e., using the parametrized pre-activation, see Eq. (4)]. The overall statistics of sampled networks does not vary significantly across dynamic steps (Fig. 9). This observation is in stark contrast to the traditional RNN training, where a deterministic weight matrix is used in all time steps. In fact, in a biological neural circuit, changes in synaptic connections are essential for the development and function of the nervous system [25,26]. In other words, the specific details of the connection pattern for a circuit may not be critical to the behavior, but rather, the parameters underlying the statistics of weight distributions become dominant factors for the behavioral output of the network (see also a recent paper demonstrating that innate face-selectivity could emerge from statistical variation of the feedforward projections in hierarchical neural networks [27]). Therefore, our current ensemble perspective of training RNNs offers a promising framework to achieve neural computation with a dynamic network, which is likely to be consistent with biologically plausible computations with fluctuating dendritic spines, e.g., dendritic spines can undergo morphological remodeling in adaptation to sensory stimuli or in learning [26].

B. Multisensory integration task
Multisensory integration is a fundamental ability of the brain to combine cues from multiple senses to form robust perception of signals in the noisy world [16,17]. In typical cognitive experiments, two separate sources of information  FIG. 10: Schematic illustration of the RNN performing a multi-sensory integration task. In one trial, the task consists of three consecutive stages (of lengths 300, 900, and 400 ms, respectively) including a fixation period (t0 ∼ 0), a stimulus period (0 ∼ t1), and a decision period (t1 ∼ t2). Note that in a test trial, the lengths of three stages are 500, 1000, and 300 ms, respectively. During the fixation period, the network must maintain a low output value, ensuring that the network will be sensitive to the stimulus signals. During the stimulus period, either or both (or multisensory) sources of signals with the same frequency f , including visual signal and auditory signal, act as the input to the network. The frequency f ranges from 9 to 16. Each type of sensory input contains positively tuned and negatively tuned signals (only the positively tuned one is shown in the plot). During the decision period, there are no modality input signal, and the network must learn to discriminate the even rate of the stimulus period, i.e., the output must hold a high or low value correctly. The output at the last time step of the decision period is assigned the final decision of one trial. Other parameters for training are as follows: x0 = 0.2 for all neurons, τ = 100 ms, α = 0.2, σ = 0.15, σin = 0.01, 2 regularization strength is 10 −3 , and the initial learning rate lr0 = 0.001. ∆t = 0.5 ms for testing.
are provided to animals, e.g., rats, before the animals make a decision. The source of information can be either auditory clicks or visual flashes, or both [28]. When the task becomes difficult, the multisensory training is more effective than the unisensory one. Akin to animals trained to perform a behavior task, RNNs can also learn the input-output mapping through an optimization procedure, which is able to provide a quantitative way to understand the computational principles underlying the multisensory integration. In particular, by increasing biological plausible levels of network architectures and even learning dynamics [8,10], one may generate quantitative hypotheses about the dynamical mechanisms that may be implemented in real neural circuits to solve the same behavior task. In this section, we restrict the network setting to the simplest case, i.e., a group of neurons are reciprocally connected with continuous firing rate dynamics. To consider more biological details is straightforward, e.g., taking Dale's principle [21]. The goal in this section is to show that our ensemble perspective also works in modeling cognitive computational tasks.
In the multisensory integration experiment (MSI), either or both of auditory and visual signals are fed into the recurrent network via deterministic weights (the same reason as in pixel-by-pixel MNIST task), which is required, after a stimulus period, to report whether the frequency of the presented stimulus is above a decision boundary (12.5 events per second)(see Fig. 10). One third of network units receive only visual input, while another third receive only auditory input, and the remaining third do not receive any input. The continuous variable r(t) in our model is an activity vector indicating the firing rates of neurons, obtained through a non-linear transfer function (ReLU here) of the synaptic current input [h(t)]. The current includes both of external input and recurrent feedback. The output z(t) is a weighted readout of the neural responses in the reservoir (a binary choice for the MSI task). The RNN is used here to model the multisensory event-rate discrimination task for the rats [21,28], and is trained by our gBPTT to solve the same audiovisual integration task. The visual and auditory inputs can be either positively (increasing function of event rate) or negatively tuned (decreasing function of event rate). Showing the network both types of tunned inputs could improve the training [29]. The RNN is composed of 150 neurons, whose recurrent dynamics is required to hold a high output value if the input event rate (represented by time-dependent noisy inputs) was above the decision boundary, and hold a low output value otherwise. Neurons are reciprocally connected, and the global statistics of the topology is learned from the training trials. The benefits of multisensory inputs, as known in cognitive science [16], are reproduced by our RNN model trained by gBPTT (Fig. 11). Integrating multiple sources of information, rather than unisensory inputs, helps decision making particularly when the task becomes hard (i.e., around the decision boundary). As the training proceeds, the sparsity of the recurrent layer grows rapidly until saturation, while the sparsity of the output layer grows in the same manner but finally reaching a lower value yet with a small fluctuation (Fig. 12). The recurrent layer becomes sparser with training, demonstrating that the latent dynamics is likely low dimensional, because of existence of some unnecessary degrees of freedom in recurrent feedbacks. In contrast, the goal of the output layer is to decode the recurrent dynamics, and the output layer should therefore keep all relevant dimensions of information integrated, which requires a densely connected output layer. This behavior is consistent with the evolution of the entropy profile. The recurrent layer maintains a relatively higher level of variability, compared with that of the output layer at the end of training (see the continuous entropy computed according to Eq. (14), or the Ξ-entropy in Fig. 13). In addition, the discrete π-entropy decreases with training in the output layer, in contrast to the increasing behavior of the same type of entropy in the recurrent layer (Fig. 13). Let us then look at the distribution profile of hyperparameters (Fig. 14). In the output layer, the distribution of π is U-shaped, and Ξ shows a sharp single peak at zero, which demonstrates that a dominant part of the weight distribution reduces to the Bernoulli distribution with two peaks at 0 and m = 0 respectively. The observed less variability in weight values of the output layer is consistent with the decoding stability. In the recurrent layer, the profile of π-distribution is U-shaped, and the distribution profile of Ξ resembles an L shape. This implies that there emerge VIP and UIP connections in the network. Moreover, a certain level of variability is allowed for the weight values, making a flexible computation possible in the internal dynamics.
Next, we ask whether our training leads to the emergence of neural selectivity, which indeed exists in the prefrontal cortex of behaving animals [24]. The selectivity properties of neurons play a critical role in the complex cognitive ability. In our trained networks, we also find that neurons in the recurrent layer display different types of selectivity (Fig. 15). In other words, neurons become highly active for either of modality (visual or auditory) and choice (high or low), or both (mixed selectivity).
We then explore the detailed patterns of the hyper-parameter matrices, which conveys individual contributions of each connection to the behavioral performance. By inspecting the sparsity matrix [ Fig. 16 (b)], one can identify both unnecessary (π = 1) and important connections (π = 0). We also find that, some neurons prefer receiving or sending information during the recurrent computation. An interpretation is that, the spatio-temporal information is divided into relevant and irrelevant parts; the relevant parts are maintained through sending preference, while the irrelevant parts are blocked through receiving preference.
Target weight perturbation experiments (Fig. 17) show that the VIP weights play a fundamental role in supporting the task accuracy reached by the recurrent computation. Our method can thus provide precise temporal credit assignment to the MSI task, which the standard BPTT could not.
Finally, we remark that a dynamic network with time-dependent specified weight values is able to reach an equivalent accuracy with the network using the mean-field propagation. Note that the overall statistics of the network does not change significantly (Fig. 18), suggesting that the hyperparameters for the the weight statistics are more important than precise values of weights. In fact, dendritic spines in neural circuits, biological substrates for synaptic contacts, are also subject to fluctuation, i.e., a highly dynamic structure [25,26]. Future exploration of this interesting connection would be fruitful, as an ensemble perspective is much more abstract than a concrete topology, while the specified stationary topology is still widely used in modern machine learning. Therefore, the ensemble perspective yielding a  dynamic network in adaptation of external stimuli could shed light on our understanding of adaptive RNNs.

IV. CONCLUSION
In this study, we propose an ensemble perspective for understanding temporal credit assignment, which yields a generalized BPTT algorithm for training RNNs in practice. Our training protocol produces the hyperparameters underlying the weight distribution, or the statistics of the RNN ensemble. In contrast to the standard BPTT, gBPTT highlights the importance of network statistics, which can particularly make a dynamic network changing its specified network weights a potential substrate for recurrent computation in adaptation to sensory inputs. It is thus interesting in future works to explore the biological counterpart of neural computation in brain circuits, e.g., in terms of dendritic spine dynamics [25].
Our SaS model has three types of hyperparameters with distinct roles. m tells us the mean of the continuous Gaussian slab, while Ξ controls the variance (fluctuation) of the slab. In both computational tasks we are interested in, we find that the Ξ-distribution profile is L-shaped. In other words, the peak at zero turns the SaS distribution into a Bernoulli distribution, while the tail at finite values of variance endows a computational flexibility to each connection, which may be critical to recode the high-dimensional sensory inputs into a low-dimensional latent dynamics. It is thus interesting to establish this hypothesis, by addressing precisely how the ensemble perspective helps to clarify the mechanism of low-dimensional latent dynamics encoding relevant spatio-temporal features of inputs.
The spike probability π tells us the sparsity level of the network. We find that a sparse RNN emerges after training. In particular, the recurrent layer is sparser than the output layer, reflecting different roles of both layers. The sparseness allows the recurrent layer to remove some unnecessary degrees of freedom in the recurrent dynamics, making recurrent feedbacks maintain only relevant information for decision making. In contrast, the output layer is averaged over ten independent runs. much denser, being thus able to read out relevant information both completely and robustly (note that most of Ξ values are zero). It is worth noticing that the π-distribution profile implies the existence of VIP weights, which we show has a critical contribution to the overall performance of the network. Our method thus provides a practical way to identify the critical elements of the network topology, which could contribute to understanding temporal credit assignment underlying the behavior of the RNN.
Working at the ensemble level, the constituent neurons of our model also display neural selectivity of different nature in response to the task parameters. Some neurons show uni-selectivity, while the others show mixed selectivity. The selectivity property of neurons, shaped by the recurrent connections, is also a key factor impacting the behavioral output of the network, and thus deserves future studies about the structure basis of the neural selectivity.
Taken together, our model and the associated training method are able to provide insights toward a precise understanding of temporal credit assignment, not only in engineering applications (e.g., MNIST digit classification), but also in modeling brain dynamics (e.g., multisensory integration task).