In our model, a learning task is defined by a mapping from task variables ${\displaystyle \mathbf{x}}$ to an output ${\displaystyle f(\mathbf{x})}$, representing a target change in activity of a readout neuron, for example a Purkinje cell. The limited scope of this definition implies our results should not strongly depend on the influence of the readout neuron on downstream circuits. The readout adjusts its incoming synaptic weights from the expansion layer to better approximate this target output. For example, for an associative learning task in which sensory stimuli are classified into categories such as appetitive or aversive, the task may be represented as a mapping from inputs to two discrete firing rates corresponding to the readout’s response to stimuli of each category (Figure 2A). In contrast, for a forward model, in which the consequences of motor commands are computed using a model of movement dynamics, an input encoding the current sensorimotor state is mapped to a continuous output representing the readout neuron’s tuning to a predicted sensory variable (Figure 2B).

Figure 2 with 4 supplements see all

Download asset Open asset

To examine how properties of the expansion layer representation influence learning performance across tasks, we designed two families of tasks: one modeling categorization of random stimuli, which is often used to study the performance of expanded neural representations (Rigotti et al., 2013; Barak et al., 2013; Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017; Cayco-Gajic et al., 2017), and the other modeling learning of a continuously varying output. The former we refer to as a ‘random categorization task’ and is parameterized by the number of input pattern-to-category associations $P$ learned during training (Figure 2C). During the training phase, the network learns to associate random input patterns ${\displaystyle {\mathbf{x}}^{\mu }\in {\mathbb{R}}^D}$ for $\mu =1,\mathrm{\dots },P$ with random binary categories $y^{\mu }=\pm 1$. The elements of ${\displaystyle {\mathbf{x}}^{\mu }}$ are drawn i.i.d. from a normal distribution with mean 0 and variance $1/D$. We refer to ${\displaystyle {\mathbf{x}}^{\mu }}$ as ‘training patterns’. To assess the network’s generalization performance, it is presented with ‘test patterns’ generated by adding noise (parameterized by a noise magnitude $ϵ$; see Methods) to the training patterns. Tasks with continuous outputs (Figure 2D) are parameterized by a length scale that determines how quickly the output changes as a function of the input (specifically, input-output functions are drawn from a Gaussian process with length scale $\gamma$ for variations in ${\displaystyle f(\mathbf{x})}$ as a function of ${\displaystyle \mathbf{x}}$; see Methods). In this case, both training and test patterns are drawn uniformly on the unit sphere. Later, we will also consider tasks implemented by specific cerebellum-like systems. See Table 1 for a summary of parameters throughout this study.

We trained the readout to approximate the target output for training patterns and generalize to unseen test patterns. The network’s prediction is ${\displaystyle \hat{f}(\mathbf{x})=\mathbf{w}\cdot \mathbf{h}(\mathbf{x})}$ for tasks with continuous outputs, or ${\displaystyle \hat{f}(\mathbf{x})=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\mathbf{w}\cdot \mathbf{h}(\mathbf{x}))}$ for categorization tasks, where ${\displaystyle \mathbf{w}}$ are the synaptic weights of the readout from the expansion layer. These weights were set using least squares regression. Performance was measured as the fraction of incorrect predictions for categorization tasks, or relative mean squared error for tasks with continuous targets: ${\displaystyle \mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\frac{\mathbb{E}[(f(\mathbf{x})-\hat{f}(\mathbf{x}){)}^2]}{\mathbb{E}[f(\mathbf{x}{)}^2]}}$, where the expectation is across test patterns.

We began by examining the dependence of learning performance on the coding level of the expansion layer. For random categorization tasks, performance is maximized at low coding levels (Figure 2E), consistent with previous results (Barak et al., 2013; Babadi and Sompolinsky, 2014). The optimal coding level remains below 0.1 in the model, regardless of the number of associations $P$, the level of input noise, and the dimension $D$ (Figure 2—figure supplement 1). For continuously varying outputs, the dependence is qualitatively different (Figure 2F). The optimal coding level depends strongly on the length scale, with learning performance for slowly varying functions optimized at much higher coding levels than quickly varying functions. This dependence holds for different choices of threshold-nonlinear functions (Figure 2—figure supplement 2) or input dimension (Figure 2—figure supplement 3) and is most pronounced when the number of training patterns is limited (Figure 2—figure supplement 4). Our observations are at odds with previous theories of the role of sparse granule cell representations (Marr, 1969; Albus, 1971; Babadi and Sompolinsky, 2014; Billings et al., 2014) and show that sparse activity does not always optimize performance for this broader set of tasks.

To determine how the optimal coding level depends on the task, we begin by quantifying how the expansion layer transforms the geometry of the task subspace. Later we will address how this transformation affects the ability of the network to learn a target. For ease of analysis, we will assume for now that inputs are normalized, ${\displaystyle \Vert \mathbf{x}\Vert =1}$, so that they lie on the surface of a sphere in $D$ dimensions. The set of neurons in the expansion layer activated by an input ${\displaystyle \mathbf{x}}$ are those neurons $i$ for which the alignment of their effective weights with the input, ${\displaystyle {\mathbf{J}}_i^{\mathrm{e}\mathrm{f}\mathrm{f}}\cdot \mathbf{x}}$, exceeds the activation threshold $\theta$ (Equation 1; Figure 3A). Increasing $\theta$ reduces the size of this set of neurons and hence reduces the coding level.

Figure 3

Download asset Open asset

Different inputs activate different sets of neurons, and more similar inputs activate sets with greater overlap. As the coding level is reduced, this overlap is also reduced (Figure 3B). In fact, this reduction in overlap is greater than the reduction in number of neurons that respond to either of the individual inputs, reflecting the fact that representations with low coding levels perform ‘pattern separation’ (Figure 3B, compare purple and red or blue regions).

This effect is summarized by the ‘kernel’ of the network (Schölkopf and Smola, 2002; Rahimi and Recht, 2007), which measures overlap of representations in the expansion layer as a function of the task variables:

