Population codes enable learning from few examples by shaping inductive bias
Abstract
Learning from a limited number of experiences requires suitable inductive biases. To identify how inductive biases are implemented in and shaped by neural codes, we analyze sampleefficient learning of arbitrary stimulusresponse maps from arbitrary neural codes with biologicallyplausible readouts. We develop an analytical theory that predicts the generalization error of the readout as a function of the number of observed examples. Our theory illustrates in a mathematically precise way how the structure of population codes shapes inductive bias, and how a match between the code and the task is crucial for sampleefficient learning. It elucidates a bias to explain observed data with simple stimulusresponse maps. Using recordings from the mouse primary visual cortex, we demonstrate the existence of an efficiency bias towards lowfrequency orientation discrimination tasks for grating stimuli and low spatial frequency reconstruction tasks for natural images. We reproduce the discrimination bias in a simple model of primary visual cortex, and further show how invariances in the code to certain stimulus variations alter learning performance. We extend our methods to timedependent neural codes and predict the sample efficiency of readouts from recurrent networks. We observe that many different codes can support the same inductive bias. By analyzing recordings from the mouse primary visual cortex, we demonstrate that biological codes have lower total activity than other codes with identical bias. Finally, we discuss implications of our theory in the context of recent developments in neuroscience and artificial intelligence. Overall, our study provides a concrete method for elucidating inductive biases of the brain and promotes sampleefficient learning as a general normative coding principle.
Editor's evaluation
This important study presents a theory of generalization in neural population codes and proposes sample efficiency as a new normative principle. The theory can be used to identify the set of 'easily learnable' stimulusresponse mappings from neural data and makes strong behavioral predictions that can be evaluated experimentally. Overall, the new method for elucidating inductive biases of the brain is highly compelling and will be of interest to theoretical and experimental neuroscientists working towards understanding how the cortex works.
https://doi.org/10.7554/eLife.78606.sa0Introduction
The ability to learn quickly is crucial for survival in a complex and an everchanging environment, and the brain effectively supports this capability. Often, only a few experiences are sufficient to learn a task, whether acquiring a new word (Carey and Bartlett, 1978) or recognizing a new face (Peterson et al., 2009). Despite the importance and ubiquity of sample efficient learning, our understanding of the brain’s information encoding strategies that support this faculty remains poor (Tenenbaum et al., 2011; Lake et al., 2017; Sinz et al., 2019).
In particular, when learning and generalizing from past experiences, and especially from few experiences, the brain relies on implicit assumptions it carries about the world, or its inductive biases (Wolpert, 1996; Sinz et al., 2019). Reliance on inductive bias is not a choice: inferring a general rule from finite observations is an illposed problem which requires prior assumptions since many hypotheses can explain the same observed experiences (Hume, 1998). Consider learning a rule that maps photoreceptor responses to a prediction of whether an observed object is a threat or is neutral. Given a limited number of visual experiences of objects and their threat status, many threatdetection rules are consistent with these experiences. By choosing one of these threatdetection rules, the nervous system reveals an inductive bias. Without the right biases that suit the task at hand, successful generalization is impossible (Wolpert, 1996; Sinz et al., 2019). In order to understand why we can quickly learn to perform certain tasks accurately but not others, we must understand the brain’s inductive biases (Tenenbaum et al., 2011; Lake et al., 2017; Sinz et al., 2019).
In this paper, we study sample efficient learning and inductive biases in a general neural circuit model which comprises of a population of sensory neurons and a readout neuron learning a stimulusresponse map with a biologicallyplausible learning rule (Figure 1A). For this circuit and learning rule, inductive bias arises from the nature of the neural code for sensory stimuli, specifically its similarity structure. While different population codes can encode the same stimulus variables and allow learning of the same output with perfect performance given infinitely many samples, learning performance can depend dramatically on the code when restricted to a small number of samples, where the reliance on and the effect of inductive bias are strong (Figure 1B, C and D). Given the same sensory examples and their associated response values, the readout neuron may make drastically different predictions depending on the inductive bias set by the nature of the code, leading to successful or failing generalizations (Figure 1C and D). We say that a code and a learning rule, together, have a good inductive bias for a task if the task can be learned from a small number of examples.
In order to understand how population codes shape inductive bias and allow fast learning of certain tasks over others with a biologically plausible learning rule, we develop an analytical theory of the readout neuron’s learning performance as a function of the number of sampled examples, or sample size. We find that the readout’s performance is completely determined by the code’s kernel, a function which takes in pairs of population response vectors and outputs a representational similarity defined by the inner product of these vectors. We demonstrate that the spectral properties of the kernel introduce an inductive bias toward explaining sampled data with simple stimulusresponse maps and determine compatibility of the population code with the learning task, and hence the sampleefficiency of learning. We apply this theory to data from the mouse primary visual cortex (V1) (Stringer et al., 2021; Pachitariu et al., 2019; Stringer et al., 2018a; Stringer et al., 2018b), and show that mouse V1 responses support sampleefficient learning of low frequency orientation discrimination and low spatial frequency reconstruction tasks over high frequency ones. We demonstrate the discrimination bias in a simple model of V1 and show how response nonlinearity, sparsity, and relative proportion of simple and complex cells influence the code’s bias and performance on learning tasks, including ones that involve invariances. We extend our theory to temporal population codes, including codes generated by recurrent neural networks learning a delayed response task. We observe that many codes could support the same kernel function, however, by analyzing data from mouse primary visual cortex (V1) (Stringer et al., 2021; Pachitariu et al., 2019; Stringer et al., 2018a; Stringer et al., 2018b), we find that the biological code is metabolically more efficient than others.
Overall, our results demonstrate that for a fixed learning rule, the neural sensory representation imposes an inductive bias over the space of learning tasks, allowing some tasks to be learned by a downstream neuron more sampleefficiently than others. Our work provides a concrete method for elucidating inductive biases of populations of neurons and suggest sampleefficient learning as a novel functional role for population codes.
Results
Problem setup
We denote vectors with bold lowercase symbols $\mathbf{r}$ and matrices $\mathbf{K}$ with bold uppercase symbols. We denote an average of a function $g(\mathit{\theta})$ over random variable $\mathit{\theta}$ as $\u27e8g(\mathit{\theta})\u27e9}_{\mathit{\theta}$. Euclidean inner products between vectors are denoted either as $\mathbf{x}\cdot \mathbf{y}$ or ${\mathbf{x}}^{\top}\mathbf{y}$ and real Euclidean $n$space is denoted ${\mathbb{R}}^{n}$. Sets of variables are represented with $\{\cdot \}$.
We consider a population of $N$ neurons whose responses, $\{{r}_{1}(\mathit{\theta}),{r}_{2}(\mathit{\theta}),\mathrm{\dots},{r}_{N}(\mathit{\theta})\}$, vary with the input stimuli, which is parameterized by a vector variable $\mathit{\theta}\in {\mathbb{R}}^{d}$, such as the orientation and the phase of a grating (Figure 1A). These responses define the population code. Throughout this work, we will mostly assume that this population code is deterministic: that identical stimuli generate identical neural responses.
From the population responses, a readout neuron learns its weights $\mathbf{w}$ to approximate a stimulusresponse map, or a target function $y(\mathit{\theta})$, such as one that classifies stimuli as apetitive ($y=1$) or aversive ($y=1$), or a more smooth one that attaches intermediate values of valence. We emphasize that in our model only the readout neuron performs learning, and the population code is assumed to be static through learning. Our theory is general in its assumptions about the structure of the population code and the stimulusresponse map considered (Methods Theory of generalization), and can apply to many scenarios.
The readout neuron learns from $P$ stimulusresponse examples with the goal of generalizing to previously unseen ones. Example stimuli ${\mathit{\theta}}^{\mu}$, ($\mu =1,\mathrm{\dots},P$) are sampled from a probability distribution describing stimulus statistics $p(\mathit{\theta})$. This distribution can be natural or artificially created, for example, for a laboratory experiment (Appendix Discrete stimulus spaces: finding eigenfunctions with matrix eigendecomposition). From the set of learning examples, $\mathcal{D}=\{{\mathit{\theta}}^{\mu},y({\mathit{\theta}}^{\mu}){\}}_{\mu =1}^{P}$, the readout weights are learned with the local, biologicallyplausible deltarule, $\mathrm{\Delta}{w}_{j}=\eta {\sum}_{\mu}{r}_{j}({\mathit{\theta}}^{\mu})(y({\mathit{\theta}}^{\mu})\mathbf{r}({\mathit{\theta}}^{\mu})\cdot \mathbf{w})$,where $\eta $ is a learning rate (Figure 1A). Learning with weight decay, which privileges readouts with smaller norm, can also be accommodated in our theory as we discuss in (Appendix Weight decay and ridge regression). With or without weight decay, the learning rule converges to a unique set of weights ${\mathbf{w}}^{\ast}(\mathcal{D})$ (Appendix Convergence of the deltarule without weight decay). Generalization error with these weights is given by
which quantifies the expected error of the trained readout over the entire stimulus distribution $p(\mathit{\theta})$. This quantity will depend on the population code $\mathbf{r}(\mathit{\theta})$, the target function $y(\mathit{\theta})$ and the set of training examples $\mathcal{D}$. Our theoretical analysis of this model provides insights into how populations of neurons encode information and allow sampleefficient learning.
Kernel structure of population codes controls learning performance
First, we note that the generalization performance of the learned readout on a given task depends entirely on the inner product kernel, defined by $K(\mathit{\theta},{\mathit{\theta}}^{\prime})=\frac{1}{N}{\sum}_{i=1}^{N}{r}_{i}(\mathit{\theta}){r}_{i}({\mathit{\theta}}^{\prime})$, which quantifies the similarity of population responses to two different stimuli $\mathit{\theta}$ and $\mathit{\theta}}^{\mathrm{\prime}$. The kernel, or similarity matrix, encodes the geometry of the neural responses. Concretely, distances (in neural space) between population vectors for stimuli $\mathit{\theta},{\mathit{\theta}}^{\prime}$ can be computed from the kernel $\frac{1}{N}{\mathbf{r}(\mathit{\theta})\mathbf{r}({\mathit{\theta}}^{\prime})}^{2}=K(\mathit{\theta},\mathit{\theta})+K({\mathit{\theta}}^{\prime},{\mathit{\theta}}^{\prime})2K(\mathit{\theta},{\mathit{\theta}}^{\prime})$ (Edelman, 1998; Kriegeskorte et al., 2008; Laakso and Cottrell, 2000; Kornblith et al., 2019; Cadieu et al., 2014; Pehlevan et al., 2018). The fact that the solution to the learning problem only depends on the kernel is due to the convergence of the learning rule to a unique solution ${\mathbf{w}}^{\ast}(\mathcal{D})$ for the training set $\mathcal{D}$ (Neal, 1994; Girosi et al., 1995). The datasetdependent fixed point ${\mathbf{w}}^{\ast}(\mathcal{D})$ of the learning rule is a linear combination of the population vectors on the dataset ${\mathbf{w}}^{\ast}(\mathcal{D})=\frac{1}{N}\sum _{\mu =1}^{P}{\alpha}^{\mu}\mathbf{r}({\mathit{\theta}}^{\mu})$. Thus, the learned function computed by the readout neuron is
where the coefficient vector satisfies $\mathit{\alpha}={\mathbf{K}}^{+}\mathbf{y}$ (Appendix Convergence of the deltarule without weight decay), and the matrix $\mathbf{K}$ has entries ${K}_{\mu \nu}=K({\mathit{\theta}}^{\mu},{\mathit{\theta}}^{\nu})$ and ${y}_{\mu}=y({\mathit{\theta}}^{\mu})$. The matrix K^{+} is the pseudoinverse of $\mathbf{K}$. In these expressions the population code only appears through the kernel $K$, showing that the kernel alone controls the learned response pattern. This result applies also to nonlinear readouts (Appendix Convergence of Deltarule for nonlinear readouts), showing that the kernel can control the learned solution in a variety of cases.
Since predictions only depend on the kernel, a large set of codes achieve identical desired performance on learning tasks. This is because the kernel is invariant with respect to rotation of the population code. An orthogonal transformation $\mathbf{Q}$ applied to a population code $\mathbf{r}(\mathit{\theta})$ generates a new code $\stackrel{~}{\mathbf{r}}(\mathit{\theta})=\mathbf{Qr}(\mathit{\theta})$ with an identical kernel (Appendix Alternative neural codes with same kernel) since $\frac{1}{N}\stackrel{~}{\mathbf{r}}(\mathit{\theta})\cdot \stackrel{~}{\mathbf{r}}({\mathit{\theta}}^{\prime})=\frac{1}{N}\mathbf{r}{(\mathit{\theta})}^{\top}{\mathbf{Q}}^{\top}\mathbf{Qr}({\mathit{\theta}}^{\prime})=\frac{1}{N}\mathbf{r}(\mathit{\theta})\cdot \mathbf{r}({\mathit{\theta}}^{\prime})$. Codes $\mathbf{r}(\mathit{\theta})$ and $\stackrel{~}{\mathbf{r}}(\mathit{\theta})$ will have identical readout performance on all possible learning tasks. We illustrate this degeneracy in Figure 2 using a publicly available dataset which consists of activity recorded from ∼20,000 neurons from the primary visual cortex of a mouse while shown static gratings (Stringer et al., 2021; Pachitariu et al., 2019). An original code $\mathbf{r}(\mathit{\theta})$ is rotated to generate $\stackrel{~}{\mathbf{r}}(\mathit{\theta})$ (Figure 2A) which have the same kernels (Figure 2B) and the same performance on a learning task (Figure 2C).
Codetask alignment governs generalization
We next examine how the population code affects generalization performance of the readout. We calculated analytical expressions of the average generalization error in a task defined by the target response $y(\mathit{\theta})$ after observing $P$ stimuli using methods from statistical physics (Methods Theory of generalization). Because the relevant quantity in learning performance is the kernel, we leveraged results from our previous work studying generalization in kernel regression (Bordelon et al., 2020; Canatar et al., 2021), and approximated the generalization error averaged over all possible realizations of the training dataset composed of $P$ stimuli, $E}_{g}={\u27e8{E}_{g}(\mathcal{D})\u27e9}_{\mathcal{D}$. As $P$ increases, the variance in ${E}_{g}$ due to the composition of the dataset decreases, and our expressions become descriptive of the typical case. Our final analytical result is given in Equation (11) in Methods Theory of generalization. We provide details of our calculations in Methods Theory of generalization and Appendix Theory of generalization, and focus on their implications here.
One of our main observations is that given a population code $\mathbf{r}(\mathit{\theta})$, the singular value decomposition of the code gives the appropriate basis to analyze the inductive biases of the readouts (Figure 3A). The tuning curves for individual neurons ${r}_{i}(\mathit{\theta})$ form an $N$by$M$ matrix $\mathbf{R}$, where $M$, possibly infinite, is the number of all possible stimuli. We discuss the SVD for continuous stimulus spaces in Appendix Singular value decomposition of continuous population responses. The leftsingular vectors (or principal axes) and singular values of this matrix have been used in neuroscience for describing lower dimensional structure in the neural activity and estimating its dimensionality, see e.g. (Stopfer et al., 2003; Kato et al., 2015; Bathellier et al., 2008; Gallego et al., 2017; Sadtler et al., 2014; Stringer et al., 2018b, Stringer et al., 2021; LitwinKumar et al., 2017; Gao et al., 2017; Gao and Ganguli, 2015). We found that the function approximation properties of the code are controlled by the singular values, or rather their squares $\{{\lambda}_{k}\}$ which give variances along principal axes, indexed in decreasing order, and the corresponding right singular vectors $\{{\psi}_{k}(\mathit{\theta})\}$, which are also the kernel eigenfunctions (Methods Theory of generalization and Appendix Singular value decomposition of continuous population responses). This follows from the fact that learned response (Equation (2)) is only a function of the kernel $K$, and the eigenvalues ${\lambda}_{k}$ and orthonormal (uncorrelated) eigenfunctions ${\psi}_{k}(\mathit{\theta})$ collectively define the code’s innerproduct kernel $K(\mathit{\theta},{\mathit{\theta}}^{\prime})$ through an eigendecomposition $K(\mathit{\theta},{\mathit{\theta}}^{\prime})=\frac{1}{N}{\sum}_{i=1}^{N}{r}_{i}(\mathit{\theta}){r}_{i}({\mathit{\theta}}^{\prime})={\sum}_{k}{\lambda}_{k}{\psi}_{k}(\mathit{\theta}){\psi}_{k}({\mathit{\theta}}^{\prime})$ (Mercer, 1909) (Methods Theory of generalization and Appendix Theory of generalization).
Our analysis shows the existence of a bias in the readout towards learning certain target responses faster than others. The target response $y(\mathit{\theta})={\sum}_{k}{v}_{k}{\psi}_{k}(\mathit{\theta})$ and the learned readout response $f(\mathit{\theta})=\sum _{k}{\hat{v}}_{k}(\mathcal{D}){\psi}_{k}(\mathit{\theta})$ can be expressed in terms of these eigenfunctions ${\psi}_{k}$. Our theory shows that the readout’s generalization is better if the target function $y(\mathit{\theta})$ is aligned with the top eigenfunctions ${\psi}_{k}$, equivalent to ${v}_{k}^{2}$ decaying rapidly with $k$ (Appendix Spectral bias and codetask alignment). We formalize this notion by the following metric. Mathematically, generalization error $\u27e8{E}_{g}\u27e9$ can be decomposed into normalized estimation errors ${E}_{k}$ for the coefficients of these eigenfunctions ${\psi}_{k}$, $\u27e8{E}_{g}\u27e9}_{\mathcal{D}}=\sum _{k}{v}_{k}^{2}{E}_{k$, where $E}_{k}={\u27e8({\hat{v}}_{k}(\mathcal{D}){v}_{k}{)}^{2}\u27e9}_{\mathcal{D}}/{v}_{k}^{2$. We found that the ordering of the eigenvalues ${\lambda}_{k}$ controls the rates at which these mode errors ${E}_{k}$ decrease as $P$ increases (Methods Theory of generalization, Appendix Spectral bias and codetask alignment), (Bordelon et al., 2020): $\lambda}_{k}>{\lambda}_{\ell}\phantom{\rule{thickmathspace}{0ex}}\u27f9\phantom{\rule{thickmathspace}{0ex}}{E}_{k}<{E}_{\ell$. Hence, larger eigenvalues mean lower generalization error for those normalized mode errors ${E}_{k}$. We term this phenomenon the spectral bias of the readout. Based on this observation, we propose codetask alignment as a principle for good generalization. To quantify codetask alignment, we use a metric which was introduced in Canatar et al., 2021 to measure the compatibility of a kernel with a learning task. This is the cumulative power distribution $C(k)$ which measures the total power of the target function in the top $k$ eigenmodes, normalized by the total power (Canatar et al., 2021):
Stimulusresponse maps that have high alignment with the population code’s kernel will have quickly rising cumulative power distributions $C(k)$, since a large proportion of power is placed in the top modes. Target responses with high $C(k)$ can be learned with fewer training samples than target responses with low $C(k)$ since the mode errors ${E}_{k}$ are ordered for all $P$ (Appendix Spectral bias and codetask alignment).
Probing learning biases in neural data
Our theory can be used to probe the learning biases of neural populations. Here, we provide various examples of this using publicly available calcium imaging recordings from mouse primary visual cortex (V1). Our examples illustrate how our theory can be used to analyze neural data.
We first analyzed population responses to static grating stimuli oriented at an angle $\theta $ (Stringer et al., 2021; Pachitariu et al., 2019). We found that the kernel eigenfunctions have sinusoidal shape with differing frequency. The ordering of the eigenvalues and eigenfunctions in Figure 3A (and Figure 3—figure supplement 1) indicates a frequency bias: lower frequency functions of $\theta $ are easier to estimate at small sample sizes.
We tested this idea by constructing two different orientation discrimination tasks shown in Figure 3B and C, where we assign static grating orientations to positive or negative valence with different frequency squarewave functions of $\theta $. We trained the readout using a subset of the experimentally measured neural responses, and measured the readout’s generalization performance. We found that the cumulative power distribution for the low frequency task has a more rapidly rising $C(k)$ (Figure 3B). Using our theory of generalization, we predicted learning curves for these two tasks, which express the generalization error as a function of the number of sampled stimuli $P$. The error for the low frequency task is lower at all sample sizes than the hard, highfrequency task. The theoretical predictions and numerical experiments show perfect agreement (Figure 3B). More intuition can be gained by visualizing the projection of the neural response along the top principal axes (Figure 3C). For the lowfrequency task, the two target values are well separated along the top two axes. However, the highfrequency task is not well separated along even the top four axes (Figure 3C).
Using the same ideas, we can use our theory to get insight into tasks which the V1 population code is illsuited to learn. For the task of identifying mice and birds (Stringer et al., 2018b, Stringer et al., 2018a) the linear rise in cumulative power indicates that there is roughly equal power along all kernel eigenfunctions, indicating that the representation is poorly aligned to this task (Figure 3D).
To illustrate how our approach can be used for different learning problems, we evaluate the ability of linear readouts to reconstruct natural images from neural responses to those images (Figure 4). The ability to reconstruct sensory stimuli from a neural code is an influential normative principle for primary visual cortex (Olshausen and Field, 1997). Here, we ask which aspects of the presented natural scene stimuli are easiest to learn to reconstruct. Since mouse V1 neurons tend to be selective towards low spatial frequency bands (Niell and Stryker, 2008Bonin et al., 2011; Vreysen et al., 2012), we consider reconstruction of bandpass filtered images with spatial frequency wavevector $\mathbf{k}\in {\mathbb{R}}^{2}$ constrained to an annulus $\mathbf{k}\in \left[\sqrt{max({s}_{max}^{2}{r}^{2},0)},{s}_{max}\right]$ for $r=0.2$ (in units of ${\text{pixels}}^{1}$) and plot the cumulative power $C(k)$ associated with each choice of the upper limit ${s}_{max}$ (Figure 4C and D). The frequency cutoffs were chosen in this way to preserve the volume in Fourier space to ${V}_{\mathbf{k}}=\pi {r}^{2}$ for $r<{s}_{max}$, which quantifies the dimension of the function space. We see that the lower frequency bandlimited images are easier to reconstruct, as evidenced by their cumulative power $C(k)$ and learning curves ${E}_{g}$ (Figure 4D and E). This reflects the fact that the population preferentially encodes low spatial frequency content in the image (Figure 4F). Experiments with additional values of $r$ are provided in the Figure 4—figure supplement 1 with additional details found in the Appendix Visual scene reconstruction task.
Mechanisms of spectral bias and codetask alignment in a simple model of V1
How do population level inductive biases arise from properties of single neurons? To illustrate that a variety of mechanisms may be involved in a complex manner, we study a simple model of V1 to elucidate neural mechanisms that lead to the low frequency bias at the population level. In particular, we focus on neural nonlinearities and selectivity profiles.
We model responses of V1 neurons as photoreceptor inputs passed through Gabor filters and a subsequent experimentally motivated powerlaw nonlinearity (Adelson and Bergen, 1985; Olshausen and Field, 1997; Rumyantsev et al., 2020), modeling a population of orientation selective simple cells (Figure 5A) (see Appendix A simple feedforward model of V1). In this model, the kernel for static gratings with orientation $\theta \in [0,\pi ]$ is of the form $K(\theta ,{\theta}^{\prime})=\kappa (\theta {\theta}^{\prime})$, and, as a consequence, the eigenfunctions of the kernel in this setting are Fourier modes. The eigenvalues, and hence the strength of the spectral bias, are determined by the nonlinearity as we discuss in Appendix Gabor model spectral bias and fit to V1 data. We numerically fit the parameters of the nonlinearity to the V1 responses and use these parameters our investigations in Figure 5—figure supplement 1.
Next, to further illustrate the importance of codetask alignment, we study how invariances in the code to stimulus variations may affect the learning performance. We introduce complex cells in addition to simple cells in our model with proportion $s\in [0,1]$ of simple cells (Appendix Gabor model spectral bias and fit to V1 data; Figure 5A), and allow phase, $\varphi $, variations in static gratings. We use the energy model (Adelson and Bergen, 1985; Simoncelli and Heeger, 1998) to capture the phase invariant complex cell responses (Appendix Phase variation, complex cells and invariance, complex cell populations are phase invariant). We reason that in tasks that do not depend on phase information, complex cells should improve sample efficiency.
In this model, the kernel for the V1 population is a convex combination of the kernels for the simple and complex cell populations ${K}_{V1}(\theta ,{\theta}^{\prime},\varphi ,{\varphi}^{\prime})=s{K}_{s}(\theta ,{\theta}^{\prime},\varphi ,{\varphi}^{\prime})+(1s){K}_{c}(\theta ,{\theta}^{\prime})$ where ${K}_{s}$ is the kernel for a pure simple cell population that depends on both orientation and phase, and ${K}_{c}$ is the kernel of a pure complex cell population that is invariant to phase (Appendix Complex cell populations are phase invariant). Figure 5C shows top kernel eigenfunctions for various values of $s$ elucidating inductive bias of the readout.
Figure 5D and E show generalization performance on tasks with varying levels of dependence on phase and orientation. On pure orientation discrimination tasks, increasing the proportion of complex cells by decreasing $s$ improves generalization. Increasing the sensitivity to the nuisance phase variable, $\varphi $, only degrades performance. The cumulative power curve is also maximized at $s=0$. However, on a task which only depends on the phase, a pure complex cell population cannot generalize, since variation in the target function due to changes in phase cannot be explained in the codes’ responses. In this setting, a pure simple cell population attains optimal performance. The cumulative power curve is maximized at $s=1$. Lastly, in a nontrivial hybrid task which requires utilization of both variables $\theta ,\varphi $, an optimal mixture $s$ exists for each sample budget $P$ which minimizes the generalization error. The cumulative power curve is maximized at different $s$ values depending on $k$, the component of the target function. This is consistent with an optimal heterogenous mix, because components of the target are learned successively with increasing sample size. V1 must code for a variety of possible tasks and we can expect a nontrivial optimal simple cell fraction $s$. We conclude that the degree of invariance required for the set of natural tasks, and the number of samples determine the optimal simple cell, complex cell mix. We also considered a more realistic model where the relative selectivity of each visual cortex neuron to phase $\varphi $, measured with the F1/F0 ratio takes on a continuum of possible values with some cells more invariant to phase and some less invariant. In (Appendix Energy model with partially phaseselective cells, Figure 5—figure supplement 3) we discuss a simple adaptation of the energy model which can interpolate between a population of entirely simple cells and a population of entirely complex cells, giving diverse selectivity for the intermediate regime. We show that this model reproduces the inductive bias of Figure 5.
Small and large sample size behaviors of generalization
Recently, Stringer et al., 2018b argued that the inputoutput differentiability of the code, governed by the asymptotic rate of spectral decay, may be enabling better generalization. Our results provide a more nuanced view of the relation between generalization and kernel spectra. First, generalization with low sample sizes crucially depend on the top eigenvalues and eigenfunctions of the code’s kernel, not the tail. Second, generalization requires alignment of the code with the task of interest. Nondifferentiable codes can generalize well if there is such an alignment. To illustrate these points, here, we provide examples where asymptotic conditions on the kernel spectrum are insufficient to describe generalization performance for small sample sizes (Figure 6, Figure 6—figure supplement 1 and Appendix Asymptotic power law scaling of learning curves), and where nondifferentiable kernels generalize better than differentiable kernels (Figure 6—figure supplement 2).
Our example demonstrates how a code allowing good generalization for large sample sizes can be disadvantageous for small sizes. In Figure 6A, we plot three different populations of neurons with smooth (infinitely differentiable) tuning curves that tile a periodic stimulus variable, such as the direction of a moving grating. The tuning width, $\sigma $, of the tuning curves strongly influences the structure of these codes: narrower widths have more high frequency content as we illustrate in a random 3D projection of the population code for $\theta \in [0,2\pi ]$ (Figure 6A). Visualization of the corresponding (von Mises) kernels and their spectra are provided in Figure 6B. The width of the tuning curves control bandwidths of the kernel spectra Figure 6B, with narrower curves having an later decay in the spectrum and higher high frequency eigenvalues. These codes can have dramatically different generalization performance, which we illustrate with a simple “bump" target response (Figure 6C). In this example, for illustration purposes, we let the network learn with a deltarule with a weight decay, leading to a regularized kernel regression solution (Appendix Weight decay and ridge regression). For a sample size of $P=10$, we observe that codes with too wide or too narrow tuning curves (and kernels) do not perform well, and there is a wellperforming code with an optimal tuning curve width $\sigma $, which is compatible with the width of the target bump, ${\sigma}_{T}$. We found that optimal $\sigma $ is different for each $P$ (Figure 6C). In the large$P$ regime, the ordering of the performance of the three codes are reversed (Figure 6C). In this regime generalization error scales in a power law (Appendix Asymptotic power law scaling of learning curves) and the narrow code, which performed worst for $P\sim 10$, performs the best. This example demonstrates that asymptotic conditions on the tail of the spectra are insufficient to understand generalization in the small sample size limit. The bulk of the kernel’s spectrum needs to match the spectral structure of the task to generalize efficiently in the lowsample size regime. However, for large sample sizes, the tail of the eigenvalue spectrum becomes important. We repeat the same exercise and draw the same conclusions for a nondifferentiable kernel (Laplace) (Figure 6—figure supplement 1) showing that these results are not an artifact of the infinite differentiability of von Mises kernels. We further provide examples where nondifferentiable kernels generalizing better than differentiable kernels in Figure 6—figure supplement 2.
Timedependent neural codes
Our framework can directly be extended to learning of arbitrary timevarying functions of timevarying inputs from an arbitrary spatiotemporal population code (Methods RNN experiment, Appendix Time dependent neural codes). In this setting, the population code $\mathbf{r}(\{\mathit{\theta}(t)\},t)$ is a function of an input stimulus sequence $\mathit{\theta}(t)$ and possibly its entire history, and time $t$. A downstream linear readout $f(\{\mathit{\theta}\},t)=\mathbf{w}\cdot \mathbf{r}(\{\mathit{\theta}\},t)$ learns a target sequence $y(\{\mathit{\theta}\},t)$ from a total of $\mathcal{P}$ examples that can come at any time during any sequence.
As a concrete example, we focus on readout from a temporal population code generated by a recurrent neural network in a task motivated by a delayed reach task (Ames et al., 2019; Figure 7A and B). In this task, the network is presented for a short time an input cue sequence coding an angular variable which is drawn randomly from a distribution (Figure 7C). The recurrent neural network must remember this angle and reproduce an output sequence which is a simple step function whose height depends on the angle which begins after a time delay from the cessation of input stimulus and lasts for a short time (Figure 7D).
The kernel induced by the spatiotemporal code is shown in Figure 7E. The high dimensional nature of the activity in the recurrent network introduces complex and rich spatiotemporal similarity structure. Figure 7F shows the kernel’s eigensystem, which consists of stimulus dependent timeseries ${\psi}_{k}(\{\mathit{\theta}\};t)$ for each eigenvalue ${\lambda}_{k}$. An interesting link can be made with this eigensystem and linear lowdimensional manifold dynamics observed in several cortical areas (Stopfer et al., 2003; Kato et al., 2015; Gallego et al., 2017; Cunningham and Yu, 2014; Sadtler et al., 2014; Gao and Ganguli, 2015; Gallego et al., 2018; Chapin and Nicolelis, 1999; Bathellier et al., 2008). The kernel eigenfunctions also define the latent variables obtained through a singular value decomposition of the neural activity $\mathbf{r}(\{\mathit{\theta}\};t)={\sum}_{k}\sqrt{{\lambda}_{k}}{\mathbf{u}}_{k}{\psi}_{k}(\{\mathit{\theta}\};t)$ (Gallego et al., 2017). With enough samples, the readout neuron can learn to output the desired angle with high fidelity (Figure 7G). Unsurprisingly, tasks involving long time delays are more difficult and exhibit lower cumulative power curves. Consequently, the generalization error for small delay tasks drops much more quickly with increasing samples $\mathcal{P}$.
Biological codes are metabolically more efficient and more selective than other codes with identical kernels
Although, the performance of linear readouts may be invariant to rotations that preserve kernels (Figure 2), metabolic efficiency may favor certain codes over others (Barlow, 1961; Atick and Redlich, 1992; Attneave, 1954; Olshausen and Field, 1997; Simoncelli and Olshausen, 2001), reducing degeneracy in the space of codes with identical kernels. To formalize this idea, we define $\mathit{\delta}$ to be the vector of spontaneous firing rates of a population of neurons, and ${\mathbf{s}}^{\mu}=\mathbf{r}({\mathit{\theta}}^{\mu})+\mathit{\delta}$ be the spiking rate vector in response to a stimulus ${\mathit{\theta}}^{\mu}$. The vector $\mathit{\delta}$ ensures that neural responses are nonnegative. The modulation with respect to the spontaneous activity, $\mathbf{r}({\mathit{\theta}}^{\mu})$, gives the population code and defines the kernel, $K({\mathit{\theta}}^{\mu},{\mathit{\theta}}^{\mu})=\frac{1}{N}\mathbf{r}({\mathit{\theta}}^{\mu})\cdot \mathbf{r}({\mathit{\theta}}^{\nu})$. To avoid confusion with $\mathbf{r}({\mathit{\theta}}^{\mu})$, we will refer to ${\mathbf{s}}^{\mu}$ as total spiking activity. We propose that population codes prefer smaller spiking activity subject to a fixed kernel. In other words, because the kernel is invariant to any change of the spontaneous firing rates and left rotations of $\mathbf{r}(\mathit{\theta})$, the orientation and shift of the population code $\mathbf{r}(\mathit{\theta})$ should be chosen such that the resulting total spike count ${\sum}_{i=1}^{N}{\sum}_{\mu =1}^{P}{s}_{i}^{\mu}$ is small.
We tested whether biological codes exhibit lower total spiking activity than others exhibiting the same kernel on mouse V1 recordings, using deconvolved calcium activity as a proxy for spiking events (Stringer et al., 2021; Pachitariu et al., 2019; Pachitariu et al., 2018) (Methods Data analysis; Figure 8). To compare the experimental total spiking activity to other codes with identical kernels, we computed random rotations of the neural responses around spontaneous activity, $\stackrel{~}{\mathbf{r}}({\mathit{\theta}}^{\mu})=\mathbf{Qr}({\mathit{\theta}}^{\mu})$, and added the $\delta $ that minimizes total spiking activity and maintains its nonnegativity (Methods Generating RROS codes). We refer to such an operation as RROS (random rotation and optimal shift), and a code generated by an RROS operation as an RROS code. The matrix $\mathbf{Q}$ is a randomly sampled orthogonal matrix (Anderson et al., 1987). In other words, we compare the true code to the most metabolically efficient realizations of its random rotations. This procedure may result in an increased or decreased total spike count in the code, and is illustrated in a synthetic dataset in Figure 8A. We conducted this procedure on subsets of various sizes of mouse V1 neuron populations, as our proposal should hold for any subset of neurons (Methods Generating RROS codes), and found that the true V1 code is much more metabolically efficient than randomly rotated versions of the code (Figure 8B and C). This finding holds for both responses to static gratings and to natural images as we show in Figure 8B and C respectively.
To further explore metabolic efficiency, we posed an optimization problem which identifies the most efficient code with the same kernel as the biological V1 code. This problem searches over rotation matrices $\mathbf{Q}$ and finds the $\mathbf{Q}$ matrix and offset vector $\mathit{\delta}$ which gives the lowest cost ${\sum}_{i\mu}{s}_{i}^{\mu}$ (Methods Comparing sparsity of population codes) (Figure 8). Although the local optimum identified with the algorithm is lower in cost than the biological code, both the optimal and biological codes are significantly displaced from the distribution of random codes with same kernel. Our findings do not change when data is preprocessed with an alternative strategy, an upper bound on neural responses is imposed on rotated codes, or subsets of stimuli are considered (Figure 8—figure supplement 1). We further verified these results on electrophysiological recordings of mouse visual cortex from the Allen Institute Brain Observatory (de Vries et al., 2020), (Figure 8—figure supplement 2). Overall, the large disparity in total spiking activity between the true and randomly generated codes with identical kernels suggests that metabolic constraints may favor the biological code over others that realize the same kernel.
The disparity between the true biological code and the RROS code is not only manifested in terms of total activity level, but also in terms of single neuron and single stimulus sparseness measures, specifically lifetime and population sparseness distributions (Methods Lifetime and population sparseness) (Willmore and Tolhurst, 2001; Lehky et al., 2005; Treves and Rolls, 1991; Pehlevan and Sompolinsky, 2014). In Figure 8D, we compare the lifetime and population sparseness distributions of the true biological code with a RROS version of the same code, revealing biological neurons have significantly higher lifetime sparseness. In Appendix Necessary conditions for optimally sparse codes, we provide analytical arguments which suggest that tuning curves of optimally sparse nonnegative codes with fullrank kernels will have selective tuning.
Discussion
Elucidating inductive biases of the brain is fundamentally important for understanding natural intelligence (Tenenbaum et al., 2011; Lake et al., 2017; Sinz et al., 2019; Zador, 2019). These biases are coded into the brain by the dynamics of its neurons, the architecture of its networks, its representations and plasticity rules. Finding ways to extract the inductive biases from neuroscience datasets requires a deeper theoretical understanding of how all these factors shape the biases, and is an open problem. In this work, we attempted to take a step towards filling this gap by focusing on how the structure of static neural population codes shape inductive biases for learning of a linear readout neuron under a biologically plausible learning rule. If the readout neuron’s output is correlated with behavior, and that correlation is known, then our theory could possibly be modified to predict what behavioral tasks can be learned faster.
Under the delta rule, the generalization performance of the readout is entirely dependent on the code’s inner product kernel; the kernel is a determinant of inductive bias. In its finite dimensional form, the kernel is an example of a representational similarity matrix and is a commonly used tool to study neural representations (Edelman, 1998; Kriegeskorte et al., 2008; Laakso and Cottrell, 2000; Kornblith et al., 2019; Cadieu et al., 2014; Pehlevan et al., 2018). Our work elucidates a concrete link between this experimentally measurable mathematical object, and sampleefficient learning.
We derived an analytical expression for the generalization error as a function of samplesize under very general conditions, for an arbitrary stimulus distribution, arbitrary population code and an arbitrary target stimulusresponse map. We used our findings in both theoretical and experimental analysis of primary visual cortex, and temporal codes in a delayed reach task. This generality of our theory is a particular strength.
Our analysis elucidated two principles that define the inductive bias. The first one is spectral bias: kernel eigenfunctions with large eigenvalues can be estimated using a smaller number of samples. The second principle is the codetask alignment: target functions with most of their power in top kernel eigenfunctions can be estimated efficiently and are compatible with the code. The cumulative power distribution, $C(k)$ (Canatar et al., 2021), provides a measure of this alignment. These findings define a notion of ‘simplicity’ bias in learning from examples, and provides a solution to the question of what stimulusresponse maps are easier to learn. A similar simplicity bias has been also observed in training deep neural networks (Rahaman et al., 2019; Xu et al., 2019; Kalimeris et al., 2019). Due to a correspondence between gradientdescent trained neural networks in the infinitewidth limit and kernel machines (Jacot et al., 2018), results on the spectral bias of kernel machines may shed light onto these findings (Bordelon et al., 2020; Canatar et al., 2021). Though our present analysis focused on learning a single layer weight vector with the biologically plausible deltarule, future work could explore the learning curves of other learning rules for deep networks (Bordelon and Pehlevan, 2022a), such as feedback alignment (Lillicrap et al., 2016) or perturbation methods (Jabri and Flower, 1992). Such analysis could explore how inductive bias is also shaped by choice of learning rule, as well as the structure of the initial population code.
We applied our findings in both theoretical and experimental analysis of mouse primary visual cortex. We demonstrated a bias of neural populations towards low frequency orientation discrimination and low spatial frequency reconstruction tasks. The latter finding is consistent with the finding that mouse visual cortex neurons are selective for low spatial frequency (Niell and Stryker, 2008; Vreysen et al., 2012). The toy model of the visual cortex as a mixture of simple and complex cells demonstrated how invariances, specifically the phase invariance of the complex cells, in the population code can facilitate learning some tasks involving phase invariant responses at the expense of performance on others. The role of invariances in learning with kernel methods and deep networks have recently been investigated in machine learning literature, showing that invariant representations can improve capacity (Farrell et al., 2021) and sample efficiency for invariant learning problems (Mei et al., 2021; Li et al., 2019; Xiao and Pennington, 2022).
A recent proposal considered the possibility that the brain acts as an overparameterized interpolator (Hasson et al., 2020). Suitable inductive biases are crucial to prevent overfitting and generalize well in such a regime (Belkin et al., 2019). Our theory can explain these inductive biases since, when the kernel is fullrank, which typically is the case when there are more neurons in the population than the number of learning examples, the delta rule without weight decay converges to an interpolator of the learning examples. Modern deep learning architectures also operate in an overparameterized regime, but generalize well (Zhang et al., 2016; Belkin et al., 2019), and an inductive bias towards simple functions has been proposed as an explanation (Bordelon et al., 2020; Canatar et al., 2021; Kalimeris et al., 2019; VallePerez et al., 2018). However, we also showed that interpolation can be harmful to prediction accuracy when the target function has some variance unexplained by the neural code or if the neural responses are significantly noisy, motivating use of explicit regularization.
Our work promotes sample efficiency as a general coding principle for neural populations, relating neural representations to the kinds of problems they are well suited to solve. These codes may be shaped through evolution or themselves be learned through prior experience (Zador, 2019). Prior related work in this area demonstrated the dependence of sampleefficient learning of a twoangle estimation task on the width of the individual neural tuning curves (Meier et al., 2020) and on additive function approximation properties of sparsely connected random networks (Harris, 2019).
A sample efficiency approach to population coding differs from the classical efficient coding theories (Attneave, 1954; Barlow, 1961; Atick and Redlich, 1992; Srinivasan et al., 1982; van Hateren, 1992; Rao and Ballard, 1999; Olshausen and Field, 1997; Chalk et al., 2018), which postulate that populations of neurons optimize information content of their code subject to metabolic constraints or noise. While these theories emphasize different aspects of the code’s information content (such as reduced redundancy, predictive power, or sparsity), they do not address sample efficiency demands on learning. Further, recent studies demonstrated hallmarks of redundancy and correlation in population responses (Chapin and Nicolelis, 1999; Bathellier et al., 2008; Pitkow and Meister, 2012; Gao and Ganguli, 2015; AbbasiAsl et al., 2016; Gallego et al., 2018; Stringer et al., 2018b), violating a generic prediction of efficient coding theories that responses of different neurons should be uncorrelated across input stimuli in high signaltonoise regimes to reduce redundancy in the code and maximize information content (Barlow, 1961; Atick and Redlich, 1992; Srinivasan et al., 1982; van Hateren, 1992; Haft and van Hemmen, 1998; Huang and Rao, 2011). In our theory, the structured correlations of neural responses correspond to the decay in the spectrum of the kernel, and play a key role in biasing learned readouts towards simple functions.
In recent related studies, the asymptotic decay rate of the kernel’s eigenspectrum was argued to be important for generalization (Stringer et al., 2018b) and robustness (Nassar et al., 2020). The spectral decay rate in the mouse V1 was found to be consistent with a high dimensional (power law) but smooth (differentiable) code, and smoothness was argued to be an enabler of generalization (Stringer et al., 2018b). While we also identify power law spectral decays, we show that sampleefficient learning requires more than smoothness conditions in the form of asymptotic decay rates on the kernel’s spectrum. The interplay between the stimulus distribution, target response and the code gives rise to sample efficient learning. Because of spectral bias, the top eigenvalues govern the small sample size behavior. The tail of the spectrum becomes important at large sample sizes.
Though the kernel is degenerate with respect to rotations of the code in the neural activity space, we demonstrated that the true V1 code has much lower average activity than random codes with the same kernel, suggesting that evolution and learning may be selecting neural codes with low average spike rates which preserve sampleefficiency demands for downstream learning tasks. We predict that metabolic efficiency may be a determinant in the orientation and placement of the ubiquitously observed lowdimensional coding manifolds in neural activity space in other parts of the brain (Gallego et al., 2018). The demand of metabolic efficiency is consistent with prior sparse coding theories (Niven and Laughlin, 2008; Olshausen and Field, 1997; Simoncelli and Olshausen, 2001; Hromádka et al., 2008), however, our theory emphasizes sampleefficient learning as the primary normative objective for the code. As a note of caution, while our analysis holds under the assumption that the neural code is deterministic, real neurons exhibit variability in their responses to repeated stimuli. Such noisy population codes do not generally achieve identical generalization performance under RROS transformations. For example, if each neuron is constrained to produce i.i.d. Poisson noise, then simple shifts of the baseline firing rate reduce the information content of the code. However, if the neural noise is Gaussian (even with stimulus dependent noise covariance), then the generalization error is conserved under RROS operations (Appendix Effect of noise on RROS symmetry). Further studies could focus on revealing the space of codes with equivalent inductive biases under realistic noise models.
Our work constitutes a first step towards understanding inductive biases in neuronal circuits. To achieve this, we focused on a linear, deltarule readout of a static population code. More work is need to study other factors that affect inductive bias. Importantly, sensory neuron tuning curves can adapt during perceptual learning tasks (Gilbert, 1994; Goltstein et al., 2021; Ghose et al., 2002; Schoups et al., 2001) with the strength of adaptation dependent on brain area (Yang and Maunsell, 2004; Adab et al., 2014; Op de Beeck et al., 2007). In many experiments, these changes to tuning in sensory areas are small (Schoups et al., 2001; Ghose et al., 2002), satisfying the assumptions of our theory. For example monkeys trained on noisy visual motion detection exhibit changes in sensorymotor (LIP) but not sensory areas (MT), consistent with a model of readout from a static sensory population code (Law and Gold, 2008; Shadlen and Newsome, 2001). However, other perceptual learning tasks and other brain areas can exhibit significant changes in neural tuning (Recanzone et al., 1993; Pleger et al., 2003; Furmanski et al., 2004). This diversity of results motivates more general analysis of learning in multilayer networks, where representations in each layer can adapt flexibly to task structure (Shan and Sompolinsky, 2021; Mastrogiuseppe et al., 2022; Bordelon and Pehlevan, 2022b; Ahissar and Hochstein, 2004). Alternatively, our current analysis of inductive bias can still be consistent with multilayer learning if the network is sufficiently overparameterized and tuning curves change very little (Jacot et al., 2018; Lee et al., 2018; Shan and Sompolinsky, 2021). In this case, network training is equivalent to kernel learning with a kernel that depends on the learning rule and architecture (Bordelon and Pehlevan, 2022a). However, in the regime of neural network training where tuning curves change significantly, more sophisticated analytical tools are needed to predict generalization (Flesch et al., 2021; Yang and Hu, 2021; Bordelon and Pehlevan, 2022b). Although our work focused on linear readouts, arbitrary nonlinear readouts which generate convex learning objectives have been recently studied in the high dimensional limit, giving qualitatively similar learning curves which depend on kernel eigenvalues and task model alignment (Loureiro et al., 2021b; Cui et al., 2022) (see Appendix Typical case analysis of nonlinear readouts).
Our work focused on how signal correlations influence inductive bias (Averbeck et al., 2006; Cohen and Kohn, 2011). However, since real neurons do exhibit variability in their responses to identical stimuli, one should consider the effect of neural noise and noise correlations in learning. We provide a preliminary analysis of learning with neural noise in Appendix Impact of neural noise and unlearnable targets on learning, where we show that neural noise can lead to irreducible asymptotic error which depends on the geometry of the signal and noise correlations. Further, if the target function is not fully expressible as linear combinations of neural responses, overfitting peaks in the learning curves are possible, but can be mitigated with regularization implemented by a weight decay in the learning rule (see Appendix 1—figure 1). Future work could extend our analysis to study how signal and noise correlations interact to shape inductive bias and generalization performance in the case where the noise correlation matrices are nonisotropic, including the role of differential correlations (MorenoBote et al., 2014). Overall, future work could build on the present analysis to incorporate a greater degree of realism in a theory of inductive bias.
Finally, we discuss possible applications of our work to experimental neuroscience. Our theory has potential implications for experimental studies of task learning. First, in cases where the population selective to stimuli can be measured directly, an experimenter could design easy or difficult tasks for an animal to learn from few examples, under a hypothesis that the behavioral output is a linear function of the observed neurons. Second, in cases where it is unclear which neural population contributes to learning, one could utilize our theory to solve the inverse problem of inferring the relevant kernel from observed learning curves on different tasks (Wilson et al., 2015). From these tasks, the experimenter could compare the inferred kernel to those of different recorded populations. For instance, one could compare the kernels from separate populations to the inferred kernel obtained from learning curves on certain visual learning tasks. This could provide new ways to test theories of perceptual learning (Gilbert, 1994). Lastly, extensions of our framework could quantify the role of neural variability on task learning and the limitation it imposes on accuracy and sample efficiency.
Methods
Generating example codes (Figure 1)
The two codes in Figure 1 were constructed to produce two different kernels for $\theta \in {S}^{1}$:
An infinite number of codes could generate either of these kernels. After diagonalizing the kernel into its eigenfunctions on a grid of 120 points, $\mathbf{K}}_{1}={\mathbf{\Psi}}_{1}{\mathbf{\Lambda}}_{1}{\mathbf{\Psi}}_{1}^{\mathrm{\top}},{\mathbf{K}}_{2}={\mathbf{\Psi}}_{2}{\mathbf{\Lambda}}_{2}{\mathbf{\Psi}}_{2}^{\mathrm{\top}$, we used a random rotation matrix $\mathbf{Q}\in {\mathbb{R}}^{N\times N}$ (which satisfies ${\mathbf{QQ}}^{\top}={\mathbf{Q}}^{\top}\mathbf{Q}=\mathbf{I}$) to generate a valid code
This construction guarantees that ${\mathbf{R}}_{1}^{\top}{\mathbf{R}}_{1}={\mathbf{K}}_{1}$ and ${\mathbf{R}}_{2}^{\top}{\mathbf{R}}_{2}={\mathbf{K}}_{2}$. We plot the tuning curves for the first three neurons. The target function in the first experiment is $y=\mathrm{cos}(\theta )0.6\mathrm{cos}(4\theta )$, while the second experiment used $y=\mathrm{cos}(6\theta )\mathrm{cos}(8\theta )$.
Theory of generalization
Recent work has established analytic results that predict the average case generalization error for kernel regression
where $E}_{g}(\mathcal{D})={\u27e8[f(\mathit{\theta},\mathcal{D})y(\mathit{\theta}){]}^{2}\u27e9}_{\mathit{\theta}$ is the generalization error for a certain sample $\mathcal{D}$ of size $P$ and $f(\mathit{\theta},\mathcal{D})=\mathbf{w}\cdot \mathbf{r}(\mathit{\theta})$ is the kernel regression solution for $\mathcal{D}$ (Appendix Convergence of the deltarule without weight decay) (Bordelon et al., 2020; Canatar et al., 2021). The typical or average case error ${E}_{g}$ is obtained by averaging over all possible datasets of size $P$. This average case generalization error is determined solely by the decomposition of the target function $y(\mathbf{x})$ along the eigenbasis of the kernel and the eigenspectrum of the kernel. This continuous diagonalization again takes the form (Appendix Singular value decomposition of continuous population responses) (Rasmussen and Williams, 2005)
Our theory is also applicable to discrete stimuli if $p(\mathit{\theta})$ is a Dirac measure as we describe in (Appendix Discrete stimulus spaces: finding eigenfunctions with matrix eigendecomposition). Since the eigenfunctions form a complete set of square integrable functions (Rasmussen and Williams, 2005), we expand both the target function $y(\mathit{\theta})$ and the learned function $f(\mathit{\theta},\mathcal{D})$ in this basis $y(\mathit{\theta})=\sum _{k}{v}_{k}{\psi}_{k}(\mathit{\theta})\text{},\text{}f(\mathit{\theta},\mathcal{D})=\sum _{k}{\hat{v}}_{k}{\psi}_{k}(\mathit{\theta})$, where ${\widehat{v}}_{k}$ are understood to be functions of the dataset $\mathcal{D}$. The eigenfunctions are orthonormal $\int d\mathit{\theta}p(\mathit{\theta}){\psi}_{k}(\mathit{\theta}){\psi}_{\mathrm{\ell}}(\mathit{\theta})={\delta}_{k,\mathrm{\ell}}$, which implies that the generalization error for any set of coefficients $\widehat{\mathbf{v}}$ is
We now introduce the equivalent training error, or empirical loss, written directly in terms of eigenfunction coefficients $\widehat{\mathbf{v}}$, which depends on the random dataset $\mathcal{D}=\{({\mathit{\theta}}^{\mu},{y}^{\mu}){\}}_{\mu =1}^{P}$
This loss function is minimized by delta rule updates with weight decay constant $\lambda $. It is straightforward to verify that the $H$minimizing coefficients are $\hat{\mathbf{v}}}^{\ast}=(\mathbf{\Psi}{\mathbf{\Psi}}^{\mathrm{\top}}+\lambda {\mathbf{\Lambda}}^{1}{)}^{1}\mathbf{\Psi}{\mathbf{\Psi}}^{\mathrm{\top}}\mathbf{v$, giving the learned function $f(\mathit{\theta},\mathcal{D})={\hat{\mathbf{v}}}^{\ast}\cdot \mathit{\psi}(\mathit{\theta})=\mathbf{k}(\mathit{\theta}{)}^{\mathrm{\top}}(\mathbf{K}+\lambda \mathbf{I}{)}^{1}\mathbf{y}$ where the vectors $\mathbf{k}$ and $\mathbf{y}$ have entries $[\mathbf{k}(\mathit{\theta}){]}_{\mu}=K(\mathit{\theta},{\mathit{\theta}}^{\mu})$ and $[\mathbf{y}{]}_{\mu}=y({\mathit{\theta}}^{\mu})$ for each training stimulus $\mathit{\theta}}^{\mu}\in \mathcal{D$. The $P\times P$ kernel gram matrix $\mathbf{K}$ has entries $[\mathbf{K}{]}_{\mu \nu}=K({\mathit{\theta}}^{\mu},{\mathit{\theta}}^{\nu})$. The $\lambda \to 0$ limit of the minimizer of $H$ coincides with kernel interpolation. This allows us to characterize generalization without reference to learned readout weights $\mathbf{w}$. The generalization error for this optimal function is
We note that the dependence on the randomly sampled dataset $\mathcal{D}$ only appears through the matrix $\mathbf{G}(\mathcal{D})$. Thus to compute the typical generalization error we need to average this matrix over realizations of datasets, i.e.$\u27e8\mathbf{G}(\mathcal{D})\u27e9}_{\mathcal{D}$. There are multiple strategies to perform such an average and we will study one here based on a partial differential equation which was introduced in Sollich, 1998; Sollich, 2002 and studied further in Bordelon et al., 2020. We describe in detail one method for performing such an average in Appendix Computation of learning curves. After this computation, we find that the generalization error can be approximated at large $P$ as
where $\gamma =P{\sum}_{k}\frac{{\lambda}_{k}^{2}}{{({\lambda}_{k}p+\kappa )}^{2}}$, giving the desired result. We note that (11) defines the function $\kappa $ implicitly in terms of the sample size $P$. Taking $\lambda \to 0$ gives the generalization error of the minimum norm interpolant, which desribes the generalization error of the solution. This result was recently reproduced using the replica method from statistical mechanics in an asymptotic limit where the number of neurons and samples are large (Bordelon et al., 2020; Canatar et al., 2021). Other recent works have verified our theoretical expressions on a variety of kernels and datasets (Loureiro et al., 2021b; Simon et al., 2021).
Additional intuition for the spectral bias phenomenon can be gained from the expected learned function ${\u27e8f(\mathit{\theta},\mathcal{D})\u27e9}_{\mathcal{D}}=\sum _{k}\frac{{\lambda}_{k}P}{{\lambda}_{k}P+\kappa}{v}_{k}{\psi}_{k}(\mathit{\theta})$, which is the average readout prediction over possible datasets $\mathcal{D}$. The function $\kappa (P)$ is defined implicitly as $\kappa =\lambda +\kappa {\sum}_{k}\frac{{\lambda}_{k}}{{\lambda}_{k}P+\kappa}$ and decreases with $P$ from $\kappa (0)=\lambda +{\sum}_{k}{\lambda}_{k}$ to its asymptotic value ${lim}_{P\to \mathrm{\infty}}\kappa (P)=\lambda $. The coefficient for the $k$th eigenfunction $\frac{{\lambda}_{k}P}{{\lambda}_{k}P+\kappa}{v}_{k}$ approaches the true coefficient v_{k} as $P\to \mathrm{\infty}$. The $k$th eigenfunction ${\psi}_{k}$ is effectively learned when $P\gg \frac{\kappa}{{\lambda}_{k}}$. For large eigenvalues ${\lambda}_{k}$, fewer samples $P$ are needed to satisfy this condition, while small eigenvalue modes will require more samples.
RNN experiment
For the simulations in Figure 7, we integrated a ratebased recurrent network model with $N=6000$ neurons, time constant $\tau =0.05$ and gain $g=1.5$. Each of the $P=80$ randomly chosen angles ${\gamma}^{\mu}$ generates a trajectory over $T=100$ equally spaced points in $t\in [0,3]$. The two dimensional input sequence is simply $\mathit{\theta}(t)=H(t)H(1t){[\mathrm{cos}({\gamma}^{\mu}),\mathrm{sin}({\gamma}^{\mu})]}^{\top}\in {\mathbb{R}}^{2}$. Target function for a delay $d$ is $\mathbf{y}({\theta}^{\mu},t)=H(1.5+dt)H(td1){[\mathrm{cos}({\gamma}^{\mu}),\mathrm{sin}({\gamma}^{\mu})]}^{\top}$ which is nonzero for times $t\in [1+d,1.5+d]$. In each simulation, the activity in the network is initialized to $\mathbf{u}(0)=\mathbf{0}$. The kernel gram matrix $\mathbf{K}\in {\mathbb{R}}^{PT\times PT}$ is computed by taking inner products of the time varying code at for different inputs ${\gamma}^{\mu}$ and at different times. Learning curves represent the generalization error obtained by randomly sampling $\mathcal{P}$ time points from the $PT$ total time points generated in the simulation process and training readout weights $\mathbf{w}$ to convergence with gradient descent.
Data analysis
Data source and processing
Mouse V1 neuron responses to orientation gratings were obtained from a publicly available dataset (Stringer et al., 2021; Pachitariu et al., 2019). Twophoton calcium microscopy fluorescence traces were deconvolved into spike trains and spikes were counted for each stimulus, as described in Stringer et al., 2021. The presented grating angles were distributed uniformly over $[0,2\pi ]$ radians. Data preprocessing, which included zscoring against the mean and standard deviation of null stimulus responses, utilized the provided code for this experiment, which also publicly available at https://github.com/MouseLand/stringeretal2019 (Stringer, 2019). This preprocessing technique was used in all Figures in the paper. To reduce corruption of the estimated kernel from neural noise (trialtotrial variability), we first trial average responses, binning the grating stimuli oriented at different angles $\theta $ into a collection of 100 bins over the interval from $[0,2\pi ]$ and averaging over all of the available responses from each bin. Since grating angles were sampled uniformly, there is a roughly even distribution of about 45 responses in each bin. After trial averaging, SVD was performed on the response matrix $\mathbf{R}$, generating the eigenspectrum and kernel eigenfunctions as illustrated in Figure 3. Figures 2, 3 and 8, all used this data anytime responses to grating stimuli were mentioned.
In Figures 3D, 4 and 8C, the responses of mouse V1 neurons to ImageNet images (Deng et al., 2009) were obtained from a different publicly available dataset (Stringer et al., 2018a). The images were taken from 15 different classes from the Imagenet dataset with ethological relevance to mice (birds, cats, flowers, hamsters, holes, insects, mice, mushrooms, nests, pellets, snakes, wildcats, other animals, other natural, other man made). In the experiment in Figure 3D we use all images from the mice and birds category for which responses were recorded. The preprocessing code and image category information were obtained from the publicly available code base at https://github.com/MouseLand/stringerpachitariuetal2018b (Stringer, 2018c). Again, spike counts were obtained from deconvolved and zscored calcium fluorescence traces. In the reconstruction experiment shown in Figure 4 we use the entire set of images for which neural responses were recorded.
Generating RROS codes
In Figure 8, the randomly rotated codes are generated by sampling a matrix $\mathbf{Q}$ from the Haar measure on the set of $N$by$N$ orthogonal matrices (Anderson et al., 1987), and chosing a $\mathit{\delta}$ by solving the following optimization problem:
which minimizes the total spike count subject to the kernel and nonnegativity of firing rates. The solution to this problem is given by ${\delta}_{i}^{*}={\mathrm{min}}_{\mu =1,\mathrm{\dots},P}{[\mathbf{Qr}({\mathit{\theta}}^{\mu})]}_{i}$.
Comparing sparsity of population codes
To explore the metabolic cost among the set of codes with the same inductive biases, we estimate the distribution of average spike counts of codes with the same inner product kernel as the biological code. These codes are generated in the form ${\mathbf{s}}^{\mu}={\mathbf{Qr}}^{\mu}+\mathit{\delta}$ where $\mathit{\delta}$ solves the optimization problem
To quantify the distribution of such codes, we randomly sample $\mathbf{Q}$ from the invariant (Haar) measure for $N\times N$ orthogonal matrices and compute the optimal $\mathit{\delta}$ as described above. This generates the aqua colored distribution in Figure 8B and C.
We also attempt to characterize the most efficient code with the same inner product kernel
Since this optimization problem is nonconvex in $\mathbf{Q}$, there is no theoretical guarantee that minima are unique. Nonetheless, we attempt to optimize the code by starting $\mathbf{Q}$ at the identity matrix and conduct gradient descent over orthogonal matrices (Plumbley, 2004). Such updates take the form
where $\mathrm{exp}(\cdot )$ is the matrix exponential. To make the loss function differentiable, we incorporate the nonnegativity constraint with a softminimum:
where $Z=\sum _{\nu}\mathrm{exp}(\beta {a}^{\nu})$ is a normalizing constant and $\mathbf{Q}=[{\mathbf{q}}_{1},\mathrm{\dots}{\mathbf{q}}_{N}]$. In the $\beta \to \mathrm{\infty}$ limit, this cost function converges to the exact optimization problem with nonnegativity constraint. Finite $\beta $, however, allows learning with gradient descent. Gradients are computed with automatic differentiation in JAX (Bradbury et al., 2018). This optimization routine is run until convergence and the optimal value is plotted as dashed red lines labeled ‘opt’. in Figure 8.
We show that our result is robust to different preprocessing techniques and to imposing bounds on neural firing rates in the Figure 8—figure supplement 1. To demonstrate that our result is not an artifact of zscoring the deconvolved signals against the spontaneous baseline activity level, we also conduct the random rotation experiment on the raw deconvolved signals. In addition, we show that imposing realistic constraints on the upper bound of the each neuron’s responses does not change our findings. We used a subset of $N=100$ neurons and computed random rotations. However, we only accepted a code as valid if it’s maximum value was less than some upper bound u_{b}. Subsets of $N=100$ neurons in the biological code achieve maxima in the range between 3.2 and 4.7. We performed this experiment for ${u}_{b}\in \{3,4,5\}$ so that the artificial codes would have maxima that lie in the same range as the biological code.
Lifetime and population sparseness
We compute two more refined measures of sparseness in a population code. For each neuron $i$ we compute the lifetime sparseness $L{S}_{i}$ (also known as selectivity) and for each stimulus $\mathit{\theta}$ we compute the population sparseness $P{S}_{\mathit{\theta}}$ which are defined as the following two ratios (Willmore and Tolhurst, 2001; Lehky et al., 2005; Treves and Rolls, 1991; Pehlevan and Sompolinsky, 2014)
The normalization factors ensure that these quantities lie in the interval between $(0,1)$. Intuitively, lifetime sparseness quantifies the variability of each neuron’s responses over the full set of stimuli, whereas population sparseness quantifies the variability of responses in the code for a given stimulus $\mathit{\theta}$.
Fitting a Gabor model to mouse V1 kernel
Under the assumption of translation symmetry in the kernel $K(\theta ,{\theta}^{\prime})$, we averaged the elements of the over rows of the empirical mouse V1 kernel (Pachitariu et al., 2019)
where angular addition is taken mod $\pi $. This generates the black dots in Figure 5B. We aimed to fit a thresholdpower law nonlinearity of the form ${g}_{q,a}(z)=\mathrm{max}{\{0,za\}}^{q}$ to the kernel. Based on the Gabor model discussed above, we parameterized tuning curves as
where ${\theta}_{i}$ is the preferred angle of the $i$th neuron’s tuning curve. Rather than attempting to perform a fit of ${\sigma}^{2},a,q,{\{{\theta}_{i}\}}_{i=1}^{N}$ of this form to the responses of each of the $\sim 20$k neurons, we instead simply attempt to fit to the population kernel by optimizing over $({\sigma}^{2},a,q)$ under the assumption of uniform tiling of ${\theta}_{i}$. However, we noticed that two of these variables ${\sigma}^{2},a$ are constrained by the sparsity level of the code. If each neuron, on average, fires for only a fraction $f$ of the uniformly sampled angles $\theta $, then the following relationship holds between ${\sigma}^{2}$ and
Calculation of the coding level $f$ for the recorded responses allowed us to infer $a$ from $\sigma}^{2$ during optimization. This reduced the free parameter set to $({\sigma}^{2},q)$. We then solve the following optimization problem
where integration over ${\theta}_{i}$ is performed numerically. Using the Scipy TrustRegion constrained optimizer, we found $(q,{\sigma}^{2},a)=(1.7,5.0,0.2)$ which we use as the fit parameters in Figure 5.
Lead contact
Requests for information should be directed to the lead contact, Cengiz Pehlevan (cpehlevan@seas.harvard.edu).
Data and code availability
Mouse V1 neuron responses to orientation gratings and preprocessing code were obtained from a publicly available dataset: https://github.com/MouseLand/stringeretal2019, (Stringer et al., 2021; Pachitariu et al., 2019).
Responses to ImageNet images and preprocessing code were obtained from another publicly available dataset, https://github.com/MouseLand/stringerpachitariuetal2018b (Stringer et al., 2018a).
The code generated by the authors for this paper is also available https://github.com/PehlevanGroup/sample_efficient_pop_codes (PehlevanGroup, 2022).
Appendix 1
Singular value decomposition of continuous population responses
SVD of population responses is usually evaluated with respect to a discrete and finite set of stimuli. In the main paper, we implicitly assumed that a generalization of SVD to a continuum of stimuli. In this section we provide an explicit construction of this generalized SVD using techniques from functional analysis. Our construction is an example of the quasimatrix SVD defined in Townsend and Trefethen, 2015 and justifies our use of SVD in the main text.
For our construction, we note that Mercer’s theorem guarantees the existence of an eigendecomposition of any inner product kernel $K(\mathit{\theta},{\mathit{\theta}}^{\prime})$ in terms of a complete orthonormal set of functions ${\{{\psi}_{k}\}}_{k=1}^{\mathrm{\infty}}$ (Rasmussen and Williams, 2005). In particular, there exist a nonnegative (but possibly zero) summable eigenvalues ${\{{\lambda}_{k}\}}_{k=1}^{\mathrm{\infty}}$ and a corresponding set of orthonormal eigenfunctions such that
For a stimulus distribution $p(\mathit{\theta})$, the set of functions ${\{{\psi}_{k}\}}_{k=1}^{\mathrm{\infty}}$ are orthonormal and form a complete basis for square integrable functions L_{2} which means
Next, we use this basis to construct the SVD. Each of the tuning curves ${r}_{i}\in {L}_{2}$ (assumed to be square integrable) can be expressed in this basis with the top $N$ of the functions in the set
where we introduced a matrix $\mathbf{A}\in {\mathbb{R}}^{N\times N}$ of expansion coefficients. Note that $\text{rank}(\mathbf{A})\le N$. We compute the singular value decomposition of the finite matrix $\mathbf{A}$
We note that the signal correlation matrix for this population code can be computed in closed form
due to the orthonormality of $\{{\psi}_{k}\}$. Thus the principal axes ${\mathbf{u}}_{k}$ of the neural correlations are the left singular vectors of $\mathbf{A}$. We may similarly express the inner product kernel in terms of the eigenfunctions
The kernel eigenvalue problem demands (Rasmussen and Williams, 2005)
The ${\mathbf{v}}_{k}$ vectors must be identical to $\pm {\mathbf{e}}_{k}$, the Cartesian unit vectors, if the eigenvalues are nondegenerate. From this exercise, we find that the SVD for $\mathbf{A}$ has the form $\mathbf{A}=\sqrt{N}{\sum}_{k=1}^{\text{rank}(\mathbf{A})}\sqrt{{\lambda}_{k}}{\mathbf{u}}_{k}{\mathbf{e}}_{k}^{\top}$. With this choice, the population code admits a singular value decomposition
This singular value decomposition demonstrates the connection between neural manifold structure (principal axes ${\mathbf{u}}_{k}$) and function approximation (kernel eigenfunctions ${\psi}_{k}$). This singular value decomposition can be verified by computing the inner product kernel and the correlation matrix, utilizing the orthonormality of $\{{\mathbf{u}}_{k}\}$ and $\{{\psi}_{k}\}$. This exercise has important consequences for the space of learnable functions, which is at most $\text{rank}(\mathbf{A})$ dimensional since linear readouts lie in $\text{span}{\{{r}_{i}(\mathit{\theta})\}}_{i=1}^{N}$.
Discrete stimulus spaces: finding eigenfunctions with matrix eigendecomposition
In our discussion so far, our notation suggested that $\mathit{\theta}$ take a continuum of values. Here we want to point that our theory still applies if $\mathit{\theta}$ take a discrete set of values. In this case, we can think of a Dirac measure $p(\mathit{\theta})={\sum}_{i=1}^{\stackrel{~}{P}}{p}_{i}\delta (\mathit{\theta}{\mathit{\theta}}^{i})$, where $i$ indexes all the $\stackrel{~}{P}$ values $\mathit{\theta}$ can take. With this choice
Demanding this equality for $\mathit{\theta}}^{\mathrm{\prime}}={\mathit{\theta}}^{i},\text{}i=1,...,\stackrel{~}{P$ generates a matrix eigenvalue problem
where $\mathbf{B}}_{ij}={\delta}_{ij}{p}_{i$. The eigenfunctions over the stimuli are identified as the columns of $\mathbf{\Psi}$ while the eigenvalues are the diagonal elements of $\mathbf{\Lambda}}_{k\ell}={\lambda}_{k}{\delta}_{k\ell$.
Experimental considerations
In an experimental setting, a finite number of stimuli are presented and the SVD is calculated over this finite set regardless of the support of $p(\mathit{\theta})$. This raises the question of the interpretation of this SVD and its relation to the inductive bias theory we presented. Here we provide two interpretations.
In the first interpretation, we think of the empirical SVD as providing an estimate of the SVD over the full distribution $p(\mathit{\theta})$. To formalize this notion, we can introduce a MonteCarlo estimate of the integral eigenvalue problem
For this interpretation to work, the experimenter must sample the stimuli from $p(\mathit{\theta})$, which could be the natural stimulus distribution. Measuring responses to a larger number of stimuli gives a more accurate approximation of the integral above, which will provide a better estimate of generalization performance on the true distribution $p(\mathit{\theta})$.
In the second interpretation, we construct an empirical measure on $\stackrel{~}{P}$ experimental stimulus values $\widehat{p}(\mathit{\theta})=\frac{1}{\stackrel{~}{P}}{\sum}_{\mu =1}^{\stackrel{~}{P}}\delta (\mathit{\theta}{\mathit{\theta}}^{\mu})$, and consider learning and generalization over this distribution. This allows the application of our theory to an experimental setting where $\widehat{p}(\mathit{\theta})$ is designed by an experimenter. For example, the experimenter could procure a complicated set of $\stackrel{~}{P}$ videos, to which an associated function $y(\theta )$ must be learned. After showing these videos to the animal and measuring neural responses, the experimenter could compute, with our theory, generalization error for a uniform distribution over this full set of $\stackrel{~}{P}$ videos. Our theory would predict generalization over this distribution after providing supervisory feedback for only a strict subset of $P<\stackrel{~}{P}$ videos. Under this interpretation, the relationship between the integral eigenvalue problem and matrix eigenvalue problem is exact rather than approximate
Demanding either of (32) or (33) equalities for ${\mathit{\theta}}^{\prime}={\mathit{\theta}}^{\nu},\nu =1,\mathrm{\dots},P$ generates a matrix eigenvalue problem
The eigenfunctions restricted to $\{{\mathit{\theta}}^{\mu}\}$ are identified as the columns of $\mathbf{\Psi}$ while the eigenvalues are the diagonal elements of $\mathbf{\Lambda}}_{k\ell}={\lambda}_{k}{\delta}_{k\ell$. For the case where $N$ and $P$ are finite, the spectrum obtained through eigendecomposition of the kernel $\mathbf{K}$ is the same as would be obtained through the finite $N$ signal correlation matrix ${\mathbf{\Sigma}}_{s}$, since they are inner and outer products of trial averaged population response matrices $\mathbf{R}$.
Translation invariant kernels
For the special case where the data distribution $p(\mathit{\theta})=\frac{1}{V}$ is uniform over volume $V$ and the kernel is translation invariant $K(\mathit{\theta},{\mathit{\theta}}^{\prime})=\kappa (\mathit{\theta}{\mathit{\theta}}^{\prime})$, the kernel can be diagonalized in the basis of plane waves
The eigenvalues are the Fourier components of the Kernel ${\lambda}_{\mathbf{k}}=\frac{1}{V}\widehat{\kappa}(\mathbf{k})=\frac{1}{V}\int d\mathit{\theta}{e}^{i\mathbf{k}\cdot \mathit{\theta}}\kappa (\mathit{\theta})$ while the eigenfunctions are plane waves ${\psi}_{\mathbf{k}}(\mathit{\theta})={e}^{i\mathbf{k}\cdot \mathit{\theta}}$. The set of admissible momenta $\mathcal{S}}_{\mathbf{k}}=\{{\mathbf{k}}_{0},\pm {\mathbf{k}}_{1},\pm {\mathbf{k}}_{2},...\$ are determined by the boundary conditions. The diagonalized representation of the kernel is therefore
For example, if the space is the torus ${\mathbb{T}}^{n}={S}^{1}\times {S}^{1}\times \mathrm{\dots}\times {S}^{1}$, then the space of admissable momenta are the points on the integer lattice $\mathcal{S}}_{\mathbf{k}}={\mathbb{Z}}^{n}=\{\mathbf{k}\in {\mathbb{R}}^{n}{k}_{i}\in \mathbb{Z}\text{}\mathrm{\forall}i=1,...,n\$. Reality and symmetry of the kernel demand that $\text{Im}({\lambda}_{\mathbf{k}})=0$ and ${\lambda}_{\mathbf{k}}={\lambda}_{\mathbf{k}}\ge 0$. Most of the models in this paper consider $\theta \sim \text{Unif}\left({S}^{1}\right)$, where the kernel has the following Fourier/Mercer decomposition
where we invoked the simple trigonometric identity $\mathrm{cos}(ab)=\mathrm{cos}(a)\mathrm{cos}(b)+\mathrm{sin}(a)\mathrm{sin}(b)$. By recognizing that ${\{\sqrt{2}\mathrm{cos}(k\theta ),\sqrt{2}\mathrm{sin}(k\theta )\}}_{k=0}^{\mathrm{\infty}}$ form a complete orthonormal set of functions with respect to $\text{Unif}\left({S}^{1}\right)$, we have identified this as the collection of kernel eigenfunctions.
Invariant kernels possess invariant eigenfunctions
Suppose the kernel $K(\mathit{\theta},{\mathit{\theta}}^{\prime})$ is invariant to some set of transformations $t\in \mathcal{T}$, by which we mean that
We will now show that any eigenfunction of such a kernel with nonzero eigenvalue must be an invariant function. Let ${\psi}_{k}(\mathit{\theta})$ be an eigenfunction with eigenvalue ${\lambda}_{k}>0$, then
This establishes that all functions which depend on $\mathcal{T}$ transformations must necessarily lie in the nullspace of $K$.
Theory of generalization
Convergence of the deltarule without weight decay
In this section, we discuss the deltarule convergence when weight decay parameter is set to $\lambda =0$. The next section considers the simpler case where $\lambda >0$. Gradient descent training of readout weights $\mathbf{w}$ on a finite sample of size $P$ converges to the kernel regression solution (Bartlett et al., 2020; Hastie et al., 2020). Let $\mathcal{D}=\{{\mathit{\theta}}^{\mu},{y}^{\mu}{\}}_{\mu =1}^{P}$ be the dataset with samples ${\mathit{\theta}}^{\mu}$ and target values ${y}^{\mu}$. We introduce a shorthand ${\mathbf{r}}^{\mu}=\mathbf{r}({\mathit{\theta}}^{\mu})$ for convenience. The empirical loss we aim to minimize is a sum of the squared losses of each data point in the training set
Performing gradient descent updates
recovers the delta rule that we discussed in the main text (Widrow and Hoff, 1960; Hertz et al., 1991). Letting the empirical response matrix $\mathbf{R}=[{\mathbf{r}}^{1},\mathrm{\dots},{\mathbf{r}}^{P}]\in {\mathbb{R}}^{N\times P}$ have a SVD $\mathbf{R}=\sum _{k}\sqrt{{\hat{\lambda}}_{k}}{\hat{\mathbf{u}}}_{k}{\hat{\mathit{\psi}}}_{k}^{\mathrm{\top}}$, and expanding the weights ${\mathbf{w}}_{t}={\sum}_{k}{w}_{t,k}{\widehat{\mathbf{u}}}_{k}$ and labels $\mathbf{y}=\sum _{k}{\hat{v}}_{k}{\hat{\mathit{\psi}}}_{k}$ in their respective SVD bases, we find
For all directions with ${\hat{\lambda}}_{k}>0$, the dynamics converge to the unique fixed point ${w}_{k}^{*}=\frac{{\widehat{v}}_{k}}{\sqrt{{\widehat{\lambda}}_{k}}}$, while for all modes with ${\widehat{\lambda}}_{k}=0$, the weights remain at ${w}_{k}^{*}=0$. Thus
where K^{+} is the MoorePenrose inverse of the kernel matrix ${K}_{\mu ,\nu}=K({\mathit{\theta}}^{\mu},{\mathit{\theta}}^{\nu})$. The predictions of the learned function are given by $f={\mathbf{w}}^{*}\cdot \mathbf{r}(\mathit{\theta})$ which can be expressed as
The fact that the solution can be written in terms of a linear combination of ${\{K(\mathit{\theta},{\mathit{\theta}}^{\mu})\}}_{\mu =1}^{P}$ is known as the representer theorem (Schölkopf et al., 2001; Rasmussen and Williams, 2005). A similar analysis for nonlinear readouts where $f(\mathit{\theta})=g\left(\mathbf{w}\cdot \mathbf{r}(\mathit{\theta})\right)$ is provided in Appendix Convergence of deltarule for nonlinear readouts.
Weight decay and ridge regression
We can introduce a regularization term in our learning problem which penalizes the size of the readout weights. This leads to a modified learning objective of the form
Inclusion of this regularization alters the learning rule through weight decay${\mathbf{w}}_{t+1}=(1\eta \lambda ){\mathbf{w}}_{t}\eta \sum _{\mu}{\mathbf{r}}^{\mu}({\mathbf{r}}^{\mu}\cdot {\mathbf{w}}_{t}{y}^{\mu})$, which multiplies the existing weight value by a factor of $1\eta \lambda $ before adding the data dependent update. The fixed point of these dynamics is $\mathbf{w}={({\mathbf{RR}}^{\top}+\lambda \mathbf{I})}^{1}\mathbf{Ry}$. This learning problem and gradient descent dynamics have a closed form solution
The generalization benefits of explicit regularization through weight decay is known to be related to the noise statistics in the learning problem (Canatar et al., 2021). This is visible in the Appendix 1—figure 1 , where unlearnable target functions demand nonzero optimal regularization. We simulate weight decay only in Figure 6C, where we use $\lambda =0.01{\sum}_{k}{\lambda}_{k}$ to improve numerical stability at large $P$.
Computation of learning curves
Recent work has established analytic results that predict the average case generalization error for kernel regression
where $E}_{g}(\mathcal{D})={\u27e8(f(\mathit{\theta},\mathcal{D})y(\mathit{\theta}){)}^{2}\u27e9}_{\mathit{\theta}$ is the generalization error for a certain sample $\mathcal{D}$ of size $P$ and $f(\mathit{\theta},\mathcal{D})$ is the kernel regression solution for $\mathcal{D}$ (Bordelon et al., 2020; Canatar et al., 2021). The typical or average case error ${E}_{g}$ is obtained by averaging over all possible datasets of size $P$. This average case generalization error is determined solely by the decomposition of the target function $y(\mathbf{x})$ along the eigenbasis of the kernel and the eigenspectrum of the kernel. This diagonalization takes the form
Since the eigenfunctions form a complete set of square integrable functions, we expand both the target function $y(\mathit{\theta})$ and the learned function $f(\mathit{\theta})$ in this basis
Due to the orthonormality of the kernel eigenfunctions $\{{\psi}_{k}\}$, the generalization error for any set of coefficients $\widehat{v}$ is
We now introduce training error, or empirical loss, which depends on the disorder in the dataset $\mathcal{D}=\{({\mathit{\theta}}^{\mu},{y}^{\mu}){\}}_{\mu =1}^{P}$
It is straightforward to verify that the optimal ${\widehat{\mathbf{v}}}^{*}$ which minimizes $H(\hat{\mathbf{v}},\mathcal{D})$ is the kernel regression solution for kernel with eigenvalues $\{{\lambda}_{k}\}$ when $\lambda \to 0$. The optimal weights $\widehat{\mathbf{v}}$ can be identified through the first order condition $\mathrm{\nabla}H(\hat{\mathbf{v}},\mathcal{D})=0$ which gives
where ${\mathrm{\Psi}}_{k,\mu}={\psi}_{k}({\mathit{\theta}}^{\mu})$ are the eigenfunctions evaluated on the training data and ${\mathrm{\Lambda}}_{k,\mathrm{\ell}}={\delta}_{k,\mathrm{\ell}}{\lambda}_{k}$ is a a diagonal matrix containing the kernel eigenvalues. The generalization error for this optimal solution is
We note that the dependence on the randomly sampled dataset $\mathcal{D}$ only appears through the matrix $\mathbf{G}(\mathcal{D})$. Thus to compute the typical generalization error we need to average over this matrix $\u27e8\mathbf{G}(\mathcal{D})\u27e9}_{\mathcal{D}$. There are multiple strategies to perform such an average and we will study one here based on a partial differential equation which was introduced in Sollich, 1998; Sollich, 2002 and studied further in Bordelon et al., 2020; Canatar et al., 2021. In this setting, we denote the average matrix $\mathbf{G}(P)={\u27e8\mathbf{G}(\mathcal{D})\u27e9}_{\mathcal{D}=P}$ for a dataset of size $P$. We first will derive a recursion relationship using the Sherman Morrison formula for a rank1 update to an inverse matrix. We imagine adding a new sampled feature vector $\varphi $ to a dataset $\psi $ with size $P$. The average matrix $\mathbf{G}(P+1)$ at $P+1$ samples can be related to $\mathbf{G}(P)$ through the Sherman Morrison rule
where in the last step we approximated the average of the ratio with the ratio of averages. This operation, is of course, unjustified theoretically, but has been shown to produce accurate learning curves (Sollich, 2002; Bordelon et al., 2020). Since the chosen basis of kernel eigenfunctions are orthonormal, the average over the new sample is trivial $\u27e8\psi {\psi}^{\mathrm{\top}}\u27e9}_{\psi}=\mathbf{I$. We thus arrive at the following recursion relationship for $\mathbf{G}$
By introducing an additional source $J$ so that $\mathbf{G}(\mathcal{D},J{)}^{1}=\frac{1}{\lambda}\mathbf{\Psi}{\mathbf{\Psi}}^{\mathrm{\top}}+{\mathbf{\Lambda}}^{1}+J\mathbf{I}$, we can relate $\mathbf{G}(\mathcal{D},J)$’s first and second moments through differentiation
Thus the recursion relation simplifies to
where we approximated the finite difference in $P$ as a derivative, treating $P$ as a continuous variable. Taking the trace of both sides and defining $\kappa (P,J)=\lambda +\text{Tr}\mathbf{G}(P,J)$ we arrive at the following quasilinear PDE
with the initial condition $\kappa (0,J)=\lambda +\text{Tr}({\mathbf{\Lambda}}^{1}+J\mathbf{I}{)}^{1}$. Using the method of characteristics, we arrive at the solution $\kappa (P,J)=\lambda +\text{Tr}{\left({\mathbf{\Lambda}}^{1}+(v+\frac{P}{\kappa (P,J)})\mathbf{I}\right)}^{1}$. Using this solution to $\kappa $, we can identify the solution to $\mathbf{G}$
The generalization error, therefore can be written as
where $\gamma =P{\sum}_{k}\frac{{\lambda}_{k}^{2}}{{({\lambda}_{k}P+\kappa )}^{2}}$, giving the desired result. Note that $\kappa $ depends on $J$ implicitly, which is the source of the $\frac{1}{1\gamma}$ factor. This result was recently reproduced using techniques from statistical mechanics (Bordelon et al., 2020; Canatar et al., 2021).
Spectral bias and codetask alignment
Through implicit differentiation it is straightforward to verify that the ordering of the mode errors ${E}_{k}=\frac{{\kappa}^{2}}{1\gamma}{\left({\lambda}_{k}P+\kappa \right)}^{2}$ matches the ordering of the eigenvalues (Canatar et al., 2021). Let $\lambda}_{k}>{\lambda}_{\ell$, then we have
Since $\lambda}_{\ell}<{\lambda}_{k$, the first bracket must be negative and the second bracket must be positive. Further, it is straightforward to compute that ${\kappa}^{\mathrm{\prime}}(P)=\frac{\kappa \gamma}{P(1+\gamma )}<0$. Therefore $\lambda}_{k}>{\lambda}_{\ell$ implies $\frac{d}{dP}\mathrm{log}\left(\frac{{E}_{k}}{{E}_{\ell}}\right)<0$ for all $P$. Since $\mathrm{log}\left(\frac{{E}_{k}}{{E}_{\mathrm{\ell}}}\right)=0$ at $P=0$ we therefore have that $\mathrm{log}({E}_{k}/{E}_{\ell})<0$ for all $P$ and consequently $E}_{k}<{E}_{\ell$. Modes with larger eigenvalues ${\lambda}_{k}$ have lower normalized mode errors ${E}_{k}$. This observation can be used to prove that target functions acting on the same data distribution with higher cumulative power distributions $C(k)$ for all $k$ will have lower generalization error normalized by total target power, ${E}_{g}(P)/{E}_{g}(0)$, for all $P$. Proof can be found in Canatar et al., 2021.
Asymptotic power law scaling of learning curves
Exponential spectral decays:
First, we will study the setting relevant to the vonMises kernel where ${\lambda}_{k}\sim {\beta}^{k}$ and ${v}_{k}^{2}\sim {\alpha}^{k}$ where $\alpha ,\beta <1$. This exponential behavior accounts for differences in bandwidth between kernels which modulates the base $\beta $ of the exponential scaling of ${\lambda}_{k}$ with $k$. We will approximate the sum over all mode errors with an integral
If we include a regularization parameter $\lambda $, then $\kappa \sim \lambda $ as $P\to \mathrm{\infty}$. Making the change of variables $u=P{\beta}^{k}/\lambda$, we transform the above integral into
The remaining integral over $u$ is either dominated near $u\approx 0$ if $\frac{\mathrm{ln}(1/\alpha )}{\mathrm{ln}(1/\beta )}<2$ and behaves as a $P$independent constant or else is dominated near $u\approx P/\lambda$, in which case the integral scales as $\sim {P}^{\frac{\mathrm{ln}(1/\alpha )}{\mathrm{ln}(1/\beta )}2}$. Multiplying these resulting functions with the prefactor, we find the following scaling laws for generalization.
Thus, we obtain a power law scaling of the learning curve ${E}_{g}$ which is dominated at large $P$ by ${E}_{g}\sim {P}^{\mathrm{min}(2,\frac{\mathrm{ln}(1/\alpha )}{\mathrm{ln}(1/\beta )})}$. For the vonMises kernel we can approximate the spectra with ${\lambda}_{k}\sim {\sigma}^{2k}$ and ${v}_{k}^{2}\sim {\sigma}_{T}^{2k}$ giving rise to a generalization scaling scaling ${E}_{g}\sim {P}^{\mathrm{min}(2,\frac{\mathrm{ln}{\sigma}_{T}}{\mathrm{ln}\sigma})}$.
Power law spectral decays
The same arguments can be applied for power law kernels ${\lambda}_{k}\sim {k}^{b}$ and power law targets ${v}_{k}^{2}\sim {k}^{a}$, which is of interest due to its connection to nonlinear rectified neural populations. In this setting, the generalization error is
We see that there are two possible power law scalings for ${E}_{g}$ with the exponents $(a1)/b$ and 2. At large $P$ this formula will be dominated by the term with minimum exponent so ${E}_{g}\sim {P}^{\mathrm{min}(a1,2b)/b}$.
Laplace kernel generalization
We calculate similar learning curves as we did for the vonMises kernel but with Laplace kernels to show that our results is not an artifact of the infinite differentiability of the Von Mises kernel. Each of these Laplace kernels has the same asymptotic power law spectrum ${\lambda}_{k}\sim o({k}^{2})$, exhibiting a discontinuous first derivative (Figure 6A). Despite having the same spectral scaling at large $k$, these kernels can give dramatically different performance in learning tasks, again indicating the influence of the top eigenvalues on generalization at small $P$ (Figure 6). Again, the trend for which kernels perform best at low $P$ can be reversed at large $P$. In this case, all generalization errors scale with ${E}_{g}\sim {P}^{2}$ (Figure 6B). More generally, our theory shows that if the task power spectrum and kernel eigenspectrum are both falling as power laws with exponents $a$ and $b$ respectively, then the generalization error asymptotically falls with a power law, ${E}_{g}\sim {P}^{\mathrm{min}(a1,2b)/b}$ (Methods) (Bordelon et al., 2020). This decay is fastest when $b\ge \frac{a1}{2}$ for which ${E}_{g}\sim {P}^{2}$. Therefore, the tail of the kernel’s eigenvalue spectrum determines the large sample size behavior of the generalization error for power law kernels. Small sample size limit is still governed by the bulk of the spectrum.
Learning with multiple output channels
Our theory is not limited to scalar target functions but rather can be easily extended to multiple output functions ${y}_{1},\mathrm{\dots},{y}_{C}$ from the same data, if for example the task requires computing class membership for $C$ categories. In this setting, each data point has the form $({\mathit{\theta}}^{\mu},{\mathbf{y}}^{\mu})$ with ${\mathbf{y}}^{\mu}\in {\mathbb{R}}^{C}$. For these $C$ classes, the generalization error takes the form
We therefore find that the generalization error in the multiclass setting is the same as the ${E}_{g}$ obtained for a single scalar target function with power spectrum $v}_{k}^{2}=\sum _{c}{\u27e8{y}_{c}(\mathit{\theta}){\varphi}_{k}(\mathit{\theta})\u27e9}^{2$(Bordelon et al., 2020; Canatar et al., 2021). The relevant cumulative power distribution measures the fraction of total output variance captured by the first $k$ eigenfunctions of the population code
Convergence of Deltarule for nonlinear readouts
In this section, we consider gradient descent dynamics of an arbitrary convex loss function. For instance, we can consider a binary classification problem where $y\in \{\pm 1\}$ by outputting a prediction of $\widehat{y}=\text{sign}(\mathbf{w}\cdot \mathbf{r})$. We could, for example, train a model using the hinge loss $\mathrm{\ell}(\mathbf{w}\cdot \mathbf{r},y)=\text{max}(0,1\mathbf{w}\cdot \mathbf{r}y)$ so that the classifier will converge to a kernel support vector machine (SVM) (Schölkopf et al., 2002). The generalization of the classifier would be the error rate of $\widehat{y}(\mathit{\theta})=\text{sign}(\mathbf{w}\cdot \mathbf{r}(\mathit{\theta}))$ compared to the ground truth $y(\mathit{\theta})$.
Let $\mathcal{D}=\{{\mathit{\theta}}^{\mu},{y}^{\mu}{\}}_{\mu =1}^{P}$ be the dataset with samples ${\mathit{\theta}}^{\mu}$ and target values ${y}^{\mu}$. We introduce a shorthand ${\mathbf{r}}^{\mu}=\mathbf{r}({\mathit{\theta}}^{\mu})$ for convenience. The loss we aim to minimize is the sum of the losses of each data point in the training set with an additional weight decay parameter
For convex $\mathrm{\ell}$ and nonzero $\lambda $, the above objective is strongly convex, indicating the existence of a unique minimizer which can be found from simple first order learning rules. For $\lambda >0$ the initial condition for $\mathbf{w}$ does not influence the final result.
We will now show that the dynamics will converge to a function which only depends on the code $\mathbf{r}(\mathit{\theta})$ through the kernel $K(\mathit{\theta},{\mathit{\theta}}^{\prime})$. To simplify the argument, we consider starting from an initial condition of ${\mathbf{w}}_{t=0}=\mathbf{0}$ and performing gradient descent updates. Under such an assumption, the weights ${\mathbf{w}}_{t}$ will always be in the span of the population vectors on the training set $\{{\mathbf{r}}^{\mu}{\}}_{\mu =1}^{P}$ since
The derivative in the final term is taken with respect to the first argument ${\mathrm{\ell}}^{\prime}(f,y)=\frac{\partial \mathrm{\ell}(f,y)}{\partial f}$. This update is still local and recovers the delta rule that we discussed in the main text for $\mathrm{\ell}(\mathbf{w}\cdot \mathbf{r},y)=\frac{1}{2}{(\mathbf{w}\cdot \mathbf{r}y)}^{2}$(Widrow and Hoff, 1960; Hertz et al., 1991). We can express ${\mathbf{w}}_{t}$ in terms of the population vectors $\mathbf{w}}_{t}=\sum _{\mu =1}^{P}{\alpha}_{t}^{\mu}\text{}{\mathbf{r}}^{\mu}=\mathbf{R}{\mathit{\alpha}}_{t$ so that ${\mathit{\alpha}}_{t}\in {\mathbb{R}}^{P}$ defines the linear weighting of each sample. The dynamics of these coefficients are
where $\mathbf{K}={\mathbf{R}}^{\top}\mathbf{R}\in {\mathbb{R}}^{P\times P}$ is the kernel Gram matrix evaluated on the training points. We multiply both sides of this equation by ${\mathbf{R}}^{\top}$, and define ${\mathit{\beta}}_{t}=\mathbf{K}{\mathit{\alpha}}_{t}$, which satisfy the following simplified dynamics
where K^{+} is the pseudoinverse of $\mathbf{K}$. The nonlinear fixed point condition is $\mathit{\beta}=\frac{1}{\lambda}\mathbf{K}{\mathrm{\ell}}^{\prime}(\mathit{\beta},\mathbf{y})$, which transparently only depends on the kernel $\mathbf{K}$ rather than the full code $\mathbf{R}$. The above equation recovers the correct linear equation $\mathit{\beta}=\mathbf{K}{\left(\mathbf{K}+\lambda \mathbf{I}\right)}^{1}\mathbf{y}$ for square loss. For an arbitrary test point $\mathit{\theta}$, the model makes prediction using $f(\mathit{\theta})=\mathbf{r}(\mathit{\theta})\cdot \mathbf{w}=\mathbf{r}(\mathit{\theta})\cdot [{\mathbf{RK}}^{+}\mathit{\beta}]=\mathbf{k}(\mathit{\theta})\cdot \mathit{\alpha}$, which also only depends on the kernel $[\mathbf{k}(\mathit{\theta}){]}_{\mu}=K(\mathit{\theta},{\mathit{\theta}}^{\mu})$ on test point $\mathit{\theta}$ and train points ${\mathit{\theta}}_{\mu}$, as well as the kernel gram matrix $[\mathbf{K}{]}_{\mu \nu}=K({\mathit{\theta}}^{\mu},{\mathit{\theta}}^{\nu})$.
To illustrate a specific case with a square error and nonlinear readout, consider output neurons which produce activity $g(\mathbf{w}\cdot \mathbf{r}(\mathit{\theta}))$ for invertible nonlinear function $g$ with nonvanishing gradient, and gradient based learning on $L={\sum}_{\mu}{(g(\mathbf{w}\cdot {\mathbf{r}}^{\mu})y({\mathit{\theta}}^{\mu}))}^{2}$. This gives $\mathrm{\Delta}\mathbf{w}\propto {\sum}_{\mu}{\mathbf{r}}^{\mu}{g}^{\prime}(\mathbf{w}\cdot {\mathbf{r}}^{\mu})({y}^{\mu}g(\mathbf{w}\cdot {\mathbf{r}}^{\mu}))\in \text{span}{\{{\mathbf{r}}^{\mu}\}}_{\mu =1}^{P}$, which is still a local learning rule. Thus the weights at convergence can be written as $\mathbf{w}=\sum _{\mu}{\alpha}^{\mu}{\mathbf{r}}^{\mu}$ and the learned function can be written as $f(\mathit{\theta})=g\left(\sum _{\mu}{\alpha}^{\mu}K(\mathit{\theta},{\mathit{\theta}}^{\mu})\right)$ is given by $\mathit{\alpha}={\mathbf{K}}^{+}{g}^{1}(\mathbf{y})$. To see this, first note that ${\mathbf{w}}_{t}\in \text{span}{\{{\mathbf{r}}^{\mu}\}}_{\mu =1}^{P}$ for all $t$ so that ${\mathbf{w}}^{*}=\mathbf{R}{\mathit{\alpha}}^{*}$. The interpolation condition can be expressed as $g\left({\mathbf{R}}^{\top}{\mathbf{w}}^{*}\right)=g\left(\mathbf{K}\mathit{\alpha}\right)=\mathbf{y}$, giving the desired result ${\mathit{\alpha}}^{*}={\mathbf{K}}^{+}{g}^{1}(\mathbf{y})$. The predictions of the model on a test stimulus $\mathit{\theta}$ are given by $f(\mathit{\theta})=g\left(\sum _{\mu =1}^{P}{\alpha}^{\ast \mu}K(\mathit{\theta},{\mathit{\theta}}^{\mu})\right)$. We see that this solution only depends on the kernel (directly and indirectly through ${\alpha}^{*}$), rather than the full code.
Typical case analysis of nonlinear readouts
The analysis of typical case generalization can be extended to nonlinear predictors and loss functions described by (68) which depend on the scalar prediction variable $\mathbf{w}\cdot \mathbf{r}(\mathit{\theta})$ (Loureiro et al., 2021a). Thanks Further, if $\mathbf{r}(\mathit{\theta})$ is well approximated as a Gaussian process, then the generalization performance can still be characterized using statistical mechanics methods (Loureiro et al., 2021a). Many qualitative features of our results continue to hold, including that the kernel’s diagonalization governs training and generalization and that improvements in code task alignment lead to improvements in generalization.
In a later work by Cui et al., 2022, SVM and ridge classifiers trained on codes and tasks with power law spectra were analyzed asymptotically, showing power law generalization error decay rates ${E}_{g}\sim {P}^{\beta}$. These classification learning curves for power law spectra were shown to follow power laws with exponents $\beta $ which are qualitatively similar to the exponents obtained with the square loss which we describe in our section titled Small and Large Sample Size Behaviors of Generalization. Just as in our theory, decay rate exponents $\beta $ are larger for codes which are well aligned to the task and are smaller for codes which are nonaligned.
Visual scene reconstruction task
Reconstruction of natural scenes from neural responses
Using the mouse V1 responses to natural scenes, we attempt to reconstruct original images from the neural codes using different numbers of images. The presented natural scenes are taken from ten classes of imagenet which can be downloaded from https://github.com/MouseLand/stringerpachitariuetal2018b. Let ${\mathit{\theta}}^{\mu}\in {\mathbb{R}}^{D}$ be a $D$dimensional flattened vector containing the pixel values of the μth image and let ${\mathbf{r}}^{\mu}\in {\mathbb{R}}^{N}$ represent the neural response to the μth image. The goal in the problem is to learn a collection of weights $\mathbf{W}\in {\mathbb{R}}^{D\times N}$ which map neural responses ${\mathbf{r}}^{\mu}$ to images ${\mathit{\theta}}^{\mu}$
The generalization error ${E}_{g}$ again measures the average error on all points, averaged over all possible datasets $\mathcal{D}=\{({\mathit{\theta}}^{\mu},{\mathbf{r}}^{\mu}){\}}_{\mu =1}^{P}$ of size $P$. If the optimal weights for dataset $\mathcal{D}$ is $\mathbf{W}(\mathcal{D})$ then the generalization error is
After identifying eigenfunctions ${\varphi}_{k}(\mathit{\theta})$, we expand the images in this basis $\mathit{\theta}=\sum _{k}{\mathbf{v}}_{k}{\psi}_{k}(\mathit{\theta})$ where ${\mathbf{v}}_{k}\in {\mathbb{R}}^{D}$. The generalization error is therefore ${E}_{g}={\sum}_{k}{{\mathbf{v}}_{k}}^{2}{E}_{k}(P)$ and the cumulative power is $C(k)=\frac{\sum _{\ell <k}{\mathbf{v}}_{\ell}{}^{2}}{\sum _{\ell =1}^{\mathrm{\infty}}{\mathbf{v}}_{\ell}{}^{2}}$. We perform this reconstruction task on many filtered versions of the natural scenes. To construct a filter, we first compute the Fourier transform of the image. Let $\mathbf{M}(\mathit{\theta})\in {\mathbb{R}}^{\sqrt{D}\times \sqrt{D}}$ represent the nonflattened image and let $\widehat{\mathbf{M}}(\mathit{\theta})\in {\mathbb{R}}^{\sqrt{D}\times \sqrt{D}}$ represent the Fourier transform of the image, computed explicitly as
To develop the bandpass filter, we calculate $\mathbf{k}=\sqrt{{k}^{2}+{\left({k}^{\prime}\right)}^{2}}$ for each of the indices in the matrix. For a bandpass filter with parameters ${s}_{max},r$ we simply zero out the entries in $\widehat{M}$ which correspond to states with frequencies outside the appropriate band: for any $k,{k}^{\prime}$ with $\mathbf{k}\notin [\sqrt{{s}_{max}^{2}{r}^{2}},{s}_{max}^{2}]$ then ${\widehat{M}}_{k,{k}^{\prime}}\to 0$. We then perform the inverse Fourier transform on $\widehat{M}$ to obtain a filtered version of the original image.
A simple feedforward model of V1
Linear neurons
We consider a simplified but instructive model of the V1 population code as a linearnonlinear map from photoreceptor responses through Gabor filters and then nonlinearity (Adelson and Bergen, 1985; Olshausen and Field, 1997; Rumyantsev et al., 2020). Let $\mathbf{x}\in {\mathbb{R}}^{2}$ represent the twodimensional retinotopic position of photoreceptors. The firing rates of the photoreceptor at position $\mathbf{x}$ to a static grating stimulus oriented at angle $\theta $ is
We model each V1 neuron’s receptive field as a Gabor filter of the receptor responses $h(\mathbf{x},\theta )$. The $i$th V1 neuron has preferred wavevector ${\mathbf{k}}_{i}$, generating the following set of weights between photoreceptors and the $i$th V1 neuron
The V1 population code is obtained by filtering the photoreceptor responses. By approximating the resulting sum over all retinal photoreceptors with an integral, we find the response of neuron $i$ to grating stimulus with wavenumber $\mathbf{k}$ is
The response of neuron $i$ is computed through nonlinear rectification of this input current ${r}_{i}(\theta )=g(\mathbf{w}({\theta}_{i})\cdot h(\theta ))$. For a linear neuron $g(z)=z$, the kernel has the following form
where the kernel is normalized to have maximum value of 1. Note that this normalization of the kernel is completely legitimate since it merely rescales each eigenvalue by a constant and does not change the learning curves.
Since the kernel only depends on the difference between angles $\theta {\theta}^{\prime}$, it is said to posess translation invariance. Such translation invariant kernels admit a Mercer decomposition in terms of Fourier modes $K(\theta )={\sum}_{n}{\lambda}_{n}\mathrm{cos}(n\theta )$ since the Fourier modes diagonalize shift invariant integral operators on ${S}^{1}$. For the linear neuron, the kernel eigenvalues scale like ${\lambda}_{n}\sim \frac{{\beta}^{n}}{{2}^{n}n!}$, indicating infinite differentiability of the tuning curves. Since ${\lambda}_{n}$ decays rapidly with $n$, we find that this Gabor code has an inductive bias that favors low frequency functions of orientation $\theta $.
Nonlinear simple cells
Introducing nonlinear functions $g(z)$ that map input currents $z$ into the V1 population into firing rates, we can obtain a nonlinear kernel ${K}_{g}(\theta )$ which has the following definition
In this setting, it is convenient to restrict ${\mathbf{k}}_{i},\mathbf{k},{\mathbf{k}}^{\prime}\in {S}^{1}$ and assume that the preferred wavevectors ${\mathbf{k}}_{i}$ are uniformly distributed over the circle. In this case, it suffices to identify a decomposition of the composed function $g({\mathbf{w}}_{i}\cdot \mathbf{h}(\theta ))$ in the basis of Chebyshev polynomials ${T}_{n}(z)$ which satisfy ${T}_{n}(\mathrm{cos}(\theta ))=\mathrm{cos}(n\theta )$
which can be computed efficiently with an appropriate quadrature scheme. Once the coefficients a_{n} are determined, we can compute the kernel by first letting ${\theta}_{i}$ to be the angle between $\mathbf{k}$ and ${\mathbf{k}}_{i}$ and letting $\theta $ be the angle between $\mathbf{k}$ and ${\mathbf{k}}^{\prime}$
Thus the kernel eigenvalues are ${\lambda}_{n}=\frac{1}{2}{a}_{n}^{2}(\psi )$.
Asymptotic scaling of spectra
Activation functions that encourage sparsity have slower eigenvalue decays. If the nonlinear activation function has the form ${g}_{q,a}(z)=\mathrm{max}{\{0,za\}}^{q}$, then the spectrum decays like ${\lambda}_{n}\sim {n}^{2q2}$. A simple argument justifies this scaling: if the function $g({e}^{{\sigma}^{2}}\mathrm{cosh}({\sigma}^{2}z))$ is only $q1$ times differentiable then ${a}_{n}{n}^{q}\sim {n}^{1}$ since ${\sum}_{n}{a}_{n}{n}^{q}$ must diverge. Therefore ${\lambda}_{n}={a}_{n}^{2}\sim {n}^{2q2}$. Note that this scaling is independent of the threshold. Examples of these scalings can be found in Figure 5—figure supplements 1 and 2.
Phase variation, complex cells and invariance
We can consider a slightly more complicated model where Gabors and stimuli have phase shifts
The simple cells are generated by nonlinearity
The input currents into the simple V1 cells can be computed exactly
When $\mathbf{k}={\mathbf{k}}_{i}=1$, the simple cell tuning curves ${r}_{i}=g({\mathbf{w}}_{i}\cdot \mathbf{h})$ only depend on $\mathrm{cos}(\theta {\theta}_{i})$ and $\varphi $, allowing a Fourier decomposition
The simple cell kernel ${K}_{s}$, therefore decomposes into Fourier modes over $\theta $
where $b}_{n}(\varphi ,{\varphi}^{\mathrm{\prime}})={\u27e8{a}_{n}(\varphi ,{\varphi}_{i}){a}_{n}({\varphi}^{\mathrm{\prime}},{\varphi}_{i})\u27e9}_{{\varphi}_{i}$. It therefore suffices to solve the infinite sequence of integral eigenvalue problems over $\varphi $
With this choice it is straightforward to verify that the kernel eigenfunctions are ${v}_{n,k}(\theta ,\varphi )={e}^{in\theta}{v}_{n,k}(\varphi )$ with corresponding eigenvalue ${\lambda}_{n,k}$. Since b_{n} is not translation invariant in $\varphi {\varphi}^{\prime}$, the eigenfunctions ${v}_{n,k}$ are not necessarily Fourier modes. These eigenvalue problems for b_{n} must be solved numerically when using arbitrary nonlinearity $g$. The top eigenfunctions of the simple cell kernel depend heavily on the phase of the two grating stimuli $\varphi $. Thus, a pure orientation discrimination task which is independent of phase requires a large number of samples to learn with the simple cell population.
Complex cell populations are phase invariant
V1 also contains complex cells which possess invariance to the phase $\varphi $ of the stimulus.
Again using Gabor filters we model the complex cell responses with a quadratic nonlinearity and sum over two squared filters which are phase shifted by $\pi /2$
which we see is independent of the phase $\varphi $ of the grating stimulus. Integrating over the set of possible Gabor filters $({\mathbf{k}}_{i},{\varphi}_{i})$ with $\mathbf{k}=1$ again gives the following kernel for the complex cells
Remarkably, this kernel is independent of the phase $\varphi $ of the grating stimulus. Thus, complex cell populations possess good inductive bias for vision tasks where the target function only depends on the orientation of the stimulus rather than it’s phase. In reality, V1 is a mixture of simple and complex cells. Let $s\in [0,1]$ represent the relative proportion of neurons which are simple cells and $(1s)$ the relative proportion of complex cells. The kernel for the mixed V1 population is given by a simple convex combination of the simple and complex cell kernels
where $n$ denotes neuron type (simple vs complex, tuning etc) and ${p}_{V1}(n),{p}_{s}(n),{p}_{c}(n)$ are probability distributions over the V1 neuron identities, the simple cell identities and the complex cell identities respectively. Increasing $s$ increases the phase dependence of the code by giving greater weight to the simple cell population. Decreasing $s$ gives weight to the complex cell population, encouraging phase invariance of readouts.
Visualization of feedforward Gabor V1 model and induced kernels
Examples of the induced kernels for the Gaborbank V1 model are provided in Figure 5. We show how choice of rectifying nonlinearity $g(z)$ and sparsifying threshold $a$ influence the kernel and their spectra. Learning curves for simple orientation tasks are provided.
Gabor model spectral bias and fit to V1 data
Motivated by findings in the primary visual cortex (Hansel and van Vreeswijk, 2002; Miller and Troyer, 2002; Priebe et al., 2004; Priebe and Ferster, 2008), we studied the spectral bias induced by rectified powerlaw nonlinearities of the form $g(z)=\mathrm{max}{\{0,za\}}^{q}$. From theory, such a powerlaw activation function arises in a spiking neuron when firing is driven by input fluctuations (Hansel and van Vreeswijk, 2002; Miller and Troyer, 2002). Further, this activation is observed in intracellular recordings over the dynamic range of neurons in primary visual cortex (Priebe and Ferster, 2008). For example, in cats, the power, $q$, ranges from 2.7 to 3.9 (Priebe et al., 2004). We fit parameters of our model to the Mouse V1 kernel and compared to other parameter sets in Figure 5—figure supplement 1. Our best fit value of $q=1.7$ is lower but on par with the estimates from the cat and reproduces the observed kernel. Computation of the kernel and its eigenvalues (Appendix Nonlinear simple cells) indicates a low frequency bias: the eigenvalues for low frequency modes are higher than those for high frequency modes, indicating a strong inductive bias to learn functions of low frequency in the orientation. Decreasing sparsity (lower $a$) leads to a faster decay in the spectrum (but similar asymptotic scaling at the tail, see Figure 5—figure supplements 1 and 2) and a stronger bias towards lower frequency functions (Figure 5). The effect of the power of nonlinearity $q$ is more nuanced: increasing power may increase spectra at lower frequencies, but may also lead to a faster decay at the tail (Figure 5—figure supplements 1 and 2 ). In general, an exponent $q$ implies a powerlaw asymptotic spectral decay ${\lambda}_{k}\sim {k}^{2q2}$ as $k\to \mathrm{\infty}$ (Appendix Nonlinear simple cells). The behavior at low frequencies may have significant impact for learning with few samples. Overall, our findings show that the spectral bias of a population code can be determined in nontrivial ways by its biophysical parameters, including neural thresholds and nonlinearities.
Energy model with partially phaseselective cells
The model of the V1 population as a mixture of purely simple and purely complex cells is an idealization which fails to capture the variability in phase selectivity of cells observed in experiment. In this section, we describe a simple model which can interpolate between an invariant code and a code which has high alignment with phasedependent eigenfunctions. Further, a single scalar parameter $\alpha $ will control how strongly the population is biased towards invariance. We define ${r}_{i}(\theta ,\varphi )=g\left({z}_{i}(\theta ,\varphi )\right)$ for nonlinear function $g$ and scalar $z$ which is constructed as follows
This linear combination is inspired by the construction of simple cells in Dayan & Abbot Chapter 2 (Dayan and Abbott, 2001). If all $\beta $ are equal, then this tuning curve is invariant to phase $\varphi $. To generate variability in selectivity to phase $\varphi $, we will draw $\mathit{\beta}$ from a Dirichlet distrbution on the simplex with concentration parameter $\alpha \mathbf{1}$ so that $p(\mathit{\beta})\propto {\prod}_{j=1}^{4}{\beta}_{j}^{\alpha 1}$ with ${\sum}_{j=1}^{4}{\beta}_{j}=1$. In the $\alpha \to \mathrm{\infty}$ limit, the probability density concentrates on $\frac{1}{4}\mathbf{1}$, leading to a code comprised entirely of complex cells which are invariant to phase $\varphi $. In the $\alpha \to 0$ limit, the density is concentrated around the “edges” of the simplex such as $(1,0,0,0),(0,1,0,0)$, where only one preferred phase is present per neuron. For intermediate values, neurons are partially selective to phase. As before, the selectivity or invariance to phase is manifested in the kernel decomposition and leads to similar learning curves for the three tasks of the main paper (Orientation, Phase, Hybrid). We provide an illustration of tuning curves, F1/F0 distributions, eigenfunctions, and learning curves in Figure 5—figure supplement 3.
Time dependent neural codes
RNN model and decomposition
In this setting, the population code $\mathbf{r}(\{\mathit{\theta}(t)\},t)$ is a function of an input stimulus sequence $\mathit{\theta}(t)$ and time $t$. In general the neural code $\mathbf{r}$ at time $t$ can depend on the entire history of the stimulus input $\mathit{\theta}({t}^{\prime})$ for ${t}^{\prime}\le t$, as is the case for recurrent neural networks. We denote dependence of a function $f$ on $\mathit{\theta}(t)$ in this causal manner with the notation $f(\{\mathit{\theta}\},t)$. In a learning task, a set of readout weights $\mathbf{w}$ are chosen so that a downstream linear readout $f(\{\mathit{\theta}\},t)=\mathbf{w}\cdot \mathbf{r}(\{\mathit{\theta}\},t)$ approximates a target sequence $y(\{\mathit{\theta}\},t)$ which maps input stimulus sequences to output scalar sequences. The quantity of interest is the generalization ${E}_{g}$, which in this case is an average over both input sequences and time, ${E}_{g}={\u27e8{(y(\{\mathit{\theta}\},t)f(\{\mathit{\theta}\},t))}^{2}\u27e9}_{\mathit{\theta}(t),t}$. The average is computed over a distribution of input stimulus sequences $p(\mathit{\theta}(t))$. To train the readout, $\mathbf{w}$, the network is given a sample of $P$ stimulus sequences ${\mathit{\theta}}^{\mu}(t),\mu =1,\mathrm{\dots},P$. For the μth training input sequence, the target system $y$ is evaluated at a set of discrete time points $\mathcal{T}}_{\mu}=\{{t}_{1},{t}_{2},...,{t}_{{\mathcal{T}}_{\mu}}\$ giving a collection of target values $\{{y}_{t}^{\mu}{\}}_{t\in {\mathcal{T}}_{\mu}}$ and a total dataset of size $\mathcal{P}=\sum _{\mu =1}^{P}{\mathcal{T}}_{\mu}$. The average case generalization computes a further average of the generalization error ${E}_{g}$ over randomly sampled datasets of size $\mathcal{P}$.
Learning is again achieved through iterated weight updates with deltarule form, but now have contributions from both sequence index and time $\mathrm{\Delta}\mathbf{w}=\eta \sum _{\mu}\sum _{t\in {\mathcal{T}}_{\mu}}{\mathbf{r}}_{t}^{\mu}({y}_{t}^{\mu}{f}_{t}^{\mu})$. As before, optimization of the readout weights is equivalent to kernel regression with a kernel that computes inner products of neural population vectors at different times $t,{t}^{\prime}$ for different input sequences $\{\mathit{\theta}\},\{{\mathit{\theta}}^{\mathrm{\prime}}\}$ . This kernel depends on details of the time varying population code including its recurrent intrinsic dynamics as well as its encoding of the timevarying input stimuli. The optimization problem and delta rule described above converge to the kernel regression solution for kernel gram matrix ${K}_{t,{t}^{\prime}}^{\mu ,{\mu}^{\prime}}=\frac{1}{N}{\mathbf{r}}_{t}^{\mu}\cdot {\mathbf{r}}_{{t}^{\prime}}^{{\mu}^{\prime}}$(Dong et al., 2020; Yang, 2019; Yang, 2020). The learned function has the form $f(\{\mathit{\theta}\},t)=\sum _{\mu ,{t}^{\mathrm{\prime}}\in {\mathcal{T}}_{\mu}}{\alpha}_{t}^{\mu}K(\{\mathit{\theta}\},\{\mathit{\theta}{\}}^{\mu},t,{t}^{\mathrm{\prime}})$, where $\mathit{\alpha}={\mathbf{K}}^{+}\mathbf{y}$ for kernel gram matrix $\mathbf{K}\in {\mathbb{R}}^{\mathcal{P}\times \mathcal{P}}$ which is computed for the entire set of training sequences, and the vector $\mathbf{y}\in {\mathbb{R}}^{\mathcal{P}}$ is the vector containing the desired target outputs for each sequence. Assuming a probability distribution over sequences $\mathit{\theta}(t)$, the kernel can be diagonalized with orthonormal eigenfunctions ${\psi}_{k}(\{\mathit{\theta}\},t)$. Our theory carries over from the static case: kernels whose top eigenfunctions have high alignment with the target dynamical system $y(\{\mathit{\theta}\},t)$ will achieve the best average case generalization performance.
Alternative neural codes with same kernel
Orthogonal transformations are sufficient for linear kernelpreserving transformations
We will now show that for any linear transformation $\stackrel{~}{\mathbf{r}}=\mathbf{Ar}$ which preserves the inner product kernel $K(\mathit{\theta},{\mathit{\theta}}^{\prime})$, there exists an orthogonal matrix $\mathbf{Q}$ such that $\stackrel{~}{\mathbf{r}}=\mathbf{Qr}$.
Proof.
Let $\stackrel{~}{\mathbf{r}}(\mathit{\theta})=\mathbf{Ar}(\mathit{\theta})$ for all stimuli $\mathit{\theta}$. To preserve the kernel, we must have
Taking projections against each of the orthonormal eigenfunctions ${\psi}_{\mathrm{\ell}}(\theta )$ (see Appendix Singular value decomposition of continuous population responses), we define vectors ${\mathbf{u}}_{k}$ as $\sqrt{{\lambda}_{k}}{\mathbf{u}}_{k}={\u27e8\mathbf{r}(\mathit{\theta}){\psi}_{k}(\mathit{\theta})\u27e9}_{\mathit{\theta}}$, allowing us to express the SVD of the population code $\mathbf{r}(\mathit{\theta})={\sum}_{k}\sqrt{{\lambda}_{k}}{\mathbf{u}}_{k}{\psi}_{k}(\mathit{\theta})$. These vectors $\{{\mathbf{u}}_{k}\}$ are orthonormal ${\mathbf{u}}_{k}\cdot {\mathbf{u}}_{\mathrm{\ell}}={\delta}_{k\mathrm{\ell}}$ since, by the definition of the kernel eigenfunctions ${\psi}_{k}$,
Since $\mathbf{r}(\mathit{\theta})$ and $\stackrel{~}{\mathbf{r}}(\mathit{\theta})$ have the same inner product kernel, they must posess the same kernel eigenfunctions ${\psi}_{k}$ and kernel eigenvalues ${\lambda}_{k}$, which are identified through the eigenvalue problem
We therefore have the following two singular value decompositions for $\mathbf{r}$ and $\stackrel{~}{\mathbf{r}}$
where ${\{{\mathbf{u}}_{k}\}}_{k=1}^{N}$ and ${\{{\stackrel{~}{\mathbf{u}}}_{k}\}}_{k=1}^{N}$ are both complete sets of orthonormal vectors (the sums above run over possible zero eigenvalues). Taking the equation $\stackrel{~}{\mathbf{r}}(\mathit{\theta})=\mathbf{Ar}(\mathit{\theta})$, we multiply both sides of the equation by ${\psi}_{k}(\mathit{\theta})$ and average over $\mathit{\theta}$ giving
For an eigenmode $k$ with positive eigenvalue ${\lambda}_{k}>0$, this implies ${\stackrel{~}{\mathbf{u}}}_{k}={\mathbf{Au}}_{k}$, while there is no corresponding constraint for the null modes with ${\lambda}_{k}=0$. However, the action of $\mathbf{A}$ on the nullspace of the code has no influence on $\stackrel{~}{\mathbf{r}}$ so there is no loss in generality to restrict consideration to transformations $\mathbf{A}$ which satisfy ${\stackrel{~}{\mathbf{u}}}_{k}={\mathbf{Au}}_{k}$ for all $k\in [N]$ (rather than just the ${\lambda}_{k}>0$ modes). This choice gives $\mathbf{A}={\sum}_{k=1}^{N}{\stackrel{~}{\mathbf{u}}}_{k}{\mathbf{u}}_{k}^{\top}$. Thus, the space of codes $\stackrel{~}{\mathbf{r}}(\mathit{\theta})$ with equivalent kernels to $\mathbf{r}(\mathit{\theta})\cdot \mathbf{r}({\mathit{\theta}}^{\prime})$ generated through linear transformations is equivalent to all possible orthogonal transformations of the original code $\{\mathbf{Qr}(\mathit{\theta}):{\mathbf{QQ}}^{\top}={\mathbf{Q}}^{\top}\mathbf{Q}=\mathbf{I}\}$. ∎
Effect of noise on RROS symmetry
The random rotation and optimal shift (RROS) operations introduced in the main text preserve generalization performance under the assumption of a deterministic neural code. However, for noisy codes, the presence of RROS symmetry is dependent on the noise distribution. Below we discuss two commonly analyzed distributions: the Gaussian distribution and the Poisson distribution. For Gaussian noise, the RROS operations preserve the generalization performance and the local Fisher information. However, if noise is constrained to be Poisson then RROS operations do not preserve generalization or Fisher information.
First, we will analyze stimulus dependent Gaussian noise, where generalization performance is preserved under rotations and baseline shifts. Note that if the code at $\mathit{\theta}$ obeyed $\mathbf{r}(\mathit{\theta})\sim \mathcal{N}(\overline{\mathbf{r}}(\mathit{\theta}),{\mathbf{\Sigma}}_{n}(\mathit{\theta}))$, then the rotated and shifted code follows $\mathbf{Q}\mathbf{r}(\mathit{\theta})+\mathit{\delta}\sim \mathcal{N}(\mathbf{Q}\overline{\mathbf{r}}(\mathit{\theta})+\mathit{\delta},\mathbf{Q}{\mathbf{\Sigma}}_{n}(\mathit{\theta}){\mathbf{Q}}^{\mathrm{\top}})$. This rotated and shifted code $\mathbf{Qr}(\mathit{\theta})+\mathit{\delta}$, when centered, will exhibit identical generalization performance as the original code. This is true both for learning from a trial averaged or nontrial averaged code. In the case of Gaussian noise on a centered code, the dataset transforms under a rotation as $\mathcal{D}=\{{\mathbf{r}}_{\mu},{y}_{\mu}\}\to {\mathcal{D}}^{\mathrm{\prime}}=\{\mathbf{Q}{\mathbf{r}}_{\mu},{y}_{\mu}\}$. The optimal weights for a linear model similarly transform as $\mathbf{w}(\mathcal{D})\to \mathbf{Q}\mathbf{w}(\mathcal{D})$. Under these transformations the predictor on test point $\mathit{\theta}$ is unchanged since
Further, the local Fisher information matrix is $\mathbf{I}(\mathit{\theta})=\frac{\partial \overline{\mathbf{r}}{(\mathit{\theta})}^{\top}}{\partial \mathit{\theta}}{\mathbf{\Sigma}}_{n}^{1}(\mathit{\theta})\frac{\partial \overline{\mathbf{r}}(\mathit{\theta})}{\partial {\mathit{\theta}}^{\top}}+\frac{1}{2}\text{Tr}{\mathbf{\Sigma}}_{n}^{1}(\mathit{\theta})\frac{\partial {\mathbf{\Sigma}}_{n}(\mathit{\theta})}{\partial \mathit{\theta}}{\mathbf{\Sigma}}_{n}^{1}(\mathit{\theta})\frac{\partial {\mathbf{\Sigma}}_{n}(\mathit{\theta})}{\partial {\mathit{\theta}}^{\top}}$ is unchanged under the transformation $\mathbf{r}\to \mathbf{Qr}+\mathit{\delta}$. Under this transformation, the covariance simply transforms linearly ${\mathbf{\Sigma}}_{n}\to \mathbf{Q}{\mathbf{\Sigma}}_{n}(\mathit{\theta}){\mathbf{Q}}^{\top}$ and the $\mathbf{Q}$ matrices will annihilate under the trace. This shows that, for some noise models, our assumption that rotations and baseline shifts preserve generalization performance will be valid.
However, for Poisson noise, where the variance is tied to the mean firing rate, the RROS operations will not preserve noise structure or information content. The Fisher information at scalar stimulus $\theta $ for a Poisson neuron is $I(\theta )=\frac{{\overline{r}}^{\prime}{(\theta )}^{2}}{\overline{r}(\theta )}$. A baseline shift $r\to r+\delta $ to the tuning curve will not change the numerator since the derivative of the tuning curve is invariant to this transformation, but it will increase the denominator.
Necessary conditions for optimally sparse codes
Next we argue why optimally sparse codes should be lifetime and population selective. We consider the following optimization problem: find a nonnegative neural responses $\mathbf{S}\in {\mathbb{R}}^{N\times P}$ and baseline vector $\mathit{\delta}\in {\mathbb{R}}^{N}$ so that baseline subtracted responses $\mathbf{R}=\mathbf{S}\mathit{\delta}{\mathbf{1}}^{\top}$ realize a desired inner product kernel $\mathbf{K}\in {\mathbb{R}}^{P\times P}$ and have minimal total firing. This is equivalent to finding the most metabolically efficient code among the space of codes with equivalent inductive bias. Mathematically, we formulate this problem as
To enforce the constraints for the definition of the kernel and the nonnegativity of the responses, we introduce the following Lagrangian
where 1 is the vector containing all ones, the Lagrange multiplier matrix $\mathbf{A}$ enforces the definition of the kernel and the KKT multiplier matrix $\mathbf{V}$ enforces the nonnegativity constraints for each element of $\mathbf{S}$. The KKT conditions require that any local optimum of the objective would have to satisfy the following equations (Kuhn and Tucker, 2014)
where $\odot $ denotes the elementwise Hadamard product. Using the complementary slackness condition $\mathbf{S}\odot \mathbf{V}=\mathbf{0}$, and the first optimality condition $\frac{\partial \mathcal{L}}{\partial \mathbf{S}}=\mathbf{0}$, we have
Therefore, for any neuronstimulus pair $(i,\mu )$, either ${S}_{i\mu}=0$ or ${\sum}_{\nu \in [P]}({S}_{i\nu}{\delta}_{i}){A}_{\nu \mu}=1$. Further, under the condition that K is full rank, we conclude that for any stimulus μ, ${\sum}_{\nu \in [P]}{A}_{\mu \nu}=0$ from the equation $\frac{\partial \mathcal{L}}{\partial \mathit{\delta}}=\mathbf{0}$. Let $\mathcal{I}}_{i}=\{\mu \in [P]:{S}_{i\mu}>0\$ represent the set of stimuli for which neuron $i$ fires. We will call this the receptive field set for neuron $i$. Let $\mathbf{B}}_{(i)}\in {\mathbb{R}}^{P\times P$ have entries
where the matrix ${\mathbf{A}}_{(i)}$ is the ${\mathcal{I}}_{i}\times {\mathcal{I}}_{i}$ minor of $\mathbf{A}$ obtained by taking all rows and columns with indices $\mu ,\nu \in {\mathcal{I}}_{i}$, and ${\mathbf{A}}^{+}$ denotes pseudoinverse of $\mathbf{A}$. Then the $i$th neuron’s tuning curve is a function of the index set ${\mathcal{I}}_{i}$ the baseline ${\delta}_{i}$ and the neuronindependent $P\times P$ matrix $\mathbf{A}$. The nonnegativity constraint for neuron $i$’s tuning curve implies that ${S}_{i\mu}=\sum _{\nu \in {\mathcal{I}}_{i}}{B}_{(i),\mu \nu}[{\delta}_{i}\sum _{\gamma \in [P]}{A}_{\nu \gamma}+1]>0$ for all $\mu \in {\mathcal{I}}_{i}$. To satisfy the definition of the kernel, we have the following constraint on the matrix $\mathbf{A}$, the index sets ${\mathcal{I}}_{i}$ and baselines ${\delta}_{i}$
This equation implictly defines the index sets ${\mathcal{I}}_{i}$ the baselines ${\delta}_{i}$ and the KKT matrix $\mathbf{A}$. We see that, in order to fit an arbitrary kernel, the receptive field sets $\{{\mathcal{I}}_{i}\}$ and baselines ${\delta}_{i}$ for each neuron must be sufficiently diverse since otherwise only a low rank kernel matrix can be achieved from the optimally sparse code. As a concrete example, suppose that ${\mathcal{I}}_{i}=\mathcal{I}$ so that ${\mathbf{V}}_{(i)}=\mathbf{V}$ and ${\delta}_{i}=\delta $ for all $i$. For example, this could occur if each neuron fired for every possible stimulus. In this case, the kernel would be rank one: $\mathbf{K}=N(\mathbf{s}(\mathcal{I},\delta ,\mathbf{A})\delta \mathbf{1}){(\mathbf{s}(\mathcal{I},\delta ,\mathbf{A})\delta \mathbf{1})}^{\top}$. In order to achieve a higher rank code there must be sufficient diversity of the receptive fields ${\mathcal{I}}_{i}$. Thus the only way for optimally sparse codes to realize high rank kernels $\mathbf{K}$ is to have neurons to have different receptive field sets ${\mathcal{I}}_{i}$. The necessary optimality conditions thus reveal a preference for sparse neural tuning curves to have high lifetime sparseness; to achieve diverse index sets ${\mathcal{I}}_{i}$, any given neuron will fire only for a unique subset of the possible stimuli.
Impact of neural noise and unlearnable targets on learning
While our analysis so far has focused on deterministic population codes, our theory can be extended to neural populations which exhibit variability in responses to identical stimuli. For each stimulus $\mathit{\theta}$, we let the population response $\mathbf{r}(\mathit{\theta})$ be a random vector with mean $\overline{\mathbf{r}}(\mathit{\theta})={\u27e8\mathbf{r}(\mathit{\theta})\u27e9}_{\mathbf{r}\mathit{\theta}}$ and covariance ${\mathbf{\Sigma}}_{n}(\mathit{\theta})={\u27e8(\mathbf{r}(\mathit{\theta})\overline{\mathbf{r}}(\mathit{\theta})){(\mathbf{r}(\mathit{\theta})\overline{\mathbf{r}}(\mathit{\theta}))}^{\top}\u27e9}_{\mathbf{r}\mathit{\theta}}$.
The (deterministic) target function can be decomposed in terms of the mean response as $y(\mathit{\theta})={\mathbf{w}}^{*}\cdot \overline{\mathbf{r}}(\mathit{\theta})$ (the usual decomposition $y={\mathbf{w}}^{*}\cdot \mathbf{r}(\mathit{\theta})$ gives an unphysical target function which fluctuates with the variability in neural responses). For a given configuration of weights $\mathbf{w}$, the generalization error (which is an average over the joint distribution of $\mathbf{r},\mathit{\theta}$) is determined only by the signal ${\mathbf{\Sigma}}_{s}={\u27e8\overline{\mathbf{r}}(\mathit{\theta})\overline{\mathbf{r}}{(\mathit{\theta})}^{\top}\u27e9}_{\mathit{\theta}}$ and noise ${\mathbf{\Sigma}}_{n}={\u27e8{\mathbf{\Sigma}}_{n}(\mathit{\theta})\u27e9}_{\mathit{\theta}}$ correlation matrices:
where we utilized the fact that ${\u27e8\mathbf{r}(\mathit{\theta})\overline{\mathbf{r}}(\mathit{\theta})\u27e9}_{\mathbf{r}\theta}=0$ to eliminate the crossterm. The two terms in the final expression can be thought of as a biasvariance decomposition over the noise in neural responses. The minimum achievable loss can be obtained by differentiation of the generalization error expression with respect to $\mathbf{w}$, giving ${E}_{g}^{*}={\mathbf{w}}^{*}{\mathbf{\Sigma}}_{n}{({\mathbf{\Sigma}}_{s}+{\mathbf{\Sigma}}_{n})}^{1}{\mathbf{\Sigma}}_{s}{\mathbf{w}}^{*}$. We note that any noise correlation matrix with noise orthogonal to coding direction ${\mathbf{\Sigma}}_{n}{\mathbf{w}}^{*}=0$ will give the minimal (zero) asymptotic error. Alignment of the noise ${\mathbf{\Sigma}}_{n}$ with ${\mathbf{w}}^{*}$ gives higher asymptotic error.
In addition to the irreducible error, the presence of neural noise can alter the learning curve at finite $P$. An analytical study of this is difficult, which we leave for future work. We numerically study the effect of neural variability on generalization performance in the orientation discrimination tasks for nontrialaveraged Mouse V1 code in Appendix 1—figure 1 . We note that the generalization error is worse at each finite value of $P$ when compared to trial averaged (noise free) learning curves. We varied the regularization parameter and did not find an obvious nonzero optimal weight decay $\lambda $, consistent with small noise levels.
Neural noise is not the only phenomenon that can degrade task learning. Codes which are incapable of expressing the target function through linear readouts are also susceptible to overfitting. As explained in Canatar et al., 2021, the components of the target function that are inexpressible act as a source of noise on the learning process which can overfit this noise. Such a scenario can occur, for example, when the readout neuron only gets input from a sparse subset of the coding neural population (Seeman et al., 2018). We show in Appendix 1—figure 1CD that using subsampled populations of size $N$ can indeed lead to a regime where more data can hurt performance leading to an overfitting error peak, a subsequent nonvanishing asymptotic error, and an optimal weight decay parameter $\lambda $. This phenomenon is known as double descent in machine learning literature (Belkin et al., 2019; Mei and Montanari, 2020; Canatar et al., 2021). At small $N$, these codes are not sufficiently expressive to learn the target function through linear readout. The overfitting peak occurs near the interpolation threshold, the largest value of $P$ where all training sets could be perfectly fit in the $\lambda \to 0$ limit (Canatar et al., 2021). At infinite $P$, generalization error asymptotes to the amount of unexplained variance in the target function.
Data availability
Mouse V1 neuron responses to orientation gratings and preprocessing code were obtained from a publicly available dataset: https://github.com/MouseLand/stringeretal2019. Responses to ImageNet images and preprocessing code were obtained from another publicly available dataset, https://github.com/MouseLand/stringerpachitariuetal2018b. The code generated by the authors for this paper is also available https://github.com/PehlevanGroup/sample_efficient_pop_codes, (copy archived at swh:1:rev:6cd4f0fe7043ae214dd682a9dc035a497ffa2d61).

FigshareRecordings of ten thousand neurons in visual cortex in response to 2,800 natural images.https://doi.org/10.25378/janelia.6845348.v4

FigshareRecordings of ~20,000 neurons from V1 in response to oriented stimuli.https://doi.org/10.25378/janelia.8279387.v3
References

ConferenceDo retinal ganglion cells project natural scenes to their principal subspace and whiten them?2016 50th Asilomar Conference on Signals, Systems and Computers.https://doi.org/10.1109/ACSSC.2016.7869658

Perceptual learning of simple stimuli modifies stimulus representations in posterior inferior temporal cortexJournal of Cognitive Neuroscience 26:2187–2200.https://doi.org/10.1162/jocn_a_00641

Spatiotemporal energy models for the perception of motionJournal of the Optical Society of America. A, Optics and Image Science 2:284–299.https://doi.org/10.1364/josaa.2.000284

The reverse hierarchy theory of visual perceptual learningTrends in Cognitive Sciences 8:457–464.https://doi.org/10.1016/j.tics.2004.08.011

Generation of random orthogonal matricesSIAM Journal on Scientific and Statistical Computing 8:625–629.https://doi.org/10.1137/0908055

What does the retina know about natural scenes?Neural Computation 4:196–210.https://doi.org/10.1162/neco.1992.4.2.196

Some informational aspects of visual perceptionPsychological Review 61:183–193.https://doi.org/10.1037/h0054663

Neural correlations, population coding and computationNature Reviews. Neuroscience 7:358–366.https://doi.org/10.1038/nrn1888

BookPossible Principles Underlying the Transformation of Sensory MessagesCambridge University.

Benign overfitting in linear regressionPNAS 117:30063–30070.https://doi.org/10.1073/pnas.1907378117

Local diversity and finescale organization of receptive fields in mouse visual cortexThe Journal of Neuroscience 31:18506–18521.https://doi.org/10.1523/JNEUROSCI.297411.2011

ConferenceSpectrum dependent learning curves in kernel regression and wide neural networksProceedings of the 37th International Conference on Machine Learning of Proceedings of Machine Learning Research. pp. 1024–1034.

ConferenceSelfConsistent Dynamical Field Theory of Kernel Evolution in Wide Neural NetworksAdvances In Neural Information Processing Systems.

Deep neural networks rival the representation of primate it cortex for core visual object recognitionPLOS Computational Biology 10:e1003963.https://doi.org/10.1371/journal.pcbi.1003963

Principal component analysis of neuronal ensemble activity reveals multidimensional somatosensory representationsJournal of Neuroscience Methods 94:121–140.https://doi.org/10.1016/s01650270(99)001302

Measuring and interpreting neuronal correlationsNature Neuroscience 14:811–819.https://doi.org/10.1038/nn.2842

Dimensionality reduction for largescale neural recordingsNature Neuroscience 17:1500–1509.https://doi.org/10.1038/nn.3776

BookTheoretical Neuroscience: Computational and Mathematical Modeling of Neural SystemsThe MIT Press.

ConferenceImagenet: A largescale hierarchical image databaseIn 2009 IEEE conference on computer vision and pattern recognition.

ConferenceReservoir computing meets recurrent kernels and structured transformsNeurIPS Proceedings.

Representation is representation of similaritiesThe Behavioral and Brain Sciences 21:449–467.https://doi.org/10.1017/s0140525x98001253

On simplicity and complexity in the brave new world of largescale neuroscienceCurrent Opinion in Neurobiology 32:148–155.https://doi.org/10.1016/j.conb.2015.04.003

Physiological correlates of perceptual learning in monkey V1 and V2Journal of Neurophysiology 87:1867–1888.https://doi.org/10.1152/jn.00690.2001

Regularization theory and neural networks architecturesNeural Computation 7:219–269.https://doi.org/10.1162/neco.1995.7.2.219

How noise contributes to contrast invariance of orientation tuning in cat visual cortexThe Journal of Neuroscience 22:5118–5128.https://doi.org/10.1523/JNEUROSCI.221205118.2002

Predictive codingWiley Interdisciplinary Reviews. Cognitive Science 2:580–593.https://doi.org/10.1002/wcs.142

Weight perturbation: an optimal architecture and learning technique for analog vlsi feedforward and recurrent multilayer networksIEEE Transactions on Neural Networks 3:154–157.https://doi.org/10.1109/72.105429

Advances in Neural Information Processing SystemsNeural tangent kernel: convergence and generalization in neural networks, Advances in Neural Information Processing Systems, Curran Associates, Inc.

ConferenceSGD on neural networks learns functions of increasing complexityIn Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019.

Representational similarity analysisconnecting the branches of systems neuroscienceFrontiers in Systems Neuroscience 2:4.https://doi.org/10.3389/neuro.06.004.2008

BookNonlinear programmingIn: Kuhn HW, editors. Traces and Emergence of Nonlinear Programming. Springer. pp. 1–4.https://doi.org/10.1007/9783034804394

Content and cluster analysis: assessing representational similarity in neural systemsPhilosophical Psychology 13:47–76.https://doi.org/10.1080/09515080050002726

Building machines that learn and think like peopleThe Behavioral and Brain Sciences 40:e253.https://doi.org/10.1017/S0140525X16001837

ConferenceDeep neural networks as gaussian processesIn International Conference on Learning Representations.

Random synaptic feedback weights support error backpropagation for deep learningNature Communications 7:13276.https://doi.org/10.1038/ncomms13276

Learning curves of generic features maps for realistic datasets with a teacherstudent modelAdvances in Neural Information Processing Systems 34:18137–18151.https://doi.org/10.1088/17425468/ac9825

The generalization error of random features regression: precise asymptotics and the double descent curveCommunications on Pure and Applied Mathematics 75:667–766.https://doi.org/10.1002/cpa.22008

Adaptive tuning curve widths improve sample efficient learningFrontiers in Computational Neuroscience 14:12.https://doi.org/10.3389/fncom.2020.00012

XVI. Functions of positive and negative type, and their connection the theory of integral equationsPhilosophical Transactions of the Royal Society of London. Series A, Containing Papers of A Mathematical or Physical Character 209:415–446.https://doi.org/10.1098/rsta.1909.0016

Neural noise can explain expansive, powerlaw nonlinearities in neural response functionsJournal of Neurophysiology 87:653–659.https://doi.org/10.1152/jn.00425.2001

ConferenceOn 1/n neural representation and robustnessAdvances in Neural Information Processing Systems 33.

Highly selective receptive fields in mouse visual cortexThe Journal of Neuroscience 28:7520–7536.https://doi.org/10.1523/JNEUROSCI.062308.2008

Energy limitation as a selective pressure on the evolution of sensory systemsThe Journal of Experimental Biology 211:1792–1804.https://doi.org/10.1242/jeb.017574

Robustness of spike deconvolution for neuronal calcium imagingThe Journal of Neuroscience 38:7976–7985.https://doi.org/10.1523/JNEUROSCI.333917.2018

BookRecordings of 20,000 Neurons from V1 in Response to Oriented StimuliAmerican Physiological Society.

Why do similarity matching objectives lead to hebbian/antihebbian networks?Neural Computation 30:84–124.https://doi.org/10.1162/neco_a_01018

Decorrelation and efficient coding by retinal ganglion cellsNature Neuroscience 15:628–635.https://doi.org/10.1038/nn.3064

ConferenceLie group methods for optimization with orthogonality constraintsIn International Conference on Independent Component Analysis and Signal Separation.

The contribution of spike threshold to the dichotomy of cortical simple and complex cellsNature Neuroscience 7:1113–1122.https://doi.org/10.1038/nn1310

ConferenceOn the spectral bias of neural networksIn International Conference on Machine Learning.

ConferenceA generalized representer theoremIn Proceedings of the 14th Annual Conference on Computational Learning Theory and 5th European Conference on Computational Learning Theory, COLT ’01/EuroCOLT ’01.

BookLearning with Kernels: Support Vector Machines, Regularization, Optimization, and BeyondMIT press.

Neural basis of a perceptual decision in the parietal cortex (area lip) of the rhesus monkeyJournal of Neurophysiology 86:1916–1936.https://doi.org/10.1152/jn.2001.86.4.1916

ConferenceNeural tangent kernel eigenvalues accurately predict generalizationICLR 2022 Conference.

A model of neuronal responses in visual area MTVision Research 38:743–761.https://doi.org/10.1016/s00426989(97)001831

Natural image statistics and neural representationAnnual Review of Neuroscience 24:1193–1216.https://doi.org/10.1146/annurev.neuro.24.1.1193

ConferenceApproximate learning curves for Gaussian processes9th International Conference on Artificial Neural Networks.https://doi.org/10.1049/cp:19991148

BookGaussian process regression with mismatched modelsIn: Dietterich T, Becker S, Ghahramani Z, editors. Advances in Neural Information Processing Systems. MIT Press. pp. 1–2.

Predictive coding: a fresh view of inhibition in the retinaProceedings of the Royal Society of London. Series B, Biological Sciences 216:427–459.https://doi.org/10.1098/rspb.1982.0085

Continuous analogues of matrix factorizationsProceedings of the Royal Society A 471:20140585.https://doi.org/10.1098/rspa.2014.0585

ConferenceDeep learning generalizes because the parameterfunction map is biased towards simple functionsIn International Conference on Learning Representations.

Dynamics of spatial frequency tuning in mouse visual cortexJournal of Neurophysiology 107:2937–2949.https://doi.org/10.1152/jn.00022.2012

Characterizing the sparseness of neural codesNetwork: Computation in Neural Systems 12:255.https://doi.org/10.1080/net.12.3.255.270

The lack of a priori distinctions between learning algorithmsNeural Computation 8:1341–1390.https://doi.org/10.1162/neco.1996.8.7.1341

The effect of perceptual learning on neuronal responses in monkey visual area v4The Journal of Neuroscience 24:1617–1626.https://doi.org/10.1523/JNEUROSCI.444203.2004

ConferenceTensor programs iv: Feature learning in infinitewidth neural networksIn International Conference on Machine Learning.

ConferenceUnderstanding deep learning requires rethinking generalizationIn 5th Int. Conf. on Learning Representations (ICLR 2017).
Decision letter

Thomas SerreReviewing Editor; Brown University, United States

Michael J FrankSenior Editor; Brown University, United States

Fabio AnselmiReviewer; Baylor College of Medicine, United States

Jeff BeckReviewer; Duke University, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Population codes enable learning from few examples by shaping inductive bias" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Michael Frank as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Fabio Anselmi (Reviewer #1); Jeff Beck (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
The consensus among all 3 reviewers is that the manuscript describes an insightful and novel mathematical framework for evaluating the ability of a downstream linear decoder to learn new stimulusresponse mappings from relatively few training examples. However, there is also broad consensus among the 3 reviewers that (1) the manuscript is too technical and written in a way that may not be palatable to the eLife readership (to paraphrase one of the reviewers "I don't think that many experimentalists will immediately see the sheer utility of having a tool like this in their arsenal" and that 2) some of the assumptions made are overly simplistic and the theory needs to be extended with more realistic neural assumptions regarding noise, readout and kernels.
Wrt (1), before the manuscript can be accepted for publication, it is necessary that the authors expand on possible applications to neural and behavioral experiments. The authors also need to distill the math better, provide more intuitive explanations, and reduce the jargon. Wrt (2) R2 provided a detailed lists of comments with suggestions for more realistic assumptions.
Reviewer #1 (Recommendations for the authors):
The authors provide a simple and clear way to understand an aspect of the implicit bias of a neural population code linking it with wellknown machine learning methods and concepts such as kernel regression, sample complexity and efficiency.
Although the mathematical results the authors employ are not novel, the way they apply them to the problem of neural coding is novel and interesting to a broad audience.
In particular, the computational neuroscience community can benefit from this work being it is one of the few dealing with the impact of the model implicit bias in explaining real data.
Reviewer #2 (Recommendations for the authors):
It is my opinion that the principle utility of this approach lies in its ability to identify the set of 'easily learnable' stimulusresponse mappings from neural data which makes strong behavioral predictions that can be easily evaluated. I envision a simple experiment in which empirically obtained kernel functions are used to rank stimulusresponse mappings according to their learnability which can then be plotted against measures of performance like the observed learning rate and saturated performance. Because kernel functions are empirically obtained, there is even the potential for meaningful crossspecies comparisons. If behaviorally validated, one could also use this approach to label cortical populations by the set of easily learned stimulusresponse mappings for that population. This allows for the identification of taskrelevant neurons or regions which can be subsequently manipulated to enhance or degrade learning rates.
Of course, any theoretical approach is only as good as the underlying assumptions and so while the primary strength is the simplicity and generality of this approach, the primary weakness is its neglect of some very real and very relevant aspects of neural data in particular and statistical learning in general. In particular, the three principle limitations of this work are tied to its reliance on the assumptions that (1) neurons are noiseless, (2) decoders are linear, and (3) learned weights are unbiased.
(1) Within this framework, a realistic stimulusdependent noise model can be easily introduced and its effects on the kernel and set of easily learned stimulusresponse mappings investigated. So while the kernel would be substantially altered via the addition of a realistic noise model, the applications of the approach outlined above would not be affected. The same cannot be said for the efficient coding application described in this manuscript. There, the authors note that rotations and constant shifts of neural activity do not affect the kernel and thus do not affect the generalization error. This kernel invariance is not present when a nontrivial (i.e. nonisotropic) noise model is added. For example, suppose that neurons are independent and Poisson so that noise scales with the mean of the neural response. In this case, adding a baseline firing rate to a population of unimodal neurons representing orientation necessarily reduces the information content of the population while rotations can affect the fidelity with which certain stimulus values are represented. It is important to note, however, that while this particular efficiency result is not compelling, I believe that it is possible to perform a similar analysis that takes into account realistic noise models and focuses on a broad set of 'biologically plausible' kernels instead of particular invariant ones. For example, one could consider noise covariance structures with differential correlations (MorenoBote 2014). Since the magnitude of differential correlations controls the redundancy of the population code this would enable an analysis of the role of redundancy in suppressing (or enhancing) generalization error.
(2) Similarly, the linearity assumption is somewhat restrictive. Global linear decoders of neural activity are known to be highly inefficient and completely fail when decoding orientation in the primary visual cortex in the presence of contrast fluctuations. This is because contrast modulates the amplitude of the neural response and doubling the amplitude means doubling an estimate obtained from a linear decoder even when the underlying orientation has not changed. While the contrast issue could be partially addressed by simply considering normalized neural responses, it is not yet clear how to extend this approach to account for other sources of neural variability and covariability that cause global linear decoders to fail so badly.
(3) This analysis relies on the assumption that decoder weights learned in the presence of finite data are efficient and unbiased. This assumption is problematic particularly when it comes to inductive bias and generalization error. This is because a standard way to reduce generalization error is to introduce bias into the learned decoder weights through a penalization scheme that privileges decoder weights with small magnitudes. This kind of regularization is particularly important when neurons are noisy. Fortunately, this issue could be addressed by parameterizing changes in the kernel function by the degree and type of regularization potentially leading to a more general result.
Finally, I would like to conclude by explicitly stating that while the limitations imposed by the assumptions listed above temper my enthusiasm in regards to conclusions drawn in this work, I do not believe there is some fundamental problem with the general theoretical framework. Indeed, items 1 and 3 above can be easily addressed through straightforward extensions of the authors approach and I look forward to their implementation. Item 2 is a bit more troublesome, but my intuition tells me that an informationtheoretic extension based upon Fisher information may be capable of eliminating all three of these limiting assumptions by exploiting the relationship between FI(\theta) and FI(y=f(\theta)).
Ultimately, all I felt the need to say is in the public part of the review. But I wanted to use this space to insure that it was clear that I feel that this approach has a lot of potential and that I very much hope to see it extended/generalized.
If I were to make any suggestions they would be:
(1) a discussion of the consequences of biologically plausible noise models. I believe that this only requires simply augmenting K by a stimulus dependent noise covariance matrix.
(2) using this addition of a noise model to enable a discussion of the role of redundancy in generalization. This could be accomplished by considering perturbations to the noise covariance matrix that introduce differential correlations (MorenoBote 2014) of varying magnitude. More differential correlations means more redundancy and, I suspect, better generalization.
(3) the addition of a paragraph or two outlining some of the other ways that this learnability/generalization measure could be of use to physiologists and behavioral scientists. Training a mouse to do anything more complicated simple orientation detection is challenging to say the least. So having any tool that can be used to identify functions $y=f(\theta)$ that a mouse has a good chance of learning quickly would be highly advantageous. Similarly, by strategically subsampling neurons by tuning properties I suspect it may be possible to identify subpopulations of neurons that are particular important for learning certain functions. This is also cool because it allows physiologists to then target those populations for manipulation.
Of these (3) is probably the most important for this manuscript, because I don't think that many experimentalists will immediately see the sheer utility of having a tool like this in their arsenal.
Reviewer #3 (Recommendations for the authors):
The manuscript presents a theory of generalization performance in deterministic population codes, that applies to the case of small numbers of training examples. The main technical result, as far as I understand, is that generalization performance (the expected classification or regression error) of a population code depends exclusively on the 'kernel', i.e. a measure of the pairwise similarity between population activity patterns corresponding to different inputs. The main conceptual results are that, using this theory, one can understand the inductive biases of the code just from analyzing the kernel, particularly the top eigenfunctions; and that sampleefficient learning (low generalization performance with few samples) depends on whether the task is aligned with the population's inductive bias, that is, whether the target function (i.e. the true map from inputs to outputs) is aligned with the top eigenfunctions of the kernel. For instance, in mouse V1 data, they show that the top eigenfunctions correspond to low frequency functions of visual orientation (i.e. functions that map a broad range of similar orientations to similar output value), and that consistent with the theory, the generalization performance for small sample sizes is better for tasks defined by low frequency target functions. In my opinion, perhaps the most significant finding from a neuroscience perspective, is that the conditions for good generalization at low samples are markedly different from those in the largesample asymptotic regime studies in Stringer et al. 2018 Nature: rather than a tradeoff between highdimensionality and differentiability proposed by Stringer et al., this manuscript shows that in the lowsample regime such codes can be disadvantageous for small sample sizes, that differentiability is not required, that the top eigenvalues matter more than the tail of the spectrum, and what matters is the alignment between the task and the top eigenfunctions. The authors propose sampleefficient learning/generalization as a new principle of neural coding, replacing or complementing efficient coding.
Overall, in my opinion this is a remarkable manuscript, presenting truly innovative theory with somewhat limited but convincing application to neural data. My main concern is that this is highly technical, dense, and long; the mathematical proofs for the theory are buried in the supplement and require knowledge of disparate techniques from statistical physics. Although some of that material on the theory of generalization is covered in previous publications by the authors, it was not clear to me if that is true for all of the technical results or only some.
Fixed population code, learnable linear readout: the authors acknowledge in the very last sentences of the manuscript that this is a limitation, given that neural tuning curves (the population neural code) are adaptable. I imagine extending the theory to both learnable codes and learnable readouts is hard and I understand it's beyond the scope of this paper. But perhaps the authors could motivate and discuss this choice, not just because of its mathematical convenience but also in relation to actual neural systems: when are these assumptions expected to be a good approximation of the real system?
The analysis of V1 data, showing a bias for lowfrequency functions of orientation is convincing. But it could help if the authors provided some considerations on the kind of ethological behavioral context where this is relevant, or at least the design of an experimental behavioral task to probe it. Also related, it would be useful to construct and show a counterexample, a synthetic code for which the highfrequency task is easier.
Line 519, data preprocessing: related to the above, is it possible that binning together the V1 responses to gratings with different orientations (a range of 3.6 deg per bin, if I understood correctly) influences the finding of a lowfrequency bias?
I found the study of invariances interesting, where the theory provides a normative prediction for the proportion of simple and complex cells. However, I would suggest the authors attempt to bring this analysis a step closer to the actual data: there are no pure simple and complex cells, usually the classification is based on responses to gratings phases (F1/F0) and real neurons take a continuum of values. Could the theory qualitatively predict that distribution?
https://doi.org/10.7554/eLife.78606.sa1Author response
Essential revisions:
The consensus among all 3 reviewers is that the manuscript describes an insightful and novel mathematical framework for evaluating the ability of a downstream linear decoder to learn new stimulusresponse mappings from relatively few training examples. However, there is also broad consensus among the 3 reviewers that (1) the manuscript is too technical and written in a way that may not be palatable to the eLife readership to paraphrase one of the reviewers "I don't think that many experimentalists will immediately see the sheer utility of having a tool like this in their arsenal" and that (2) some of the assumptions made are overly simplistic and the theory needs to be extended with more realistic neural assumptions regarding noise, readout and kernels.
Wrt (1), before the manuscript can be accepted for publication, it is necessary that the authors expand on possible applications to neural and behavioral experiments. The authors also need to distill the math better, provide more intuitive explanations, and reduce the jargon. Wrt (2) R2 provided a detailed lists of comments with suggestions for more realistic assumptions.
We thank the editor for these points of advice. We respond to these points in detail below, in our pointbypoint rebuttal to the reviewers. Here, we provide a summary of the measures taken to address each of these problems.
– Writing style too technical: We attempted to eliminate mathematical jargon and keep only the crucial mathematical components in the main text. The detailed discussion of the RNN model in the main text was shortened and the reference to Lie groups was removed. We added a paragraph at the beginning of the Results section where we define our mathematical notation explicitly. We also removed unnecessary detail about the V1 Gabor model from the main text. We added more intuition about what a kernel is and how it relates to distances between points in neural space. We also tried adding more intuition about the spectral bias result by giving the simpler expression for the average learned function in Methods Theory of Generalization starting at line 617.
– Contextualize Results for Experimental Audience: To describe how our theory of inductive bias may be useful for experimental neuroscience and could lead to new neural and behavioral experiments, we added a new paragraph at the end of the discussion.
– Assumptions Overly Simplistic: To address our previous model’s failure to model neural noise, we added a Appendix Section Impact of Neural Noise and Unlearnable Targets on Learning, where we discuss effects of noise on the learning curve. We allow for L2 regularization of the readout weights with a ridge penalty λ and show that there can be an optimal λ to prevent overfitting. We discuss nonlinear readouts in greater detail in App Typical Case Analysis of Nonlinear Readouts, showing that the kernel still governs the learned function for a wider set of learning problems. We discuss recent work, which computes learning curves for nonlinear readouts with a similar technique. Lastly, we provide a more realistic model of V1 tuning curves where neurons are neither completely simple or complex, but have a distribution of F1/F0 values (see App Energy Model with Partially PhaseSelective Cells). We discuss limitations of our work brought up in the reviews in our new discussion. Overall, we view extensions to address these present limitations as promising areas of future research.
Reviewer #1 (Recommendations for the authors):
The authors provide a simple and clear way to understand an aspect of the implicit bias of a neural population code linking it with wellknown machine learning methods and concepts such as kernel regression, sample complexity and efficiency.
Although the mathematical results the authors employ are not novel, the way they apply them to the problem of neural coding is novel and interesting to a broad audience.
In particular, the computational neuroscience community can benefit from this work being it is one of the few dealing with the impact of the model implicit bias in explaining real data.
We thank the reviewer for insightful comments and the support. While addressing these comments, we had the chance to clarify some subtle points.
Reviewer #2 (Recommendations for the authors):
It is my opinion that the principle utility of this approach lies in its ability to identify the set of 'easily learnable' stimulusresponse mappings from neural data which makes strong behavioral predictions that can be easily evaluated. I envision a simple experiment in which empirically obtained kernel functions are used to rank stimulusresponse mappings according to their learnability which can then be plotted against measures of performance like the observed learning rate and saturated performance. Because kernel functions are empirically obtained, there is even the potential for meaningful crossspecies comparisons. If behaviorally validated, one could also use this approach to label cortical populations by the set of easily learned stimulusresponse mappings for that population. This allows for the identification of taskrelevant neurons or regions which can be subsequently manipulated to enhance or degrade learning rates.
We thank the reviewer for appreciating the possible applications of our work for future empirical studies on stimulusresponse learning and inductive bias in neural systems. We are also excited by these possibilities for neuroscience experiments. We added a this as a discussion item where we suggest some ways our framework could be used in future empirical work. We discuss this addition further below.
Of course, any theoretical approach is only as good as the underlying assumptions and so while the primary strength is the simplicity and generality of this approach, the primary weakness is its neglect of some very real and very relevant aspects of neural data in particular and statistical learning in general. In particular, the three principle limitations of this work are tied to its reliance on the assumptions that (1) neurons are noiseless, (2) decoders are linear, and (3) learned weights are unbiased.
We thank this reviewer for encouraging us to explore additional realism in our model of sample efficient learning and acknowledge the limitations of our current analysis. Below we summarize each of the changes we made in response.
1. We updated the manuscript to allow for biased weights through ridge regression (regression with an additional L2 norm penalty λw^{2} on the weights) with ridge parameter λ. The effect of this parameter on generalization is easy to identify as it merely alters the implicit definition of $k\text{}to\text{}k=\text{}\lambda +k\sum _{k}^{}\frac{\text{\lambda k}}{{\lambda}_{k}P+k}$ (see Equation 11).
2. We provide a preliminary analysis of neuron noise in the new Appendix Section Impact of Neural Noise and Unlearnable Targets on Learning. Although we currently do not have a dataaveraged theory of generalization (which is involves a very complicated computation), we did attempt empirically to explore the effect of neural variability on task learning. In the new Appendix Figure 1, we show the performance of the nontrial averaged V1 code on the orientation discrimination tasks. While the trial averaged code has much better sample efficiency, the optimal weight decay parameter λ for the noisy code is still small.
3. We also studied the effect of training readouts on subsamples of N neurons from the total V1 code. For sufficiently small N, the target function becomes unexpressible as a linear combination of the N selected neurons and the learning rule is likely to overfit. This leads to an optimal weight decay parameter λ. We mention how this is related to fitting noisy target functions [37].
4. For the behavior of nonlinear readouts we discuss the results of some recent works which allow for arbitrary nonlinear readout functions in the Discussion and Appendix Convergence of Δ Rule for Nonlinear Readouts. For classification, the power law exponents for the generalization convergence may change, but the phenonenon of codetask alignment is very qualitatively similar [2].
5. We add a Discussion section item about how the presence of neural noise impacts our metabolic efficiency argument. If the neural noise were Gaussian (perhaps with stimulus dependent covariance), then the rotations and baselines would still be legitimate and our argument still obtains. However, for many distributions, such as the Poisson distribution, the variance at a given stimulus is tied to the mean value at that stimulus so simple transformations such as shifts by a baseline or rotation (while preserving the mean response manifold) alter the noise structure and thus the information content of the code. This item is added to the Discussion section.
We will now proceed to provide more detail on each of these topics.
(1) Within this framework, a realistic stimulusdependent noise model can be easily introduced and its effects on the kernel and set of easily learned stimulusresponse mappings investigated. So while the kernel would be substantially altered via the addition of a realistic noise model, the applications of the approach outlined above would not be affected. The same cannot be said for the efficient coding application described in this manuscript. There, the authors note that rotations and constant shifts of neural activity do not affect the kernel and thus do not affect the generalization error. This kernel invariance is not present when a nontrivial (i.e. nonisotropic) noise model is added. For example, suppose that neurons are independent and Poisson so that noise scales with the mean of the neural response. In this case, adding a baseline firing rate to a population of unimodal neurons representing orientation necessarily reduces the information content of the population while rotations can affect the fidelity with which certain stimulus values are represented. It is important to note, however, that while this particular efficiency result is not compelling, I believe that it is possible to perform a similar analysis that takes into account realistic noise models and focuses on a broad set of 'biologically plausible' kernels instead of particular invariant ones. For example, one could consider noise covariance structures with differential correlations (MorenoBote 2014). Since the magnitude of differential correlations controls the redundancy of the population code this would enable an analysis of the role of redundancy in suppressing (or enhancing) generalization error.
Impact of Neural Noise on Learning Curve
We agree with the reviewer that any model of learning in a realistic neural system will have to consider the effect of noise (trialtotrial variability) on the accuracy of learned decoders/predictions. We now mention this explicitly at the beginning of the paper. At the beginning of the Results section on line 93, we added
“These responses define the population code. Throughout this work, we will mostly assume that this population code is determinstic: that identical stimuli generate identical neural responses."
The reviewer is indeed correct that an analysis of generalization can be performed for noisy population codes. Interestingly, even if the noise is stimulus dependent the generalization error can still be expressed in terms of only two correlation matrices, rather than the Fisher information at every value of the stimulus θ.
To justify this claim, first let r(θ) be stochastic with mean $\overline{r}\left(\theta \right)={\u27e8r(\theta \text{})\u27e9}_{\text{r\theta}}$ and covariance ${\mathrm{\Sigma}}_{n}\left(\theta \right)=\text{}\u27e8(r(\theta )\overline{r}(\theta ))(r(\theta )\overline{r}(\theta ){)}^{T}\u27e9$ For a given configuration of weights w, the generalization error (which is an average over the joint distribution of r,θ) is determined only by the signal ${\mathrm{\Sigma}}_{s}\text{}\left(w\right)=\text{}\u27e8\overline{r}(\theta )\overline{r}{(\theta )}^{T}\u27e9$ and (average) noise $\mathrm{\Sigma}}_{n}=\text{}\u27e8{\mathrm{\Sigma}}_{n}\text{}(\theta ){\u27e9}_{\theta$ correlation matrices:
where the (deterministic) target function is given by $y\left(\theta \right)=\text{}{w}^{\ast}\bullet \text{}\overline{r}(\theta )$. The above two expressions can be thought of as a biasvariance decomposition over the noise in neural responses. As before, the challenge is determining the average/typical case behavior of E_{g} when w is learned from the δ rule on P samples, this time with noisy responses r(θ_{µ}). The minimum achievable loss can be obtained by differentiation of the population risk expression, giving $E}_{g}^{\ast}={w}^{\ast}{\mathrm{\Sigma}}_{n}({\mathrm{\Sigma}}_{s}+{\mathrm{\Sigma}}_{n}{)}^{1}{\mathrm{\Sigma}}_{s}{w}^{\ast$.
While we have not yet computed the full learning curve for arbitrary noise correlation structure Σ_{n} (which is quite involved), we hope to pursue this calculation for future work. We want to extend our generalization analysis to arbitrary signal + arbitrary noise correlation structure, including an analysis of the role of differential correlations, but the same computation is made much more difficult by the fact that Σ_{n} and Σ_{s} are not generally codiagonalizable.
Further, we empirically tested the effect of neural variability on task learning by training readouts from the V1 code without trial averaging. We studied both the low and high frequency orientation discrimination tasks. The result is provided in R.3 and included in the paper as Appendix Figure 1. We found that the error for the trialaveraged code (dashed) is much lower than the error of the noisy code. However, the error tends to increase with weight decay parameter λ without any obvious optimal regularization choice.
We also tested subsampling the neural population, keeping only N neurons to model a sparsely connected readout. The result in Figure R.3 (CD) shows that for small N, the learning curves can exhibit overfitting peaks. Further, the generalization error is minimized for an optimal choice of the weight decay parameter.
We address and discuss these points by a new Appendix Section Impact of Neural Noise and Unlearnable Targets on Learning starting on line 1861 together with Figure R.3 (included in the paper as Appendix Figure 1). This section reads:
“While our analysis so far has focused on deterministic population codes, our theory can be extended to neural populations which exhibit variability in responses to identical stimuli. […] At infinite P, generalization error asymptotes to the amount of unexplained variance in the target function.”
We added the following to the Discussion section (line 546)
“Our work focused on how signal correlations influence inductive bias [3, 4]. However, since real neurons do exhibit variability in their responses to identical stimuli, one should consider the effect of neural noise and noise correlations in learning. We provide a preliminary analysis of learning with neural noise in Appendix Impact of Neural Noise and Unlearnable Targets on Learning, where we show that neural noise can lead to irreducible asymptotic error which depends on the geometry of the signal and noise correlations. Further, if the target function is not fully expressible as linear combinations of neural responses, overfitting peaks in the learning curves are possible, but can be mitigated with regularization implemented by a weight decay in the learning rule (see Appendix Figure 1). Future work could extend our analysis to study how signal and noise correlations interact to shape inductive bias and generalization performance in the case where the noise correlation matrices are nonisotropic, including the role of differential correlations [41]. Overall, future work could build on the present analysis to incorporate a greater degree of realism in a theory of inductive bias. ”
Impact of Neural Noise Metabolic Efficiency Results
The reviewer brings up a very good point about how neural noise can prevent decoding accuracy from being preserved under rotations and baseline shifts. To respond to this critique, we will first point out an example (stimulus dependent Gaussian noise), where generalization performance is preserved under rotations and baseline shifts. Note that if the code at θ obeyed r(θ) ∼ N(¯r(θ),Σ_{n}(θ)), then the rotated and shifted code follows Qr(θ) + δ ∼ N(Q¯r(θ) + δ,QΣ_{n}(θ)Q^{>}). This rotated and shifted code Qr(θ) + δ, when centered, will exhibit identical generalization performance as the original code. This is true both for learning from a trial averaged or nontrial averaged code. In the case of Gaussian noise on a centered code, the dataset transforms under a rotation as
D = {r_{µ},y_{µ}} → D^{0} = {Qr_{µ},y_{µ}}. The optimal weights for a linear model similarly transform as w(D) → Qw(D). Under these transformations the predictor on test point θ is unchanged since
Further, the local Fisher information matrix
is unchanged under the transformation r → Qr + δ since Σ_{n}(θ) → QΣ_{n}(θ)Q. This shows that, for some noise models, our assumption that rotations and baseline shifts preserve generalization performance will be valid.
However, we agree with the reviewer, that for Poisson noise, where the variance is tied to the mean firing rate, the RROS operations will not preserve noise structure or information content. The Fisher information at scalar stimulus θ for a Poisson neuron is $I\left(\theta \right)=\frac{{r}^{\prime}(\theta {)}^{2}}{r(\theta )}$
A baseline shift r → r + δ to the tuning curve will not change the numerator since the derivative of the tuning curve is invariant to this transformation, but it will increase the denominator.
We provide a detailed discussion of these points in Appendix Effect of Noise on RROS Symmetry (line 1772) and added a paragraph toward the end of our Discussion section, line 511, which reads:
“As a note of caution, while this analysis holds under the assumption that the neural code is deterministic, real neurons exhibit variability in their responses to repeated stimuli. Such noisy population codes do not, in general, have identical generalization performance under RROS transformations. For example, if each neuron is constrained to produce i.i.d. Poisson noise, then simple shifts of the baseline firing rate reduce the information content of the code. However, if the neural noise is Gaussian (even with stimulus dependent noise covariance), then the generalization error is conserved under RROS operations (App. Effect of Noise on RROS Symmetry). Further studies could focus on revealing the space of codes with equivalent inductive biases under realistic noise models."
(2) Similarly, the linearity assumption is somewhat restrictive. Global linear decoders of neural activity are known to be highly inefficient and completely fail when decoding orientation in the primary visual cortex in the presence of contrast fluctuations. This is because contrast modulates the amplitude of the neural response and doubling the amplitude means doubling an estimate obtained from a linear decoder even when the underlying orientation has not changed. While the contrast issue could be partially addressed by simply considering normalized neural responses, it is not yet clear how to extend this approach to account for other sources of neural variability and covariability that cause global linear decoders to fail so badly.
We appreciate the reviewer’s critique of global linear decoders. Our intention here is not to design the best decoder, but rather study a decoder’s performance set by the inductive bias of the population code. From this point of view, a linear decoder is a great choice for its analytical tractability. However, we agree that we should at least acknowledge the restrictions that come due to focusing on a linear decoder.
We discuss below some recent works which have extended our methods to nonlinear readouts such as binary classifiers which are invariant to scale transformations r → c r. Other possible nonlinear functions of w · r can also be handled within this framework under certain conditions. We have added an acknowledgement of this limitation of our present work and ways it can be extended to the discussion and the new Appendix Typical Case Analysis of Nonlinear Readouts (line 1471). The main text starting on line 541 reads:
“Though our work focused on linear readouts, arbitrary nonlinear readouts which generate convex learning objectives have been recently studied in the high dimensional limit, giving qualitatively similar learning curves which depend on kernel eigenvalues and task model alignment [1, 2] (see App. Typical Case Analysis of Nonlinear Readouts).”
Next, we provide an overview of what is discussed in Appendix Typical Case Analysis of Nonlinear Readouts. The analysis of typical case generalization can be extended to nonlinear predictors and loss functions which depend on the scalar prediction variable f(θ) = w · r(θ) [1], however the resulting equations are less interpretable. For instance, we can consider binary classification problem where y ∈ {±1} by outputting a prediction of yˆ = sign(w·r) which is invariant to simple rescaling of r → cr. We can train a model using the hinge loss `(w · r,y) = [1 − w · ry]_{+} so that the classifier will converge to a kernel support vector machine (SVM) [42]. In general, we can describe a loss function which depends on f(θ) and the ground truth function y(θ) as `(f,y).
Provided that l is convex in w · r, gradient descent will still converge to the optimum. Further, if r(θ) is approximately Gaussian, then the generalization performance can still be characterized using statistical mechanics methods [1]. While this set of results is more general, the solution for arbitrary code and task structure has some shortcomings when compared to least squares minimization including the requirement of Montecarlo sampling to compute test loss and the nondecomposability of the error into separate mode errors.
However, many qualitative features of our results continue to hold, including that the kernel’s diagonalization governs training (Appendix Typical Case Analysis of Nonlinear Readouts) and generalization and that improvements in code task alignment lead to improvements in generalization. In the followup work by Cui et al. [2], codes and tasks with power law spectra were analyzed asymptotically, showing power law generalization error decay rates E_{g} ∼ P^{−β}. The classification learning curves for power law spectra were shown to follow power laws with exponents β that are qualitatively similar to the exponents obtained with the square loss which we describe in our section titled Small and Large Sample Size Behaviors of Generalization. Just as in our theory, decay rate exponents β are larger for codes which are well aligned to the task and are smaller for codes which are nonaligned.
(3) This analysis relies on the assumption that decoder weights learned in the presence of finite data are efficient and unbiased. This assumption is problematic particularly when it comes to inductive bias and generalization error. This is because a standard way to reduce generalization error is to introduce bias into the learned decoder weights through a penalization scheme that privileges decoder weights with small magnitudes. This kind of regularization is particularly important when neurons are noisy. Fortunately, this issue could be addressed by parameterizing changes in the kernel function by the degree and type of regularization potentially leading to a more general result.
This is a great point. We now discuss inclusion of a ridge parameter λ in our calculation which privileges w with low L2 norm. The λ → 0^{+} limit recovers our previous result (see Equations 9 to 11 in Methods, and Appendix Weight Decay and Ridge Regression starting on line 1235).
We explored the effect of neural noise on the learning curve for the V1 population in orientation discrimination (see Figure R.3). The error of the noisy code (solid) is higher than the trial averaged code (dashed). For this code and tasks, the error tends to increase monotonically with regularization strength, indicating a weak effect of noise. However, in (CD), we study the learning curves for subsamples of N neurons from the entire neural population. In this case, overfitting peaks are visible at certain values of P. Similarly, there exist nonzero optimal values of weight decay λ.
Finally, I would like to conclude by explicitly stating that while the limitations imposed by the assumptions listed above temper my enthusiasm in regards to conclusions drawn in this work, I do not believe there is some fundamental problem with the general theoretical framework. Indeed, items 1 and 3 above can be easily addressed through straightforward extensions of the authors approach and I look forward to their implementation. Item 2 is a bit more troublesome, but my intuition tells me that an informationtheoretic extension based upon Fisher information may be capable of eliminating all three of these limiting assumptions by exploiting the relationship between FI(\theta) and FI(y=f(\theta)).
We are also excited by the possibility to analyze sampleefficiency of learning for neural codes with arbitrary signal and noise structure. Though we have only provided a formula for the asymptotic/irreducible error of learning with noisy neurons, we intend to analyze full learning curves in future work. As the reviewer mentions, this is a promising future direction which could even potentially benefit experimentalists in designing learning tasks which can be learned more quickly.
Ultimately, all I felt the need to say is in the public part of the review. But I wanted to use this space to insure that it was clear that I feel that this approach has a lot of potential and that I very much hope to see it extended/generalized.
If I were to make any suggestions they would be:
(1) a discussion of the consequences of biologically plausible noise models. I believe that this only requires simply augmenting K by a stimulus dependent noise covariance matrix.
(2) using this addition of a noise model to enable a discussion of the role of redundancy in generalization. This could be accomplished by considering perturbations to the noise covariance matrix that introduce differential correlations (MorenoBote 2014) of varying magnitude. More differential correlations means more redundancy and, I suspect, better generalization.
(3) the addition of a paragraph or two outlining some of the other ways that this learnability/generalization measure could be of use to physiologists and behavioral scientists. Training a mouse to do anything more complicated simple orientation detection is challenging to say the least. So having any tool that can be used to identify functions $y=f(\theta)$ that a mouse has a good chance of learning quickly would be highly advantageous. Similarly, by strategically subsampling neurons by tuning properties I suspect it may be possible to identify subpopulations of neurons that are particular important for learning certain functions. This is also cool because it allows physiologists to then target those populations for manipulation.
Of these (3) is probably the most important for this manuscript, because I don't think that many experimentalists will immediately see the sheer utility of having a tool like this in their arsenal.
In conclusion, we added the following to address these three points
1. Noise Model: We added the Appendix Section Impact of Neural Noise and Unlearnable Targets on Learning. We discussed this addition in our response above.
2. Structured Noise: We have not yet analyzed the average case learning curves for arbitrary noise correlation structure Σ_{n}(θ) (since it is much more of an involved calculation), but provide a formula for the irreducible loss. We also discuss our aspiration to pursue an analysis of differential correlations as an important future direction in the discussion.
3. Learnability for Physiologists/Behavioral Scientists: We added the possible applications to experimental neuroscience as a discussion item (see last paragraph) which we quote here:
“Finally, we discuss possible applications of our work to experimental neuroscience. Our theory has potential implications for experimental studies of task learning. First, in cases where the population selective to stimuli can be measured directly, an experimenter could design easy or difficult tasks for an animal to learn from few examples, under a hypothesis that the behavioral output is a linear function of the observed neurons. Second, in cases where it is unclear which neural population contributes to learning, one could utilize our theory to solve the inverse problem of inferring the relevant kernel from observed learning curves on different tasks [43]. From these tasks, the experimenter could compare the inferred kernel to those of different recorded populations. For instance one could compare the kernels from V1, V4, IT in visual ventral stream to the inferred kernel obtained from learning curves on certain visual learning tasks. This could provide new ways to test theories of perceptual learning [7]. Lastly, extensions of our framework could quantify the role of neural variability on task learning and the limitation it imposes on accuracy and sample efficiency.”
Reviewer #3 (Recommendations for the authors):
The manuscript presents a theory of generalization performance in deterministic population codes, that applies to the case of small numbers of training examples. The main technical result, as far as I understand, is that generalization performance (the expected classification or regression error) of a population code depends exclusively on the 'kernel', i.e. a measure of the pairwise similarity between population activity patterns corresponding to different inputs. The main conceptual results are that, using this theory, one can understand the inductive biases of the code just from analyzing the kernel, particularly the top eigenfunctions; and that sampleefficient learning (low generalization performance with few samples) depends on whether the task is aligned with the population's inductive bias, that is, whether the target function (i.e. the true map from inputs to outputs) is aligned with the top eigenfunctions of the kernel. For instance, in mouse V1 data, they show that the top eigenfunctions correspond to low frequency functions of visual orientation (i.e. functions that map a broad range of similar orientations to similar output value), and that consistent with the theory, the generalization performance for small sample sizes is better for tasks defined by low frequency target functions. In my opinion, perhaps the most significant finding from a neuroscience perspective, is that the conditions for good generalization at low samples are markedly different from those in the largesample asymptotic regime studies in Stringer et al. 2018 Nature: rather than a tradeoff between highdimensionality and differentiability proposed by Stringer et al., this manuscript shows that in the lowsample regime such codes can be disadvantageous for small sample sizes, that differentiability is not required, that the top eigenvalues matter more than the tail of the spectrum, and what matters is the alignment between the task and the top eigenfunctions. The authors propose sampleefficient learning/generalization as a new principle of neural coding, replacing or complementing efficient coding.
We thank the reviewer for such a careful reading and impressively accurate summary of our work.
Overall, in my opinion this is a remarkable manuscript, presenting truly innovative theory with somewhat limited but convincing application to neural data. My main concern is that this is highly technical, dense, and long; the mathematical proofs for the theory are buried in the supplement and require knowledge of disparate techniques from statistical physics. Although some of that material on the theory of generalization is covered in previous publications by the authors, it was not clear to me if that is true for all of the technical results or only some.
Thank you for your encouraging words! It is true that our presentation was technically dense. Here is what we did to address this:
– Throughout the paper, we tried eliminating unnecessary jargon and heavy mathematical notation, especially in the V1 model and timedependent codes sections.
– We tried providing additional intuition and simple derivations in the main text. In the problem setup, we provided a more detailed argument showing that the learned function only depends on the kernel. We also connected kernels to the idea of a distance metric in neural space. We additionally included a simple expression for the average learned function $\u27e8f(\theta ,D){\u27e9}_{D}$ in the Methods Theory of Generalization starting on line 617, showing how the coefficients for each eigenmode are learned at different rates.
– On the question of novelty, the main result on learning curve theory was derived in our prior works [44, 37], but the application of the theory to different neural codes or real neural data had not yet been performed.
Fixed population code, learnable linear readout: the authors acknowledge in the very last sentences of the manuscript that this is a limitation, given that neural tuning curves (the population neural code) are adaptable. I imagine extending the theory to both learnable codes and learnable readouts is hard and I understand it's beyond the scope of this paper. But perhaps the authors could motivate and discuss this choice, not just because of its mathematical convenience but also in relation to actual neural systems: when are these assumptions expected to be a good approximation of the real system?
This is a great point. Below we argue that our theory will apply to learning and generalization dynamics under one of the following two conditions
– Only the readout weights of a certain sensory layer (with code r(θ)) are optimized by a learning rule (such as δ rule in Figure 1). This is the simplest interpretation of our present results. It relies on an assumption that the readout weights converge more quickly than underlying low level sensory representation r(θ) evolves. We would expect this to be the case when learning a new simple task such as a Macaque learning to distinguish two new faces. The IT cortical representation will likely not need to change much in response to the classification of two new stimuli. However, if the animal is learning to classify entirely new types of stimuli over a very long period of training, it may be the case that cortical representations evolve significantly. Research in deep and recurrent networks has revealed operating regimes where the output function can still be written as f(θ) = w · ψ(θ) for a different static feature map ψ(θ) which depends on the architecture of the entire network [23, 24, 45]. For gradient flow, the kernel for this feature map is known as the neural tangent kernel [23]. Deep neural networks operating in this limit have scaling laws which are well predicted by our theory [44, 46, 47, 48]. This operating regime could be relevant for certain learning tasks in biological neural circuits.
While the above cases define the ways that our theory could accurately apply to learning in real sensory systems, there are also regimes of neural network learning where features cannot be modeled as static (in machine learning, this is known as the rich feature learning regime). Ideally a future analysis could extend our results to this setting, but it would be much more challenging. Such feature learning in a deep network is dependent on every layer of processing as well as the specific learning rule one assumes [27, 21]. Since the brain’s solution to the credit assignment problem is still unknown, it is currently unclear which learning rule should be analyzed (backpropagation, feedback alignment, global error broadcasting, gated linear networks, etc). Comparing the typical case performance and learned representations for each of these different learning rules would be a good first step towards understanding the possible inductive biases of networks in the rich regime. Provided a sufficiently refined analysis of these rules, one could perhaps even identify which learning rules are operating in real neural systems based on their measured learning dynamics.
We added the following to the discussion (starting line 520) about when we expect our theory to be accurate
“Our work constitutes a first step towards understanding inductive biases in neuronal circuits. To achieve this, we focused on a linear, δrule readout of a static population code. More work is need to study other factors that affect inductive bias. Importantly, sensory neuron tuning curves can adapt during perceptual learning tasks [7, 14, 9, 8] with the strength of adaptation dependent on brain area [10, 11, 12]. However, in many experiments, these changes to tuning in sensory areas are small [8, 9], satisfying the assumptions of our theory. For example monkeys trained on noisy visual motion detection exhibit changes in sensorymotor (LIP) but not sensory areas (MT), consistent with a model of readout from a static sensory population code [13, 15]. However, other perceptual learning tasks and other brain areas can exhibit significant changes in neural tuning [16, 17, 18]. This diversity of results motivates more general analysis of learning in multilayer networks, where the representations in each layer can adapt flexibly to task structure [19, 20, 21, 22]. Alternatively, our current analysis of inductive bias can still be consistent with multilayer learning if the network is sufficiently overparameterized and tuning curves change very little [23, 24, 19]. In this case, network training is equivalent to kernel learning with a kernel that depends on the learning rule and architecture [25]. However, in the regime of neural network training where tuning curves change significantly, more sophisticated analytical tools are needed to predict generalization [26, 27, 21]. Though our work focused on linear readouts, arbitrary nonlinear readouts which generate convex learning objectives have been recently studied in the high dimensional limit, giving qualitatively similar learning curves which depend on kernel eigenvalues and task model alignment [1, 2].”
The analysis of V1 data, showing a bias for lowfrequency functions of orientation is convincing. But it could help if the authors provided some considerations on the kind of ethological behavioral context where this is relevant, or at least the design of an experimental behavioral task to probe it. Also related, it would be useful to construct and show a counterexample, a synthetic code for which the highfrequency task is easier.
We thank the reviewer for this comment.
Applications to experimental neuroscience: We added the possible applications to experimental neuroscience as a discussion item (see last paragraph) which we quote here:
“Finally, we discuss possible applications of our work to experimental neuroscience. Our theory has potential implications for experimental studies of task learning. First, in cases where the population selective to stimuli can be measured directly, an experimenter could design easy or difficult tasks for an animal to learn from few examples, under a hypothesis that the behavioral output is a linear function of the observed neurons. Second, in cases where it is unclear which neural population contributes to learning, one could utilize our theory to solve the inverse problem of inferring the relevant kernel from observed learning curves on different tasks [43]. From these tasks, the experimenter could compare the inferred kernel to those of different recorded populations. For instance one could compare the kernels from V1, V4, IT in visual ventral stream to the inferred kernel obtained from learning curves on certain visual learning tasks. This could provide new ways to test theories of perceptual learning [7]. Lastly, extensions of our framework could quantify the role of neural variability on task learning and the limitation it imposes on accuracy and sample efficiency.”
Orientation discrimination task: We agree that this orientation discrimination task isn’t perhaps the most ethologically relevant one, but given that the data consists of responses to oriented gratings, it is a reasonable task to consider. Since neurons are selective to bars at particular orientations, this kind of task is intuitively expected to have high alignment with the V1 code. Also, because the eigenfuntions are close to Fourier modes, the inductive is more interpretable, and, thus, this task serves as a good starter example.
We also analyzed two ethologically more relevant tasks. First was a scene reconstruction task from V1 (Figure 4), and found that lowpass filtered images are easier to reconstruct, consistent with low spatial frequency selectivity of mouse V1 neurons. Second was a more ethologically natural scene categorization task (Figure 3D), but the V1 code did not generalize well, suggesting that later stages of the visual ventral stream may be better suited for categorization tasks. We suspect that other regions in visual stream may be more aligned to these ethologically relevant categorization/classification tasks. For example, primate IT codes may be more aligned with object classification.
Synthetic, counterexample code: This is a fantastic idea! Our Figure 1 actually does this. Figure 1B has two different codes, one for which a lowfrequency task is easier and another one for which a highfrequency task is easier. These codes are synthetically created by generating tuning curves from a procedure we describe in Methods Generating Example Codes.
Line 519, data preprocessing: related to the above, is it possible that binning together the V1 responses to gratings with different orientations (a range of 3.6 deg per bin, if I understood correctly) influences the finding of a lowfrequency bias?
Thanks for this comment. Based on the reviewer’s suggestion, we computed the eigendecomposition of trial averaged responses for different bin numbers bins ∈ {50,100,200} which correspond to angular windows of $\left\{{3.6}^{{}^{\circ}},{1.8}^{{}^{\circ}},{0.9}^{{}^{\circ}}\right\}$. The spectrum and eigenfunctions (especially the top modes) are very consistent across different bin sizes. The eigenfunctions become noisier with smaller bins since there are fewer trials to average over.
This figure is included as Figure 3—figure supplement 1 and referred to in the Main text line number 209.
I found the study of invariances interesting, where the theory provides a normative prediction for the proportion of simple and complex cells. However, I would suggest the authors attempt to bring this analysis a step closer to the actual data: there are no pure simple and complex cells, usually the classification is based on responses to gratings phases (F1/F0) and real neurons take a continuum of values. Could the theory qualitatively predict that distribution?
We thank the reviewer for pointing out that our model of simple/complex V1 cells was overly simplified. In reality, as the reviewer mentions there is continuous variability in F1/F0 across the neural population. In response, we modified our energy based model which produces a continuum of F1/F0 values. Further, a single scalar parameter α controls how strongly the population is biased towards invariance.
Concretely, we define the response of a cell $ri\left(\theta ,\varphi \right)=({z}_{i}\left(\theta ,\varphi \right))$ for nonlinear function g and scalar z which is constructed as follows