We denote vectors with bold lower-case symbols $\mathbf{r}$ and matrices $\mathbf{K}$ with bold upper-case symbols. We denote an average of a function $g(\mathit{\theta })$ over random variable $\mathit{\theta }$ as ${\displaystyle {⟨g(\mathit{\theta })⟩}_{\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 ${ℝ}^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 {ℝ}^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 stimulus-response 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 stimulus-response map considered (Methods Theory of generalization), and can apply to many scenarios.

The readout neuron learns from $P$ stimulus-response 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, ${\displaystyle \mathcal{D}=\{{\mathit{\theta }}^{\mu },y({\mathit{\theta }}^{\mu }){\}}_{\mu =1}^P}$, the readout weights are learned with the local, biologically-plausible delta-rule, $\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 ${\displaystyle {\mathbf{w}}^{\ast }(\mathcal{D})}$ (Appendix Convergence of the delta-rule 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 ${\displaystyle \mathcal{D}}$. Our theoretical analysis of this model provides insights into how populations of neurons encode information and allow sample-efficient learning.

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 ${\displaystyle {\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 ${\displaystyle {\mathbf{w}}^{\ast }(\mathcal{D})}$ for the training set ${\displaystyle \mathcal{D}}$ (Neal, 1994; Girosi et al., 1995). The dataset-dependent fixed point ${\displaystyle {\mathbf{w}}^{\ast }(\mathcal{D})}$ of the learning rule is a linear combination of the population vectors on the dataset ${\displaystyle {\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 delta-rule without weight decay), and the matrix $\mathbf{K}$ has entries $K_{\mu \nu }=K({\mathit{\theta }}^{\mu },{\mathit{\theta }}^{\nu })$ and ${\displaystyle y_{\mu }=y({\mathit{\theta }}^{\mu })}$. The matrix K⁺ is the pseudo-inverse 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 Delta-rule 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 $\tilde{\mathbf{r}}(\mathit{\theta })=\mathbf{Qr}(\mathit{\theta })$ with an identical kernel (Appendix Alternative neural codes with same kernel) since $\frac{1}{N}\tilde{\mathbf{r}}(\mathit{\theta })\cdot \tilde{\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 $\tilde{\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 $\tilde{\mathbf{r}}(\mathit{\theta })$ (Figure 2A) which have the same kernels (Figure 2B) and the same performance on a learning task (Figure 2C).

Figure 2

Download asset Open asset

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, ${\displaystyle E_g={⟨E_g(\mathcal{D})⟩}_{\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 left-singular 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; Litwin-Kumar 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 inner-product 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).

Figure 3 with 1 supplement see all

Download asset Open asset

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 ${\displaystyle 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 code-task alignment). We formalize this notion by the following metric. Mathematically, generalization error ${\displaystyle ⟨E_g⟩}$ can be decomposed into normalized estimation errors $E_k$ for the coefficients of these eigenfunctions ${\psi }_k$, ${\displaystyle {⟨E_g⟩}_{\mathcal{D}}=\sum _k{v}_k^2{E}_k}$, where ${\displaystyle E_k={⟨({\hat{v}}_k(\mathcal{D})-v_k{)}^2⟩}_{\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 code-task alignment), (Bordelon et al., 2020): ${\displaystyle {\lambda }_k>{\lambda }_{\ell }⟹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 code-task alignment as a principle for good generalization. To quantify code-task 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):

Stimulus-response 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 code-task alignment).

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 square-wave 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, high-frequency 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 low-frequency task, the two target values are well separated along the top two axes. However, the high-frequency 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 ill-suited 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 band-pass filtered images with spatial frequency wave-vector $\mathbf{k}\in {ℝ}^2$ constrained to an annulus ${\displaystyle |\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 ${\displaystyle r<s_{max}}$, which quantifies the dimension of the function space. We see that the lower frequency band-limited 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.

Figure 4 with 1 supplement see all

Download asset Open asset

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 power-law 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.

Figure 5 with 3 supplements see all

Download asset Open asset

Next, to further illustrate the importance of code-task 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 })+(1-s)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 phase-selective 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.

Recently, Stringer et al., 2018b argued that the input-output 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. Non-differentiable 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 non-differentiable kernels generalize better than differentiable kernels (Figure 6—figure supplement 2).

Figure 6 with 2 supplements see all

Download asset Open asset

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 delta-rule 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 well-performing 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 low-sample 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 non-differentiable 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 non-differentiable kernels generalizing better than differentiable kernels in Figure 6—figure supplement 2.

Our framework can directly be extended to learning of arbitrary time-varying functions of time-varying 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 ${\displaystyle \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).

Figure 7

Download asset Open asset

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 time-series ${\psi }_k(\{\mathit{\theta }\};t)$ for each eigenvalue ${\lambda }_k$. An interesting link can be made with this eigensystem and linear low-dimensional 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 ${\displaystyle \mathcal{P}}$.

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 non-negative. 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, $\tilde{\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.

Figure 8 with 2 supplements see all

Download asset Open asset

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 off-set 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 non-negative codes with full-rank kernels will have selective tuning.

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, ${\displaystyle {\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 {ℝ}^{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=\cos (\theta )-0.6\cos (4\theta )$, while the second experiment used $y=\cos (6\theta )-\cos (8\theta )$.

Recent work has established analytic results that predict the average case generalization error for kernel regression

where ${\displaystyle E_g(\mathcal{D})={⟨[f(\mathit{\theta },\mathcal{D})-y(\mathit{\theta }){]}^2⟩}_{\mathit{\theta }}}$ is the generalization error for a certain sample ${\displaystyle \mathcal{D}}$ of size $P$ and ${\displaystyle f(\mathit{\theta },\mathcal{D})=\mathbf{w}\cdot \mathbf{r}(\mathit{\theta })}$ is the kernel regression solution for ${\displaystyle \mathcal{D}}$ (Appendix Convergence of the delta-rule 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 ${\displaystyle f(\mathit{\theta },\mathcal{D})}$ in this basis ${\displaystyle y(\mathit{\theta })=\sum _k{v}_k{\psi }_k(\mathit{\theta }),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 ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle {\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 ${\displaystyle 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 ${\displaystyle [\mathbf{k}(\mathit{\theta }){]}_{\mu }=K(\mathit{\theta },{\mathit{\theta }}^{\mu })}$ and ${\displaystyle [\mathbf{y}{]}_{\mu }=y({\mathit{\theta }}^{\mu })}$ for each training stimulus ${\displaystyle {\mathit{\theta }}^{\mu }\in \mathcal{D}}$. The $P\times P$ kernel gram matrix $\mathbf{K}$ has entries ${\displaystyle [\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 ${\displaystyle \mathcal{D}}$ only appears through the matrix ${\displaystyle \mathbf{G}(\mathcal{D})}$. Thus to compute the typical generalization error we need to average this matrix over realizations of datasets, i.e.${\displaystyle {⟨\mathbf{G}(\mathcal{D})⟩}_{\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 ${\displaystyle {⟨f(\mathit{\theta },\mathcal{D})⟩}_{\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 ${\displaystyle \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 ${\displaystyle 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.

For the simulations in Figure 7, we integrated a rate-based 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(1-t){[\cos ({\gamma }^{\mu }),\sin ({\gamma }^{\mu })]}^{\top }\in {ℝ}^2$. Target function for a delay $d$ is $\mathbf{y}({\theta }^{\mu },t)=H(1.5+d-t)H(t-d-1){[\cos ({\gamma }^{\mu }),\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 {ℝ}^{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 ${\displaystyle \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.

Requests for information should be directed to the lead contact, Cengiz Pehlevan (cpehlevan@seas.harvard.edu).

Mouse V1 neuron responses to orientation gratings and preprocessing code were obtained from a publicly available dataset: https://github.com/MouseLand/stringer-et-al-2019, (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/stringer-pachitariu-et-al-2018b (Stringer et al., 2018a).

The code generated by the authors for this paper is also available https://github.com/Pehlevan-Group/sample_efficient_pop_codes (Pehlevan-Group, 2022).

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 non-negative (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₂ 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 {ℝ}^{N\times N}$ of expansion coefficients. Note that ${\displaystyle \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 non-degenerate. 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$.

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}^{\tilde{P}}{p}_i\delta (\mathit{\theta }-{\mathit{\theta }}^i)$, where $i$ indexes all the $\tilde{P}$ values $\mathit{\theta }$ can take. With this choice

Demanding this equality for ${\displaystyle {\mathit{\theta }}^{\mathrm{\prime }}={\mathit{\theta }}^i,i=1,...,\tilde{P}}$ generates a matrix eigenvalue problem

where ${\displaystyle {\mathbf{B}}_{ij}={\delta }_{ij}{p}_i}$. The eigenfunctions over the stimuli are identified as the columns of ${\displaystyle \mathbf{\Psi }}$ while the eigenvalues are the diagonal elements of ${\displaystyle {\mathbf{\Lambda }}_{k\ell }={\lambda }_k{\delta }_{k\ell }}$.

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 ${\displaystyle {\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 ${\displaystyle {\mathcal{S}}_{\mathbf{k}}={\mathbb{Z}}^n=\{\mathbf{k}\in {\mathbb{R}}^n|k_i\in \mathbb{Z}\mathrm{\forall }i=1,...,n\}}$. Reality and symmetry of the kernel demand that $\text{Im}({\lambda }_{\mathbf{k}})=0$ and ${\displaystyle {\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 $\cos (a-b)=\cos (a)\cos (b)+\sin (a)\sin (b)$. By recognizing that ${\{\sqrt{2}\cos (k\theta ),\sqrt{2}\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.

Suppose the kernel $K(\mathit{\theta },{\mathit{\theta }}^{\prime })$ is invariant to some set of transformations ${\displaystyle 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 ${\displaystyle {\lambda }_k>0}$, then

This establishes that all functions which depend on ${\displaystyle \mathcal{T}}$ transformations must necessarily lie in the null-space of $K$.