Equations 1 and 2 show that the threshold $\theta$, which determines the coding level, influences the kernel through its effect on the expansion layer activity ${\displaystyle \mathbf{h}(\mathbf{x})}$. When inputs are normalized and the effective weights are Gaussian, we compute a semi-analytic expression for the kernel of the expansion layer in the limit of a large expansion ($M\to \mathrm{\infty }$; see Appendix). In this case, the kernel depends only on the overlap of the task variables, ${\displaystyle K(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})=K(\mathbf{x}\cdot {\mathbf{x}}^{\mathrm{\prime }})}$. Plotting the kernel for different choices of coding level demonstrates that representations with lower coding levels exhibit greater pattern separation (Figure 3C; Babadi and Sompolinsky, 2014). This is consistent with the observation that decreasing the coding level increases the dimension of the representation (Figure 3D).

We now relate the geometry of the expansion layer representation to performance across the tasks we have considered. Previous studies focused on high-dimensional, random categorization tasks in which inputs belong to a small number of well-separated clusters whose centers are random uncorrelated patterns. Generalization is assessed by adding noise to previously observed training patterns (Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017; Figure 4A). In this case, performance depends only on overlaps at two spatial scales: the overlap between training patterns belonging to different clusters, which is small, and the overlap between training and test patterns belonging to the same cluster, which is large (Figure 4B). For such tasks, the kernel evaluated near these two values—specifically, the behavior of ${\displaystyle K(\mathbf{x}\cdot {\mathbf{x}}^{\mathrm{\prime }})}$ near ${\displaystyle \mathbf{x}\cdot {\mathbf{x}}^{\mathrm{\prime }}=0}$ and ${\displaystyle \mathbf{x}\cdot {\mathbf{x}}^{\mathrm{\prime }}=1-\mathrm{\Delta }}$, where ${\displaystyle \mathrm{\Delta }}$ is a measure of within-cluster noise—fully determines generalization performance (Figure 4C; see Appendix). Sparse expansion layer representations reduce the overlap of patterns belonging to different clusters, increasing dimension and generalization performance (Figure 3D, Figure 2E).

Figure 4 with 2 supplements see all

Download asset Open asset

We study tasks where training patterns used for learning and test patterns used to assess generalization are both drawn according to a distribution over a low-dimensional space of task variables. While the mean overlap between pairs of random patterns remains zero regardless of dimension, fluctuations around the mean increase when the space is low dimensional, leading to a broader distribution of overlaps (Figure 4B). In this case, generalization performance depends on values of the kernel function evaluated across this entire range of overlaps. Methods from the theory of kernel regression (Sollich, 1998; Jacot et al., 2018; Bordelon et al., 2020; Canatar et al., 2021b; Simon et al., 2021) capture these effects by quantifying a network’s performance on a learning task through a decomposition of the target function into a set of basis functions (Figure 4D). Performance is assessed by summing the contribution of each mode in this decomposition to generalization error.

The decomposition expresses the kernel as a sum of eigenfunctions weighted by eigenvalues, ${\displaystyle K(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})=\sum _{\alpha }{\lambda }_{\alpha }{\psi }_{\alpha }(\mathbf{x}){\psi }_{\alpha }({\mathbf{x}}^{\mathrm{\prime }})}$. The eigenfunctions are determined by the network architecture and the distribution of inputs. As we show below, the eigenvalues ${\lambda }_{\alpha }$ determine the ease with which each corresponding eigenfunction ${\displaystyle {\psi }_{\alpha }(\mathbf{x})}$—one element of the basis function decomposition—is learned by the network. Under our present assumptions of Gaussian effective weights and uniformly distributed, normalized input patterns, the eigenfunctions are the spherical harmonic functions. These functions are ordered by increasing frequency, with higher frequencies corresponding to functions that vary more quickly as a function of the task variables. Spherical harmonics are defined for any input dimension; for example, in two dimensions they are the Fourier modes. We find that coding level substantially changes the frequency dependence of the eigenvalues associated with these eigenfunctions (Figure 4E). Higher coding levels increase the relative magnitude of the low frequency eigenvalues compared to high-frequency eigenvalues. As we will show, this results in a different inductive bias for networks with different coding levels.

To calculate learning performance for an arbitrary task, we decompose the target function in the same basis as that of the kernel:

The coefficient $c_{\alpha }$ quantifies the weight of mode $\alpha$ in the decomposition. For the Gaussian process targets, we have considered, increasing length scale corresponds to a greater relative contribution of low versus high frequency modes (Figure 4F). Using these coefficients and the eigenvalues (Figure 4E), we obtain an analytical prediction of the mean-squared generalization error (‘Error’) for learning any given task (Figure 4G; see Methods):

where C₁ and C₂ do not depend on $\alpha$ (Canatar et al., 2021b; Simon et al., 2021; see Methods). Equation 4 illustrates that for equal values of $c_{\alpha }$, modes with greater ${\lambda }_{\alpha }$ contribute less to the generalization error.

Our theory reveals that the optima observed in Figure 2F are a consequence of the difference in eigenvalues of networks with different coding levels. This reflects an inductive bias (Sollich, 1998; Jacot et al., 2018; Bordelon et al., 2020; Canatar et al., 2021b; Simon et al., 2021) of networks with low and high coding levels toward the learning of high and low frequency functions, respectively (Figure 4E, Figure 4—figure supplement 1). Thus, the coding level’s effect on a network’s inductive bias, rather than dimension alone, determines learning performance. Previous studies that focused only on random categorization tasks did not observe this dependence, since errors in such tasks are dominated by the learning of high frequency components, for which sparse activity is optimal (Figure 4—figure supplement 2).

To simplify our analysis, we have so far assumed full connectivity between input and expansion layers without a constraint on excitatory or inhibitory synaptic weights. In particular, we have assumed that the effective weight matrix ${\displaystyle {\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}$ contains independent Gaussian entries (Figure 5A, top). However, synaptic connections between mossy fibers and granule cells are sparse and excitatory (Sargent et al., 2005), with a typical in-degree of $K=4$ mossy fibers per granule cell (Figure 5A, bottom). We therefore analyzed the performance of model networks with more realistic connectivity. Surprisingly, when ${\displaystyle \mathbf{J}}$ is sparse and nonnegative, both overall generalization performance and the task-dependence of optimal coding level remain unchanged (Figure 5B).

Figure 5

Download asset Open asset

To understand this result, we examined how ${\displaystyle \mathbf{J}}$ and ${\displaystyle \mathbf{A}}$ shape the statistics of the effective weights onto the expansion layer neurons ${\displaystyle {\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}$. A desirable property of the expansion layer representation is that these effective weights sample the space of task variables uniformly (Figure 3A), increasing the heterogeneity of tuning of expansion layer neurons (Litwin-Kumar et al., 2017). This occurs when ${\displaystyle {\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}$ is a matrix of independent random Gaussian entries. If the columns of ${\displaystyle \mathbf{A}}$ are orthornormal and ${\displaystyle \mathbf{J}}$ is fully-connected with independent Gaussian entries, ${\displaystyle {\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}$ has this uniform sampling property.

However, when ${\displaystyle \mathbf{J}}$ is sparse and nonnegative, expansion layer neurons that share connections from the same input layer neurons receive correlated input currents. When $N$ is small and ${\displaystyle \mathbf{A}}$ is random, fluctuations in ${\displaystyle \mathbf{A}}$ lead to biases in the input layer’s sampling of task variables which are inherited by the expansion layer. We quantify this by computing the distribution of correlations between the effective weights for pairs of expansion layer neurons, ${\displaystyle \mathbf{C}\mathbf{o}\mathbf{r}\mathbf{r}\left({\mathbf{J}}_i^{\mathrm{e}\mathrm{f}\mathrm{f}},{\mathbf{J}}_j^{\mathrm{e}\mathrm{f}\mathrm{f}}\right)}$. This distribution indeed deviates from uniform sampling when $N$ is small (Figure 5C). However, even when $N$ is moderately large (but much less than $M$), only small deviations from uniform sampling of task variables occur for low dimensional tasks as long as ${\displaystyle D<K\ll N}$ (see Appendix). In contrast, for high-dimensional tasks ($D\sim N$), ${\displaystyle K\ll D}$ is sufficient, in agreement with previous findings (Litwin-Kumar et al., 2017). For realistic cerebellar parameters (${\displaystyle N=7,000}$ and $K=4$), the distribution is almost indistinguishable from that corresponding to uniform sampling (Figure 5C), consistent with the similar learning performance of these two cases (Figure 5B).

In the above analysis, an important assumption is that ${\displaystyle \mathbf{A}}$ is dense and random, so that the input layer forms a distributed representation in which each input layer neuron responds to a random combination of task variables (Figure 5D, top). If, on the other hand, the input layer forms a clustered representation containing groups of neurons that each encode a single task variable (Figure 5D, bottom), we may expect different results. Indeed, with a clustered representation, sparse connectivity dramatically reduces performance (Figure 5E). This is because the distribution of ${\displaystyle \mathbf{C}\mathbf{o}\mathbf{r}\mathbf{r}\left({\mathbf{J}}_i^{\mathrm{e}\mathrm{f}\mathrm{f}},{\mathbf{J}}_j^{\mathrm{e}\mathrm{f}\mathrm{f}}\right)}$ deviates substantially from that corresponding to uniform sampling (Figure 5F), even as $N\to \mathrm{\infty }$ (see Appendix). Specifically, increasing $N$ does not reduce the probability of two expansion layer neurons being connected to input layer neurons that encode the same task variables and therefore receiving highly correlated currents. As a result, expansion layer neurons do not sample task variables uniformly and performance is dramatically reduced.

Our results show that networks with small $K$, moderately large $N$, and a distributed input layer representation approach the performance of networks that sample task variables uniformly. This equivalence validates the applicability of our theory to these more realistic networks. It also argues for the importance of distributed sensorimotor representations in the cortico-cerebellar pathway, consistent with the distributed nature of representations in motor cortex (Shenoy et al., 2013; Muscinelli et al., 2023).

A history of theoretical modeling has shown a remarkable correspondence between anatomical properties of the cerebellar cortex and model parameters optimal for learning. These include the in-degree $K$ of granule cells (Marr, 1969; Litwin-Kumar et al., 2017; Cayco-Gajic et al., 2017), the expansion ratio of the granule cells to the mossy fibers $M/N$(Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017), and the distribution of synaptic weights from granule cells to Purkinje cells (Brunel et al., 2004; Clopath et al., 2012; Clopath and Brunel, 2013). In these studies, model performance was assessed using random categorization tasks. We have shown that optimal coding level is dependent on the task being learned, raising the question of whether optimal values of these architectural parameters are also task-dependent.

Sparse connectivity ($K$=4, consistent with the typical in-degree of cerebellar granule cells) has been shown to optimize learning performance in cerebellar cortex models (Litwin-Kumar et al., 2017; Cayco-Gajic et al., 2017). We examined the performance of networks with different granule cell in-degrees learning Gaussian process targets. The optimal in-degree is small for all the tasks we consider, suggesting that sparse connectivity is sufficient for high performance across a range of tasks (Figure 6A). This is consistent with the previous observation that the performance of a sparsely connected network approaches that of a fully connected network (Figure 5B).

Figure 6

Download asset Open asset

Previous studies also showed that the expansion ratio from mossy fibers to granule cells $M/N$ controls the dimension of the granule cell representation (Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017). The dimension increases with expansion ratio but saturates as expansion ratio approaches the anatomical value ($M/N\approx 30$ when ${\displaystyle f\approx 0.1}$ for the inputs presynaptic to an individual Purkinje cell). These studies assumed that mossy fiber activity is uncorrelated ($D=N$) rather than low-dimensional (${\displaystyle D<N}$). This raises the question of whether a large expansion is beneficial when $D$ is small. We find that when the number of training patterns $P$ is sufficiently large, performance still improves as $M/N$ approaches its anatomical value, showing that Purkinje cells can exploit their large number of presynaptic inputs even in the case of low-dimensional activity (Figure 6B).

Brunel et al., 2004 showed that the distribution of granule-cell-to-Purkinje cell synaptic weights is consistent with the distribution that maximizes the number of random binary input-output mappings stored. This distribution exhibits a substantial fraction of silent synapses, consistent with experiments. These results also hold for analog inputs and outputs (Clopath and Brunel, 2013) and for certain forms of correlations among binary inputs and outputs (Clopath et al., 2012). However, the case we consider, where targets are a smoothly varying function of task variables, has not been explored. We observe a similar weight distribution for these tasks (Figure 6C), with the fraction of silent synapses remaining high across coding levels (Figure 6D). The fraction of silent synapses is lower for networks learning Gaussian process targets than those learning random categorization tasks, consistent with the capacity of a given network for learning such targets being larger (Clopath et al., 2012).

Although optimal coding level is task-dependent, these analyses suggest that optimal architectural parameters are largely task-independent. Whereas coding level tunes the inductive bias of the network to favor the learning of specific tasks, these architectural parameters control properties of the representation that improve performance across tasks. In particular, sparse connectivity and a large expansion support uniform sampling of low-dimensional task variables (consistent with Figure 5C), while a large fraction of silent synapses is a consequence of a readout that maximizes learning performance (Brunel et al., 2004).

So far, we have considered analytically tractable families of tasks with parameterized input-output functions. Next, we extend our results to more realistic tasks constrained by the functions of specific cerebellum-like systems, which include both highly structured, continuous input-output mappings and random categorization tasks.

To model the cerebellum’s role in predicting the consequences of motor commands (Wolpert et al., 1998), we examined the optimal coding level for learning the dynamics of a two-joint arm (Fagg et al., 1997). Given an initial state, the network predicts the change in the future position of the arm (Figure 7A). Performance is optimized at substantially higher coding levels than for random categorization tasks, consistent with our previous results for continuous input-output mappings (Figure 2E and F).

Figure 7 with 2 supplements see all

Download asset Open asset

The mushroom body, a cerebellum-like structure in insects, is required for learning of associations between odors and appetitive or aversive valence (Modi et al., 2020). This behavior can be represented as a mapping from random representations of odors in the input layer to binary category labels (Figure 7B). The optimal coding level in a model with parameters consistent with the Drosophila mushroom body is less than 0.1, consistent with our previous results for random categorization tasks (Figure 2E) and the sparse odor-evoked responses in Drosophila Kenyon cells (Turner et al., 2008; Honegger et al., 2011; Lin et al., 2014).

The prediction and cancellation of self-generated sensory feedback has been studied extensively in mormyrid weakly electric fish and depends on the electrosensory lateral line lobe (ELL), a cerebellum-like structure (Bell et al., 2008). Granule cells in the ELL provide a temporal basis for generating negative images that are used to cancel self-generated feedback (Figure 7C). We extended a detailed model of granule cells and their inputs (Kennedy et al., 2014) to study the influence of coding level on the effectiveness of this basis. The performance of this model saturated at relatively high coding levels, and notably the coding level corresponding to biophysical parameters estimated from data coincided with the value at which further increases in performance were modest. This observation suggests that coding level is also optimized for task performance in this system.

A canonical function of the mammalian cerebellum is the adjustment of the vestibulo-ocular reflex (VOR), in which motion of the head is detected and triggers compensatory ocular motion in the opposite direction. During VOR learning, Purkinje cells are tuned to head velocity, and their tuning curves are described as piecewise linear functions (Lisberger et al., 1994; Figure 7D). Although in vivo population recordings of granule cells during VOR adaptation are not, to our knowledge, available for comparison, our model predicts that performance for learning such tuning curves is high across a range of coding levels and shows that sparse codes are sufficient (although not necessary) for such tasks (Figure 7D).

These results predict diverse coding levels across different behaviors dependent on cerebellum-like structures. The odor categorization and VOR tasks both have input-output mappings that exhibit sharp nonlinearities and can be efficiently learned using sparse representations. In contrast, the forward modeling and feedback cancellation tasks have smooth input-output mappings and exhibit denser optima. These observations are consistent with our previous finding that more structured tasks favor denser coding levels than do random categorization tasks (Figure 2E and F).

Previous studies assessed the learning performance of cerebellum-like systems with a model Purkinje cell that associates random patterns of mossy fiber activity with one of two randomly assigned categories (Marr, 1969; Albus, 1971; Brunel et al., 2004; Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017; Cayco-Gajic et al., 2017), a common benchmark for artificial learning systems (Gerace et al., 2022). In this case, a low coding level increases the dimension of the granule cell representation, permitting more associations to be stored and improving generalization to previously unseen inputs. The optimal coding level is low but not arbitrarily low, as extremely sparse representations introduce noise that hinders generalization (Barak et al., 2013; Babadi and Sompolinsky, 2014).

To examine a broader family of tasks, our learning problems extend previous studies in several ways. First, we consider inputs that may be constrained to a low-dimensional task subspace. Second, we consider input-output mappings beyond random categorization tasks. Finally, we assess generalization error for arbitrary locations on the task subspace, rather than only for noisy instances of previously presented inputs. As we have shown, these considerations require a complete analysis of the inductive bias of cerebellum-like networks (Figure 4). Our analysis generalizes previous approaches (Barak et al., 2013; Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017) that focused on dimension and noise alone. In particular, both dimension and noise for random patterns can be directly calculated from the kernel function (Figure 3C; see Appendix).

Our theory builds upon techniques that been developed for understanding properties of kernel regression (Sollich, 1998; Jacot et al., 2018; Bordelon et al., 2020; Canatar et al., 2021b; Simon et al., 2021). Kernel approximations of wide neural networks are a major area of current research providing analytically tractable theories (Rahimi and Recht, 2007; Jacot et al., 2018; Chizat et al., 2018). Prior studies have analyzed kernels corresponding to networks with zero (Cho and Saul, 2010) or mean-zero Gaussian thresholds (Basri et al., 2019; Jacot et al., 2018), which in both cases produce networks with a coding level of 0.5. Ours is the first kernel study of the effects of nonzero average thresholds. Our full characterization of the eigenvalue spectra and their decay rates as a function of the threshold extends previous work (Bach, 2017; Bietti and Bach, 2021). Furthermore, artificial neural network studies typically assume either fully-connected or convolutional layers, yet pruning connections after training barely degrades performance (Han et al., 2015; Zhang et al., 2018). Our results support the idea that sparsely connected networks may behave like dense ones if the representation is distributed (Figure 5), providing insight into the regimes in which pruning preserves performance.

Other studies have considered tasks with smooth input-output mappings and low-dimensional inputs, finding that heterogeneous Golgi cell inhibition can improve performance by diversifying individual granule cell thresholds (Spanne and Jörntell, 2013). Extending our model to include heterogeneous thresholds is an interesting direction for future work. Another proposal states that dense coding may improve generalization (Spanne and Jörntell, 2015). Our theory reveals that whether or not dense coding is beneficial depends on the task.

We have made several assumptions in our model for the sake of analytical tractability. When comparing the inductive biases of networks with different coding levels, our theory assumes that inputs are normalized and distributed uniformly in a linear subspace of the input layer activity. This allows us to decompose the target function into a basis in which we can directly compare eigenvalues, and hence learning performance, for different coding levels (Figure 4E–G). A similar analysis can be performed when inputs are not uniformly distributed, but in this case the basis is determined by an interplay between this distribution and the nonlinearity of expansion layer neurons, making the analysis more complex (see Appendix). We have also assumed that generalization is assessed for inputs drawn from the same distribution as used for learning. Recent and ongoing work on out-of-distribution generalization may permit relaxations of this assumption (Shen et al., 2021; Canatar et al., 2021a).

When analyzing properties of the granule cell layer, our theory also assumes an infinitely wide expansion. When $P$ is small enough that performance is limited by number of samples, this assumption is appropriate, but finite-size corrections to our theory are an interesting direction for future work. We also have not explicitly modeled inhibitory input provided by Golgi cells, instead assuming such input can be modeled as a change in effective threshold, as in previous studies (Billings et al., 2014; Cayco-Gajic et al., 2017; Litwin-Kumar et al., 2017). This is appropriate when considering the dimension of the granule cell representation (Litwin-Kumar et al., 2017), but more work is needed to extend our model to the case of heterogeneous inhibition.

Another key assumption concerning the granule cells is that they sample mossy fiber inputs randomly, as is typically assumed in Marr-Albus models (Marr, 1969; Albus, 1971; Litwin-Kumar et al., 2017; Cayco-Gajic et al., 2017). Other studies instead argue that granule cells sample from mossy fibers with highly similar receptive fields (Garwicz et al., 1998; Brown and Bower, 2001; Jörntell and Ekerot, 2006) defined by the tuning of mossy fiber and climbing fiber inputs to cerebellar microzones (Apps et al., 2018). This has led to an alternative hypothesis that granule cells serve to relay similarly tuned mossy fiber inputs and enhance their signal-to-noise ratio (Jörntell and Ekerot, 2006; Gilbert and Chris Miall, 2022) rather than to re-encode inputs. Another hypothesis is that granule cells enable Purkinje cells to learn piece-wise linear approximations of nonlinear functions (Spanne and Jörntell, 2013). However, several recent studies support the existence of heterogeneous connectivity and selectivity of granule cells to multiple distinct inputs at the local scale (Huang et al., 2013; Ishikawa et al., 2015). Furthermore, the deviation of the predicted dimension in models constrained by electron-microscopy data as compared to randomly wired models is modest (Nguyen et al., 2023). Thus, topographically organized connectivity at the macroscopic scale may coexist with disordered connectivity at the local scale, allowing granule cells presynaptic to an individual Purkinje cell to sample heterogeneous combinations of the subset of sensorimotor signals relevant to the tasks that Purkinje cell participates in. Finally, we note that the optimality of dense codes for learning slowly varying tasks in our theory suggests that observations of a lack of mixing (Jörntell and Ekerot, 2002) for such tasks are compatible with Marr-Albus models, as in this case nonlinear mixing is not required.

We have quantified coding level by the fraction of neurons that are above firing threshold. We focused on coding levels ${\displaystyle f<0.5}$, as extremely dense codes are rarely found in experiments (Olshausen and Field, 2004), but our theory applies for ${\displaystyle f>0.5}$ as well. In general, representations with coding levels of $f$ and $1-f$ perform similarly in our model due to the symmetry of most of their associated eigenvalues (Figure 4—figure supplement 1 and Appendix). Under the assumption that the energetic costs associated with neural activity are minimized, the ${\displaystyle f<0.5}$ region is likely the biologically plausible one. We also note that coding level is most easily defined when neurons are modeled as rate, rather than spiking units. To investigate the consistency of our results under a spiking code, we implemented a model in which granule cell spiking exhibits Poisson variability and quantify coding level as the fraction of neurons that have nonzero spike counts (Figure 7—figure supplement 1; Figure 7C). In general, increased spike count leads to improved performance as noise associated with spiking variability is reduced. Granule cells have been shown to exhibit reliable burst responses to mossy fiber stimulation (Chadderton et al., 2004), motivating models using deterministic responses or sub-Poisson spiking variability. However, further work is needed to quantitatively compare variability in model and experiment and to account for more complex biophysical properties of granule cells (Saarinen et al., 2008).

For the Purkinje cells, our model assumes that their responses to granule cell input can be modeled as an optimal linear readout. Our model therefore provides an upper bound to linear readout performance, a standard benchmark for the quality of a neural representation that does not require assumptions on the nature of climbing fiber-mediated plasticity, which is still debated. Electrophysiological studies have argued in favor of a linear approximation (Brunel et al., 2004). To improve the biological applicability of our model, we implemented an online climbing fiber-mediated learning rule and found that optimal coding levels are still task-dependent (Figure 7—figure supplement 2). We also note that although we model several timing-dependent tasks (Figure 7), our learning rule does not exploit temporal information, and we assume that temporal dynamics of granule cell responses are largely inherited from mossy fibers. Integrating temporal information into our model is an interesting direction for future investigation.

Our results predict that qualitative differences in the coding levels of cerebellum-like systems, across brain regions or across species, reflect an optimization to distinct tasks (Figure 7). However, it is also possible that differences in coding level arise from other physiological differences between systems. In the Drosophila mushroom body, which is required for associative learning of odor categories, random and sparse subsets of Kenyon cells are activated in response to odor stimulation, consistent with our model (Figure 7B; Turner et al., 2008; Honegger et al., 2011; Lin et al., 2014). In a model of the electrosensory system of the electric fish, the inferred coding level of a model constrained by the properties of granule cells is similar to that which optimizes task performance (Figure 7C). Within the cerebellar cortex, heterogeneity in granule cell firing has been observed across cerebellar lobules, associated with both differences in intrinsic properties (Heath et al., 2014) and mossy fiber input (Witter and De Zeeuw, 2015). It would be interesting to correlate such physiological heterogeneity with heterogeneity in function across the cerebellum. Our model predicts that regions involved in behaviors with substantial low-dimensional structure, for example smooth motor control tasks, may exhibit higher coding levels than regions involved in categorization or discrimination of high-dimensional stimuli.

Our model also raises the possibility that individual brain regions may exhibit different coding levels at different moments in time, depending on immediate behavioral or task demands. Multiple mechanisms could support the dynamic adjustment of coding level, including changes in mossy fiber input (Ozden et al., 2012), Golgi cell inhibition (Eccles et al., 1966; Palay and Chan-Palay, 1974), retrograde signaling from Purkinje cells (Kreitzer and Regehr, 2001), or unsupervised plasticity of mossy fiber-to-granule cell synapses (Schweighofer et al., 2001). The predictions of our model are not dependent on which of these mechanisms are active. A recent study demonstrated that local synaptic inhibition by Golgi cells controls the spiking threshold and hence the population coding level of cerebellar granule cells in mice (Fleming et al., 2022). Further, the authors observed that granule cell responses to sensory stimuli are sparse when movement-related selectivity is controlled for. This suggests that dense movement-related activity and sparse sensory-evoked activity are not incompatible.

While our analysis makes clear qualitative predictions concerning comparisons between the optimal coding levels for different tasks, in some cases it is also possible to make quantitative predictions about the location of the optimum for a single task. Doing so requires determining the appropriate time interval over which to measure coding level, which depends on the integration time constant of the readout neuron. It also requires estimates of the firing rates and biophysical properties of the expansion layer neurons. In the electrosensory system, for which a well-calibrated model exists and the learning objective is well-characterized (Kennedy et al., 2014), we found that the coding level estimated based on the data is similar to that which optimizes performance (Figure 7C).

If coding level is task-optimized, our model predicts that manipulating coding level artificially will diminish performance. In the Drosophila mushroom body, disrupting feedback inhibition from the GABAergic anterior paired lateral neuron onto Kenyon cells increases coding level and impairs odor discrimination (Lin et al., 2014). A recent study demonstrated that blocking inhibition from Golgi cells onto granule cells results in denser granule cell population activity and impairs performance on an eye-blink conditioning task (Fleming et al., 2022). These examples demonstrate that increasing coding level during sensory discrimination tasks, for which sparse activity is optimal, impairs performance. Our theory predicts that decreasing coding level during a task for which dense activity is optimal, such as smooth motor control, would also impair performance.

While dense activity has been taken as evidence against theories of combinatorial coding in cerebellar granule cells (Knogler et al., 2017; Wagner et al., 2019), our theory suggests that the two are not incompatible. Instead, the coding level of cerebellum-like regions may be determined by behavioral demands and the nature of the input to granule-like layers (Muscinelli et al., 2023). Sparse coding has also been cited as a key property of sensory representations in the cerebral cortex (Olshausen and Field, 1996). However, recent population recordings show that such regions exhibit dense movement-related activity (Musall et al., 2019), much like in cerebellum. While the theory presented in this study does not account for the highly structured recurrent interactions that characterize cerebrocortical regions, it is possible that these areas also operate using inductive biases that are shaped by coding level in a similar manner to our model.

The expansion layer activity is given by ${\displaystyle \mathbf{h}=\varphi ({\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}\mathbf{x}-\theta ),}$ where ${\displaystyle {\mathbf{J}}^{\mathrm{e}\mathrm{f}\mathrm{f}}=\mathbf{J}\mathbf{A}}$ describes the selectivity of expansion layer neurons to task variables. For most simulations, ${\displaystyle \mathbf{A}}$ is an $N\times D$ matrix sampled with random, orthonormal columns and ${\displaystyle \mathbf{J}}$ is an $M\times N$ matrix with i.i.d. unit Gaussian entries. The nonlinearity $\varphi$ is a rectified linear activation function $\varphi (u)=\max (u,0)$ applied element-wise. The input layer activity ${\displaystyle \mathbf{n}}$ is given by ${\displaystyle \mathbf{n}=\mathbf{A}\mathbf{x}.}$

We implemented a biophysical model of a planar two-joint arm (Fagg et al., 1997). The state of the arm is specified by six variables: joint angles ${\theta }_1$ and ${\theta }_2$, angular velocities ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$, and torques u₁ and u₂. The upper and lower segments of the arm have lengths l₁ and l₂ and masses m₁ and m₂, respectively. The arm has the following dynamics:

where ${\displaystyle \mathbf{M}\mathbf{(}\mathit{\theta }\mathbf{)}}$ is the inertia matrix and ${\displaystyle \mathbf{C}(\mathit{\theta },\dot{\mathit{\theta }})}$ is the matrix of centrifugal, Coriolis, and friction forces:

where ${\overline{l}}_2$ is the center of mass of the lower arm, I₁ and I₂ are moments of inertia and D₁ and D₂ are friction terms of the upper and lower arm respectively. These parameters were $m_1=3\text{kg}$, $m_2=2.5\text{kg}$, $l_1=0.3\text{m}$, $l_2=0.35\text{m}$, ${\overline{l}}_2=0.21\text{m}$, $I_1=0.1{\text{kg m}}^2$, $I_2=0.12{\text{kg m}}^2$, $D_1=0.05{\text{kg m}}^2/\text{s}$ and $D_2=0.01{\text{kg m}}^2/\text{s}$.

The task is to predict the position of the hand based on the forward dynamics of the two-joint arm system, given the arm initial condition and the applied torques. More precisely, the $P$ network inputs ${\displaystyle {\mathbf{x}}^{\mu }}$ were generated by sampling 6-dimensional Gaussian vectors with covariance matrix ${\displaystyle \mathbf{C}=\text{diag}({\sigma }_{\theta },{\sigma }_{\theta },{\sigma }_{\dot{\theta }},{\sigma }_{\dot{\theta }},{\sigma }_u,{\sigma }_u)}$, to account for the fact that angles, angular velocities and torques might vary on different scales across simulations. For our results, we used ${\sigma }_{\theta }={\sigma }_{\dot{\theta }}=0.1$ and ${\sigma }_u=1$. Each sample ${\displaystyle {\mathbf{x}}^{\mu }}$ was then normalized and used to generate initial conditions of the arm, by setting ${\theta }_1^{\mu }=\frac{\pi }{4}+x_1^{\mu }$, ${\theta }_2^{\mu }=\frac{\pi }{4}+x_2^{\mu }$, ${\dot{\theta }}_1^{\mu }=x_3^{\mu }$, and ${\dot{\theta }}_2^{\mu }=x_4^{\mu }$. Torques were generated by setting $u_1^{\mu }=x_5^{\mu }$ and $u_2^{\mu }=x_6^{\mu }$. The target was constructed by running the dynamics of the arm forward in time for a time $\delta =0.2\text{s}$, and by computing the difference in position of the “hand” (i.e. the end of the lower segment) in Cartesian coordinates. As a result, the target in this task is two-dimensional, with each target dimension corresponding the one of the two Cartesian coordinates of the hand. The overall performance is assessed by computing the error on each task separately and then averaging the errors.

We simulated 20,000 granule cells using the biophysical model of Kennedy et al., 2014. We varied the granule cell layer coding level by adjusting the spiking threshold parameter in the model. For each choice of threshold, we generated 30 different trials of spike rasters. Each trial is 160ms long with a 1ms time bin and consists of a time-locked response to an electric organ discharge command. Trial-to-trial variability in the model granule cell responses arises from noise in the mossy fiber responses. To generate training and testing data, we sampled 4 trials ($P=640$ patterns) from the 30 total trials for training and 10 trials for testing (1600 patterns). Coding level is measured as the fraction of granule cells that spike at least once in the training data. We repeated this sampling process 30 times.

The targets were smoothed broad-spike responses of 15 MG cells time-locked to an electric organ discharge command measured during experiments (Muller et al., 2019). The original data set consisted of 55 MG cells, each with a 300ms long spike raster with a 1ms time bin. The spike rasters were trial-averaged and then smoothed with a Gaussian-weighted moving average with a 10ms time window. Only MG cells whose maximum spiking probability across all time bins exceeded 0.01 after smoothing were included in the task. The same MG cell responses were used for both training and testing. To match the length of the granule cell data, we discarded MG cell data beyond 160ms and then concatenated 4 copies of the 160ms long responses for training and 10 copies for testing. We measured the ability of the model to construct MG cell targets out of granule cell activity, generalizing across noise in granule cell responses. Errors for each MG cell target were averaged across the 30 repetitions of sampling of training and testing data, and then averaged across targets. Standard error of the mean was computed across the 30 repetitions.

Recordings of Purkinje cell activity in monkeys suggest that these neurons exhibit piecewise-linear tuning to head velocity (Lisberger et al., 1994). Thus, we designed piecewise-linear target functions representing Purkinje cell firing rate as a function of head velocity $v$, a one-dimensional input:

Inputs $v$ were sampled uniformly from $[-1,1]$ 100 times. We generated 25 total target functions using all combinations of slopes m₁ and m₂ sampled from 5 equally spaced points on the interval $[-2,2]$. We set $b=0.1$ and $c=-0.2$.

Mossy fiber responses to head velocity input were modeled as exponential tuning curves:

where g_j is a gain term, $r_j\in \pm 1$ determines a mossy fiber preference for positive or negative velocities, and b_j is the baseline firing rate. We generated 24 different tuning curves from all combinations of the following parameter values: The gain g_j was sampled from 6 equally spaced points on the interval $[0.1,1]$, r_j was set to either –1 or 1, and b_j was set to either 0 or 1. Qualitative results did not depend strongly on this parameterization. Mossy fiber to granule cell weights were random zero-mean Gaussians. Errors were averaged across targets.

Previous theories (Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017) studied generalization performance for random clusters of inputs associated with binary targets, where test patterns are formed by adding noise to training patterns (Figure 4A). The readout is trained using a supervised Hebbian rule with mean-subtracted expansion layer responses, ${\displaystyle \mathbf{w}=\sum _{\mu }{y}^{\mu }({\mathbf{h}}^{\mu }-\overline{\mathbf{h}})}$, with ${\displaystyle \overline{\mathbf{h}}=\frac{1}{P}\sum _{\mu =1}^P{\mathbf{h}}^{\mu }}$. The net input to a readout in response to a test pattern ${\displaystyle {\hat{\mathbf{h}}}^{\mu }}$ from cluster µ is ${\displaystyle g^{\mu }=\mathbf{w}\cdot ({\hat{\mathbf{h}}}^{\mu }-\overline{\mathbf{h}})}$. The statistics of $g^{\mu }$ determine generalization performance. For a Hebbian readout, the error rate is expressed in terms of the signal-to-noise ratio (SNR) (Babadi and Sompolinsky, 2014):

SNR is given in terms of the mean and variance of $g^{\mu }$:

The numerator of SNR is proportional to the average overlap of the expansion layer representations of training and test patterns belonging to the same cluster, which can be expressed in terms of the kernel function $K$:

For large networks with Gaussian i.i.d. expansion weights, ${\displaystyle K({\hat{\mathbf{x}}}^{\mu },{\mathbf{x}}^{\mu })=K(t)}$, where ${\displaystyle t={\hat{\mathbf{x}}}^{\mu }\cdot {\mathbf{x}}^{\mu }}$, and the above equation reduces to ${\displaystyle MK(t_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}/\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}})-\overline{\mathbf{h}}\cdot \overline{\mathbf{h}}}$, where $t_{\mathrm{train}/\mathrm{test}}$ is the typical overlap of training and test patterns belonging to the same cluster. When ${\displaystyle {\mathbf{x}}^{\mu }\cdot {\mathbf{x}}^{\mu }=1}$, $t_{\mathrm{train}/\mathrm{test}}$ can be written as ${\displaystyle t_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}/\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}=1-\mathrm{\Delta }}$, where Δ is a measure of within-cluster noise (Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017).

Babadi and Sompolinsky, 2014 demonstrated that, for random categorization tasks and when $M$ and $D$ are large, $\mathrm{Var}\left(g^{\mu }\right)=C(\frac{1}{M}+Q^2\frac{1}{D})$ where $C$ is a constant and $Q\in [0,1]$ is given by

assuming the entries of ${\displaystyle \mathbf{x}}$ are zero-mean. ${\displaystyle Z_a=\underset{\mu }{\mathbb{E}}[\Vert {\mathbf{a}}^{\mu }-\overline{\mathbf{a}}{\Vert }^2]}$ for $a\in \{h,x\}$ normalizes the overlaps to the typical overlap of a pattern with itself. The quantity $Q^2$ is the ratio of the variance of overlaps between patterns belonging to different clusters in the expansion layer to that of the input layer. This describes the extent to which the geometry of the input layer representation is preserved in the expansion layer. When overlaps in the input layer are small, as they are for random clusters, ${\displaystyle \frac{1}{Z_h}({\mathbf{h}}^{\mu }-\overline{\mathbf{h}})\cdot ({\mathbf{h}}^{\nu }-\overline{\mathbf{h}})\approx \frac{Q}{Z_x}\cdot ({\mathbf{x}}^{\mu }\cdot {\mathbf{x}}^{\nu })}$ as $M\to \mathrm{\infty }$. This relation illustrates that, for random clusters and $M\to \mathrm{\infty }$, $Q$ is equal to the slope of the normalized kernel function $K(t)$ evaluated at $t=0$. Litwin-Kumar et al., 2017 also showed that the dimension of the expansion layer representation is equal to $\frac{C^{\prime }}{(\frac{1}{M}+Q^2\frac{1}{D})}$, where $C^{\prime }$ is a constant.

Thus, for the random categorization task studied in Babadi and Sompolinsky, 2014; Litwin-Kumar et al., 2017, dimension and readout SNR can be calculated by evaluating $K(t_{\mathrm{train}/\mathrm{test}})$ and the slope of $K(t)$ at $t=0$.

As $M\to \mathrm{\infty }$, the normalized dot product between features (Equation 2) converges pointwise to

where ${\displaystyle {\mathbf{J}}_1}$ is a row of the weight matrix ${\displaystyle \mathbf{J}}$ (without loss of generality, the first row) with entries drawn i.i.d. from a Gaussian distribution ${\displaystyle \mathcal{N}(0,1)}$. Our goal is to compute Equation A5 for a given $\theta$ and inputs drawn on the unit sphere ${\displaystyle \mathbf{x},{\mathbf{x}}^{\mathrm{\prime }}\in {\mathbb{S}}^{D-1}}$.

Because the Gaussian weight distribution is spherically symmetric, Equation A5 restricted to the unit sphere for any nonlinearity is only a function of the dot-product ${\displaystyle t:={\mathbf{x}}^T{\mathbf{x}}^{\mathrm{\prime }}}$, making the kernel a dot-product kernel ${\displaystyle K(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})=K(t)}$.

Denote by $J_i$ the entries of ${\displaystyle {\mathbf{J}}_1}$. Let $I_1={\sum }_{i=1}^D{J}_i{x}_i$ and $I_2={\sum }_{i=1}^D{J}_i{x}_i^{\prime }$ be the pre-activations for each input. Then $(I_1,I_2)$ are jointly Gaussian with mean 0, variance 1, and covariance ${\displaystyle \mathbb{E}[I_1{I}_2]=t}$. If ${\displaystyle t>0}$, we can re-parameterize these pre-activations as the sum of an independent and shared component ${\displaystyle I_i=y_i\sqrt{1-t}+z\sqrt{t}}$, where ${\displaystyle y_i\sim \mathcal{N}(0,1)}$ for ${\displaystyle i=1,2}$ and ${\displaystyle z\sim \mathcal{N}(0,1)}$. In these coordinates, Equation A5 becomes

where the second line follows from the conditional independence of $h_1|z$ and $h_2|z$ and the third from the fact that they are identically distributed. Similarly, if ${\displaystyle t<0}$, we can write $I_1=y_1\sqrt{1-t}+z\sqrt{t}$, $I_2=y_2\sqrt{1-t}-z\sqrt{t}$.

We will use Equation A6 to solve for the kernel assuming $\varphi$ is a ReLU nonlinearity. Let

Using the fact that ${\displaystyle \varphi }$ is nonzero only when ${\displaystyle y_1\sqrt{1-t}+z\sqrt{t}-\theta >0}$, i.e. for ${\displaystyle y_1>T=\frac{\theta -z\sqrt{t}}{\sqrt{1-t}}}$, we obtain

Performing a similar calculation for ${\displaystyle t<0}$ and collecting the results leads to:

Our theory of generalization requires us to work in function spaces which are natural to the problem. The spherical harmonics are the natural basis for working with dot-product kernels on the sphere. For a thorough treatment of spherical harmonics, see Atkinson and Han, 2012, whose notation we generally follow. Both our kernel and Gaussian process (GP) tasks are defined over the sphere in $D$ dimensions

A spherical harmonic $Y_{km}(\cdot )$—where $k$ indexes frequency and $m$ indexes modes of the same frequency—is a harmonic homogeneous polynomial of degree $k$ restricted to the sphere ${\displaystyle {\mathbb{S}}^{D-1}}$. For each frequency $k\in \mathbb{Z}$, there are $N(D,k)$ linearly independent polynomials, where