#  The Onset of Variance-Limited Behavior for Networks in the Lazy and Rich
Regimes

Alexander Atanasov §‡ ,  Blake Bordelon∗ †‡ ,  Sabarish Sainathan †‡ & Cengiz
Pehlevan †‡
§Department of Physics
†John A. Paulson School of Engineering and Applied Sciences
‡Center for Brain Science
Harvard University
Cambridge, MA 02138, USA
<EMAIL_ADDRESS>
These authors contributed equally.

###### Abstract

For small training set sizes $P$, the generalization error of wide neural
networks is well-approximated by the error of an infinite width neural network
(NN), either in the kernel or mean-field/feature-learning regime. However,
after a critical sample size $P^{*}$, we empirically find the finite-width
network generalization becomes worse than that of the infinite width network.
In this work, we empirically study the transition from infinite-width behavior
to this variance-limited regime as a function of sample size $P$ and network
width $N$. We find that finite-size effects can become relevant for very small
dataset sizes on the order of $P^{*}\sim\sqrt{N}$ for polynomial regression
with ReLU networks. We discuss the source of these effects using an argument
based on the variance of the NN’s final neural tangent kernel (NTK). This
transition can be pushed to larger $P$ by enhancing feature learning or by
ensemble averaging the networks. We find that the learning curve for
regression with the final NTK is an accurate approximation of the NN learning
curve. Using this, we provide a toy model which also exhibits
$P^{*}\sim\sqrt{N}$ scaling and has $P$-dependent benefits from feature
learning.

##  1 Introduction

Deep learning systems are achieving state of the art performance on a variety
of tasks (Tan & Le, 2019; Hoffmann et al., 2022). Exactly how their
generalization is controlled by network architecture, training procedure, and
task structure is still not fully understood. One promising direction for deep
learning theory in recent years is the infinite-width limit. Under a certain
parameterization, infinite-width networks yield a kernel method known as the
neural tangent kernel (NTK) (Jacot et al., 2018; Lee et al., 2019). Kernel
methods are easier to analyze, allowing for accurate prediction of the
generalization performance of wide networks in this regime (Bordelon et al.,
2020; Canatar et al., 2021; Bahri et al., 2021; Simon et al., 2021). Infinite-
width networks can also operate in the mean-field regime if network outputs
are rescaled by a small parameter $\alpha$ that enhances feature learning (Mei
et al., 2018; Chizat et al., 2019; Geiger et al., 2020b; Yang & Hu, 2020;
Bordelon & Pehlevan, 2022).

While infinite-width networks provide useful limiting cases for deep learning
theory, real networks have finite width. Analysis at finite width is more
difficult, since predictions are dependent on the initialization of
parameters. While several works have attempted to analyze feature evolution
and kernel statistics at large but finite width (Dyer & Gur-Ari, 2020; Roberts
et al., 2021), the implications of finite width on generalization are not
entirely clear. Specifically, it is unknown at what value of the training set
size $P$ the effects of finite width become relevant, what impact this
critical $P$ has on the learning curve, and how it is affected by feature
learning.

To identify the effects of finite width and feature learning on the deviation
from infinite width learning curves, we empirically study neural networks
trained across a wide range of output scales $\alpha$, widths $N$, and
training set sizes $P$ on the simple task of polynomial regression with a ReLU
neural network. Concretely, our experiments show the following:

  * •

Learning curves for polynomial regression transition exhibit significant
finite-width effects very early, around $P\sim\sqrt{N}$. Finite-width NNs at
large $\alpha$ are always outperformed by their infinite-width counterparts.
We show this gap is driven primarily by variance of the predictor over
initializations (Geiger et al., 2020a). Following prior work (Bahri et al.,
2021), we refer to this as the variance-limited regime. We compare three
distinct ensembling methods to reduce error in this regime.

  * •

Feature-learning NNs show improved generalization both before and after the
transition to the variance limited regime. Feature learning can be enhanced
through re-scaling the output of the network by a small scalar $\alpha$ or by
training on a more complex task (a higher-degree polynomial). We show that
alignment between the final NTK and the target function on test data improves
with feature learning and sample size.

  * •

We demonstrate that the learning curve for the NN is well-captured by the
learning curve for kernel regression with the final empirical NTK, eNTKf , as
has been observed in other works (Vyas et al., 2022; Geiger et al., 2020b;
Atanasov et al., 2021; Wei et al., 2022).

  * •

Using this correspondence between the NN and the final NTK, we provide a
cursory account of how fluctuations in the final NTK over random
initializations are suppressed at large width $N$ and large feature learning
strength. In a toy model, we reproduce several scaling phenomena, including
the $P\sim\sqrt{N}$ transition and the improvements due to feature learning
through an alignment effect.

We validate that these effects qualitatively persist in the realistic setting
of wide ResNets Zagoruyko & Komodakis (2017) trained on CIFAR in appendix E.

Overall, our results indicate that the onset of finite-width corrections to
generalization in neural networks become relevant when the scale of the
variance of kernel fluctuations becomes comparable to the bias component of
the generalization error in the bias-variance decomposition. The variance
contribution to generalization error can be reduced both through ensemble
averaging and through feature learning, which we show promotes higher
alignment between the final kernel and the task. We construct a model of noisy
random features which reproduces the essential aspects of our observations.

###  1.1 Related Works

Geiger et al. (2020a) analyzed the scaling of network generalization with the
number of model parameters. Since the NTK fluctuates with variance $O(N^{-1})$
for a width $N$ network (Dyer & Gur-Ari, 2020; Roberts et al., 2021), they
find that finite width networks in the lazy regime generically perform worse
than their infinite width counterparts.

The scaling laws of networks over varying $N$ and $P$ were also studied, both
empirically and theoretically by Bahri et al. (2021). They consider two types
of learning curve scalings. First, they describe resolution-limited scaling,
where either training set size or width are effectively infinite and the
scaling behavior of generalization error with the other quantity is studied.
There, the scaling laws can been obtained by the theory in Bordelon et al.
(2020). Second, they analyze variance-limited scaling where width or training
set size are fixed to a finite value and the other parameter is taken to
infinity. While that work showed for any fixed $P$ that the learning curve
converges to the infinite width curve as $O(N^{-1})$, these asymptotics do not
predict, for fixed $N$, at which value of $P$ the NN learning curve begins to
deviate from the infinite width theory. This is the focus of our work.

The contrast between rich and lazy networks has been empirically studied in
several prior works. Depending on the structure of the task, the lazy regime
can have either worse (Fort et al., 2020) or better (Ortiz-Jiménez et al.,
2021; Geiger et al., 2020b) performance than the feature learning regime. For
our setting, where the signal depends on only a small number of relevant input
directions, we expect representation learning to be useful, as discussed in
(Ghorbani et al., 2020; Paccolat et al., 2021b). Consequently, we posit and
verify that the rich network will outperform the lazy one.

Our toy model is inspired by the literature on random feature models. Analysis
of generalization for two layer networks at initialization in the limit of
high dimensional data have been carried out using techniques from random
matrix theory (Mei & Montanari, 2022; Hu & Lu, 2020; Adlam & Pennington,
2020a; Dhifallah & Lu, 2020; Adlam & Pennington, 2020b) and statistical
mechanics (Gerace et al., 2020; d’Ascoli et al., 2020; 2020). Several of these
works have identified that when $N$ is comparable to $P$, the network
generalization error has a contribution from variance over initial parameters.
Further, they provide a theoretical explanation of the benefit of ensembling
predictions of many networks trained with different initial parameters.
Recently, Ba et al. (2022) studied regression with the hidden features of a
two layer network after taking one step of gradient descent, finding
significant improvements to the learning curve due to feature learning.
Zavatone-Veth et al. (2022) analyzed linear regression for Bayesian deep
linear networks with width $N$ comparable to sample size $P$ and demonstrated
the advantage of training multiple layers compared to only training the only
last layer, finding that feature learning advantage has leading correction of
scale $(P/N)^{2}$ at small $P/N$.

##  2 Problem Setup and Notation

We consider a supervised task with a dataset
$\mathcal{D}=\\{\bm{x}^{\mu},y^{\mu}\\}_{\mu=1}^{P}$ of size $P$. The pairs of
data points are drawn from a population distribution $p(\bm{x},y)$. Our
experiments will focus on training networks to interpolate degree $k$
polynomials on the sphere (full details in Appendix A). For this task, the
infinite width network learning curves can be found analytically. In
particular at large $P$ the generalization error scales as $1/P^{2}$ (Bordelon
et al., 2020). We take a single output feed-forward NN
$\tilde{f}_{\theta}:\mathbb{R}^{D}\to\mathbb{R}$ with hidden width $N$ for
each layer. We let $\theta$ denote all trainable parameters of the network.
Using NTK parameterization (Jacot et al., 2018), the activations for an input
$\bm{x}$ are given by

$\displaystyle
h^{(\ell)}_{i}=\frac{\sigma}{\sqrt{N}}\sum_{j=1}^{N}W_{ij}^{(\ell)}\varphi(h_{j}^{(\ell-1)}),\quad\ell=2,\dots
L,\quad h^{(1)}_{i}=\frac{\sigma}{\sqrt{D}}\sum_{j=1}^{D}W^{(1)}_{ij}x_{j}.$
(1)

Here, the output of the network is $\tilde{f}_{\theta}=h_{1}^{(L)}$. We take
$\varphi$ to be a positively homogenous function, in our case a ReLU
nonlinearity, but this is not strictly necessary (Appendix C.2). At
initialization we have $W_{ij}\sim\mathcal{N}(0,1)$. Consequently, the scale
of the output at initialization is $O(\sigma^{L})$. As a consequence of the
positive homogeneity of the network, the scale of the output is given by
$\alpha=\sigma^{L}$. $\alpha$ controls the feature learning strength of a
given NN. Large $\alpha$ corresponds to a lazy network while small $\alpha$
yields a rich network with feature movement. More details on how $\alpha$
controls feature learning are given in Appendix C.1 and C.2.

In what follows, we will denote the infinite width NTK limit of this network
by NTK∞. We will denote its finite width linearization by
$\mathrm{eNTK}_{0}(\bm{x},\bm{x}^{\prime}):=\sum_{\theta}\partial_{\theta}f(\bm{x})\partial_{\theta}f(\bm{x}^{\prime})|_{\theta=\theta_{0}}$,
and we will denote its linearization around its final parameters $\theta_{f}$
by
$\mathrm{eNTK}_{f}(\bm{x},\bm{x}^{\prime}):=\sum_{\theta}\partial_{\theta}f(\bm{x})\partial_{\theta}f(\bm{x}^{\prime})|_{\theta=\theta_{f}}$.
Following other authors (Chizat et al., 2019; Adlam & Pennington, 2020a), we
will take the output to be
$f_{\theta}(\bm{x}):=\tilde{f}_{\theta}(\bm{x})-\tilde{f}_{\theta_{0}}(\bm{x})$.
Thus, at initialization the function output is $0$. We explain this choice
further in Appendix A. The parameters are then trained with full-batch
gradient descent on a mean squared error loss. We denote the final network
function starting from initialization $\theta_{0}$ on a dataset $\mathcal{D}$
by $f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})$ or $f^{*}$ for short. The
generalization error is calculated using a held-out test set and approximates
the population risk $E_{g}(f):=\left<(f(\bm{x})-y)^{2}\right>_{\bm{x},y\sim
p(\bm{x},y)}$.

##  3 Empirical Results

In this section, we will study learning curves for ReLU NNs trained on
polynomial regression tasks of varying degrees. We take our task to be
learning $y=Q_{k}(\bm{\beta}\cdot\bm{x})$ where $\bm{\beta}$ is random vector
of norm $1/D$ and $Q_{k}$ is the $k$th gegenbauer polynomial. We will
establish the following key observations, which we will set out to
theoretically explain in Section 4.

  1. 1.

Both eNTK0 and sufficiently lazy networks perform strictly worse than NTK∞ ,
but the ensembled predictors approach the NTK∞ test error.

  2. 2.

NNs in the feature learning regime of small $\alpha$ can outperform NTK∞ for
an intermediate range of $P$. Over this range, the effect of ensembling is
less notable.

  3. 3.

Even richly trained finite width NNs eventually perform worse than NTK∞ at
sufficiently large $P$. However, these small $\alpha$ feature-learning
networks become variance-limited at larger $P$ than lazy networks. Once in the
variance-limited regime, all networks benefit from ensembling over
initializations.

  4. 4.

For all networks, the transition to the variance-limited regime begins at a
$P^{*}$ that scales sub-linearly with $N$. For polynomial regression, we find
$P^{*}\sim\sqrt{N}$.

These findings support our hypothesis that finite width introduces variance in
eNTK0 over initializations, which ultimately leads to variance in the learned
predictor and higher generalization error. Although we primarily focus on
polynomial interpolation tasks in this paper, in Appendix F we provide results
for wide ResNets trained on CIFAR and observe that rich networks also
outperform lazy ones, and that lazy ones benefit more significantly from
ensembling.

###  3.1 Finite width effects cause the onset of a variance limited regime

(a) $k=2$ generalization error (b) $k=2$ $20$-fold ensemble error (c) $k=4$
generalization error (d) $k=4$ $20$-fold ensemble error

Figure 1: Generalization errors of depth $L=3$ neural networks across a range
of $\alpha$ values compared to NTK∞ . The regression for NTK∞ was calculated
using the Neural Tangents package (Novak et al., 2020). The exact scaling of
NTK∞ is known to go asymptotically as $P^{-2}$ for this task. a) Lazy networks
perform strictly worse than NTK∞ while rich networks can outperform it for an
intermediate range of $P$ before their performance is also limited. b)
Ensembling 20 networks substantially improves lazy network and eNTK0
generalization, as well as asymptotic rich network generalization. This
indicates that at sufficiently large $P$, these neural networks become limited
by variance due to initialization. The error bars in a) and c) denote the
variance due to both both training set and initialization. The error bars in
b), d) denote the variance due to the train set.

In this section, we first investigate how finite width NN learning curves
differ from infinite width NTK regression. In Figure 1 we show the
generalization error $E_{g}(f^{*}_{\theta_{0},\mathcal{D}})$ for a depth 3
network with width $N=1000$ trained on a quadratic $k=2$ and quartic $k=4$
polynomial regression task. Additional plots for other degree polynomials are
provided in Appendix F. We sweep over $P$ to show the effect of more data on
generalization, which is the main relationship we are interested in studying.
For each training set size we sweep over a grid of 20 random draws of the
train set and 20 random network initializations. This for 400 trained networks
in total at each choice of $P,k,N,\alpha$. We see that a discrepancy arises at
large enough $P$ where the neural networks begin to perform worse than NTK∞ .

We probe the source of the discrepancy between finite width NNs and NTK∞ by
ensemble averaging network predictions $\bar{f}_{\mathcal{D}}(\bm{x}):=\langle
f_{\theta_{0},\mathcal{D}}^{*}(\bm{x})\rangle_{\theta_{0}}$ over $E=20$
initializations $\theta_{0}$. In Figures 1b and 1d, we calculate the error of
$\bar{f}_{\mathcal{D}}(\bm{x})$, each trained on the same dataset. We then
plot $E_{g}(\bar{f}_{\mathcal{D}})$. This ensembled error approximates the
bias in a bias-variance decomposition (Appendix B). Thus, any gap between 1
(a) and 1 (b) is driven by variance of $f_{\theta,\mathcal{D}}$ over $\theta$.

We sharpen these observations with phase plots of NN generalization, variance
and kernel alignment over $P,\alpha$, as shown in Figure 2. In Figure 2a,
generalization for NNs in the rich regime (small $\alpha$) have lower final
$E_{g}$ than lazy networks. As the dataset grows, the fraction of $E_{g}$ due
to initialization variance (that is, the fraction removed by ensembling)
strictly increases (2 (b)). We will show why this effect occurs in section
3.2. Figure 2b shows that, at any fixed $P$, the variance is lower for small
$\alpha$. To measure the impact of feature learning on the eNTKf , we plot its
alignment with the target function, measured as
$\frac{{\bm{y}}^{\top}{\bm{K}}{\bm{y}}}{|{\bm{y}}|^{2}\text{Tr}{\bm{K}}}$ for
a test set of targets $[\bm{y}]_{\mu}$ and kernel
$[\bm{K}]_{\mu\nu}=\mathrm{eNTK}_{f}({\bm{x}}_{\mu},{\bm{x}}_{\nu})$.
Alignment of the kernel with the target function is known to be related to
good generalization (Canatar et al., 2021). In Section 4, we revisit these
effects in a simple model which relates kernel alignment and variance
reduction.

(a) $k=3$ generalization error (b) $k=3$ variance fraction (c) $k=3$ alignment

Figure 2: Phase plots in the $P,\alpha$ plane of a) The log generalization
error $\log_{10}E_{g}(f^{\star})$, b) The fraction of generalization error
removed by ensembling $1-E_{g}(\bar{f}^{\star})/E_{g}(f^{\star})$, c) Kernel-
task alignment measured by
$\frac{\bm{y}^{T}K_{f}\bm{y}}{|\bm{y}|^{2}\mathrm{Tr}K_{f}}$ where $\bm{y}$
and $K_{f}$ are evaluated on test data. We have plotted ‘x’ markers in a) to
show the points where the NNs were trained.

In addition to initialization variance, variance over dataset $\mathcal{D}$
contributes to the total generalization error. Following (Adlam & Pennington,
2020b), we discuss a symmetric decomposition of the variance in Appendix B,
showing the contribution from dataset variance and the effects of bagging. We
find that most of the variance in our experiments is due to initialization.

We show several other plots of the results of these studies in the appendix.
We show the effect of bagging (Figure 7), phase plots of different degree
target functions (Figures 10, 9), phase plots over $N,\alpha$ (Figure 11) and
a comparison of network predictions against the initial and final kernel
regressors (Figures 18, 19).

###  3.2 Final NTK Variance leads to Generalization Plateau

In this section, we show how the variance over initialization can be
interpreted as kernel variance in both the rich and lazy regimes. We also show
how this implies a plateau for the generalization error.

To begin, we demonstrate empirically that all networks have the same
generalization error as kernel regression solutions with their final eNTKs. At
large $\alpha$, the initial and the final kernel are already close, so this
follows from earlier results of Chizat et al. (2019). In the rich regime, the
properties of the eNTKf have been studied in several prior works. Several have
empirically demonstrated that the eNTKf is a good match to the final network
predictor for a trained network (Long, 2021; Vyas et al., 2022; Wei et al.,
2022) while others have given conditions under which such an effect would hold
true (Atanasov et al., 2021; Bordelon & Pehlevan, 2022). We comment on this in
appendix C.4. We show in Figure 3 how the final network generalization error
matches the generalization error of eNTKf . As a consequence, we can use eNTKf
to study the observed generalization behavior.

(a) $E_{g}^{NN}=E_{g}^{NTK_{f}}$ (b) $E_{g}^{NN}=E_{g}^{NTK_{f}}$ across
$N,\alpha$

Figure 3: Kernel regression with eNTKf reproduces the learning curves of the
NN with high fidelity. (a) Learning curves across different laziness settings
$\alpha$ in a width $1000$ network. The solid black curve is the infinite
width network. Colored curves are the NN generalizations. Stars represent the
eNTKf s, and lie on top of the corresponding NN learning curves. (b) The
agreement of generalizations between NNs and eNTKf s across different $N$ and
$\alpha$. Here the colors denote different $\alpha$ values while the dot,
triangle and star markers denote networks of $N=\\{177,421,1000\\}$
respectively.

Next, we relate the variance of the final predictor
$f_{\theta_{0},\mathcal{D}}^{*}$ to the corresponding infinite width network
$f^{\infty}_{\mathcal{D}}$. The finite size fluctuations of the kernel at
initialization have been studied in (Dyer & Gur-Ari, 2020; Hanin & Nica, 2019;
Roberts et al., 2021). The variance of the kernel elements has been shown to
scale as $1/N$. We perform the following bias-variance decomposition: Take
$f_{\theta_{0},\mathcal{D}}$ to be the eNTK0 predictor, or a sufficiently lazy
network trained to interpolation on a dataset $\mathcal{D}$. Then,

$\displaystyle\langle(f_{\theta_{0},\mathcal{D}}^{*}(\bm{x})-y)^{2}\rangle_{\theta_{0},\mathcal{D},\bm{x},y}$
$\displaystyle=\langle(f^{\infty}_{\mathcal{D}}(\bm{x})-y)^{2}\rangle_{\mathcal{D},\bm{x},y}+O(1/N).$
(2)

We demonstrate this equality using a relationship between the infinite-width
network and an infinite ensemble of finite-width networks derived in Appendix
B. There we also show that the $O(1/N)$ term is strictly positive for
sufficiently large $N$. Thus, for lazy networks of sufficiently large $N$,
finite width effects lead to strictly worse generalization error. The
decomposition in Equation 2 continues to hold for rich networks at small
$\alpha$ if $f^{\infty}$ is interpreted as the infinite-width mean field
limit. In this case one can show that ensembles of rich networks are
approximating an infinite width limit in the mean-field regime. See Appendix B
for details.

###  3.3 Feature Learning delays variance limited transition

(a) Scaling of $P_{1/2}$ with $\alpha$ (b) Scaling of $P_{1/2}$ with $N$

Figure 4: Critical sample size $P_{1/2}$ measures the onset of the variance
limited regime as a function of $\alpha$ at fixed $N$. (a) More feature
learning (small $\alpha$) delays the transition to the variance limited
regime. (b) $P_{1/2}$ as a function of $N$ for fixed $\alpha$ has roughly
$P_{1/2}\sim\sqrt{N}$ scaling.

We now consider how feature learning alters the onset of the variance limited
regime, and how this onset scales with $\alpha,N$. We define the onset of the
variance limited regime to take place at the value $P^{*}=P_{1/2}$ where over
half of the generalization error is due to variance over initializations.
Equivalently we have $E_{g}(\bar{f}^{*})/E_{g}(f^{*})=1/2$. By using an
interpolation method together with bisection, we solve for $P_{1/2}$ and plot
it in Figure 4.

Figure 4b shows that $P_{1/2}$ scales as $\sqrt{N}$ for this task. In the next
section, we shall show that this scaling is governed by the fact that
$P_{1/2}$ is close to the value where the infinite width network
generalization curve $E_{g}^{\infty}$ is equal to the variance of eNTKf . In
this case the quantities to compare are $E_{g}^{\infty}\approx P^{-2}$ and
$\mathrm{Var}\,\mathrm{eNTK}_{f}\approx N^{-1}$.

We can understand the delay of the variance limited transition, as well as the
lower value of the final plateau using a mechanistic picture similar to the
effect observed in Atanasov et al. (2021). In that setting, under small
initialization, the kernel follows a deterministic trajectory, picking up a
low rank component in the direction of the train set targets
$\bm{y}\bm{y}^{\top}$, and then changing only in scale as the network weights
grow to interpolate the dataset. In their case, for initial output scale
$\sigma^{L}$, eNTKf is deterministic up to a variance of $O(\sigma)$. In our
case, the kernel variance at initialization scales as $\sigma^{2L}/N$. As
$\sigma\to 0$ the kernel’s trajectory becomes deterministic up to a variance
term scaling with $\sigma$ as $O(\sigma)$, which implies that the final
predictor also has a variance scaling as $O(\sigma)$.

##  4 Signal plus noise correlated feature model

(a) Ensembling Methods (b) Reduction in $E_{g}$ for Each Ensembling Technique

Figure 5: The random feature model suggests three possible types of
ensembling: averaging the output function $f({\bm{x}},\theta)$, averaging
eNTKf $K({\bm{x}},{\bm{x}}^{\prime};\theta)$, and averaging the induced
features $\bm{\psi}({\bm{x}},\theta)$. We analyze these ensembling methods for
a $k=1$ task with a width $N=100$ ReLU network. (a) While all ensembling
methods improve generalization, averaging either the kernel $\left<K\right>$
or features $\left<\psi\right>$ gives a better improvement to generalization
than averaging the output function $\left<f\right>$. Computing final kernels
for many richly trained networks and performing regression with this averaged
kernel gives the best performance. (b) We plot the relative error of each
ensembling method against the single init neural network. The gap between
ensembling and the single init NN becomes evident for sufficiently large
$P\sim P_{1/2}$. For small $\alpha$, all ensembling methods perform
comparably, while for large $\alpha$ ensembling the kernel or features gives
much lower $E_{g}$ than averaging the predictors.

In Section 3.2 we have shown that in both the rich and lazy regimes, the
generalization error of the NN is well approximated by the generalization of a
kernel regression solution with eNTKf . This finding motivates an analysis of
the generalization of kernel machines which depend on network initialization
$\theta_{0}$. Unlike many analyses of random feature models which specialize
to two layer networks and focus on high dimensional Gaussian random data (Mei
& Montanari, 2022; Adlam & Pennington, 2020a; Gerace et al., 2020; Ba et al.,
2022), we propose to analyze regression with the eNTKf for more general
feature structures. This work builds on the kernel generalization theory for
kernels developed with statistical mechanics (Bordelon et al., 2020; Canatar
et al., 2021; Simon et al., 2021; Loureiro et al., 2021). We will attempt to
derive approximate learning curves in terms of the eNTKf ’s signal and noise
components, which provide some phenomenological explanations of the onset of
the variance limited regime and the benefits of feature learning. Starting
with the final NTK $K_{\theta_{0}}({\bm{x}},{\bm{x}}^{\prime})$ which depends
on the random initial parameters $\theta_{0}$, we project its square root
$K^{1/2}_{\theta_{0}}({\bm{x}},{\bm{x}}^{\prime})$ (as defined in equation 32)
on a fixed basis $\\{b_{k}(x)\\}_{k=1}^{\infty}$ orthonormal with respect to
$p({\bm{x}})$. This defines a feature map

$\displaystyle\psi_{k}({\bm{x}},\theta_{0})=\int
d{\bm{x}}^{\prime}p({\bm{x}}^{\prime})K^{1/2}_{\theta_{0}}({\bm{x}},{\bm{x}}^{\prime})b_{k}({\bm{x}}^{\prime})\
,\ k\in\\{1,...,\infty\\}.$ (3)

The kernel can be reconstructed from these features
$K_{\theta_{0}}({\bm{x}},{\bm{x}}^{\prime})=\sum_{k}\psi_{k}({\bm{x}},\theta_{0})\psi_{k}({\bm{x}}^{\prime},\theta_{0})$.
The kernel interpolation problem can be solved by performing linear regression
with features $\bm{\psi}({\bm{x}},\theta_{0})$. Here,
${\bm{w}}(\theta_{0})=\lim_{\lambda\to
0}\text{argmin}_{{\bm{w}}}\sum_{\mu=1}^{P}[{\bm{w}}\cdot\bm{\psi}({\bm{x}}_{\mu},\theta_{0})-y_{\mu}]^{2}+\lambda|{\bm{w}}|^{2}$.
The learned function
$f({\bm{x}},\theta_{0})={\bm{w}}(\theta_{0})\cdot\bm{\psi}({\bm{x}},\theta_{0})$
is the minimum norm interpolator for the kernel
$K({\bm{x}},{\bm{x}}^{\prime};\theta_{0})$ and matches the neural network
learning curve as seen in Section 3.2. In general, since the rank of $K$ is
finite for a finite size network, the $\psi_{k}({\bm{x}},\bm{\theta}_{0})$
have correlation matrix of finite rank $N_{\mathcal{H}}$. Since the target
function $y$ does not depend on the initialization $\theta_{0}$, we decompose
it in terms of a fixed set of features
$\bm{\psi}_{M}({\bm{x}})\in\mathbb{R}^{M}$ (for example, the first $M$ basis
functions $\\{b_{k}\\}_{k=1}^{M}$). In this random feature model, one can
interpret the initialization-dependent fluctuations in
$K({\bm{x}},{\bm{x}}^{\prime};\theta_{0})$ as generating fluctuations in the
features $\bm{\psi}({\bm{x}},\theta_{0})$ which induce fluctuations in the
learned network predictor $f({\bm{x}},\theta_{0})$. To illustrate the relative
improvements to generalization from denoising these three different objects,
in Figure 5, we compare averaging the final kernel $K$, averaging the induced
features $\psi$, and averaging network predictions $f$ directly. For all
$\alpha$, all ensembling methods provide improvements over training a single
NN. However, we find that averaging the kernel directly and performing
regression with this kernel exhibits the largest reduction in generalization
error. Averaging features performs comparably. However, ensemble averaging
network predictors does not perform as well as either of these other two
methods. The gap between ensembling methods is more significant in the lazy
regime (large $\alpha$) and is negligible in the rich regime (small $\alpha$).

###  4.1 Toy Models and Approximate Learning Curves

To gain insight into the role of feature noise, we characterize the test error
associated with a Gaussian covariate model in a high dimensional limit
$P,M,N_{\mathcal{H}}\to\infty$ with $\alpha=P/M,\eta=N_{\mathcal{H}}/M$.

$\displaystyle
y=\frac{1}{\sqrt{M}}\bm{\psi}_{M}\cdot{\bm{w}}^{*},f=\frac{1}{\sqrt{M}}\bm{\psi}\cdot{\bm{w}},\
\bm{\psi}={\bm{A}}(\bm{\theta}_{0})\bm{\psi}_{M}+\bm{\epsilon},\
\begin{bmatrix}\bm{\psi}_{M}\\\
\bm{\epsilon}\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}\bm{\Sigma}_{M}&0\\\
0&\bm{\Sigma}_{\epsilon}\end{bmatrix}\right)$ (4)

This model was also studied by Loureiro et al. (2021) and subsumes the classic
two layer random feature models of prior works (Hu & Lu, 2020; Adlam &
Pennington, 2020a; Mei & Montanari, 2022). The expected generalization error
for any distribution of ${\bm{A}}(\bm{\theta}_{0})$ has the form

$\displaystyle\mathbb{E}_{\bm{\theta}_{0}}E_{g}(\bm{\theta}_{0})=\mathbb{E}_{{\bm{A}}}\frac{1}{1-\gamma}\frac{1}{M}{\bm{w}}^{*}\bm{\Sigma}_{M}^{1/2}\left[{\bm{I}}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}^{2}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}\right]\bm{\Sigma}_{M}^{1/2}{\bm{w}}^{*}$
$\displaystyle{\bm{G}}=\left({\bm{I}}+\hat{q}{\bm{A}}\bm{\Sigma}_{M}{\bm{A}}^{\top}+\hat{q}\bm{\Sigma}_{\epsilon}\right)^{-1}\
,\ \hat{q}=\frac{\alpha}{\lambda+q}\ ,\
q=\text{Tr}{\bm{G}}[{\bm{A}}\bm{\Sigma}_{M}{\bm{A}}^{\top}+\bm{\Sigma}_{\epsilon}],$
(5)

where $\alpha=P/M$ and
$\gamma=\frac{\alpha}{(\lambda+q)^{2}}\text{Tr}{\bm{G}}^{2}[{\bm{A}}\bm{\Sigma}_{M}{\bm{A}}^{\top}+\bm{\Sigma}_{\epsilon}]^{2}$.
Details of the calculation can be found in Appendix D. We also provide
experiments showing the predictive accuracy of the theory in Figure 6. In
general, we do not know the induced distribution of ${\bm{A}}(\theta_{0})$
over disorder $\theta_{0}$. In Appendix D.5, we compute explicit learning
curves for a simple toy model where ${\bm{A}}(\bm{\theta}_{0})^{\prime}s$
entries as i.i.d. Gaussian over the random initialization $\bm{\theta}_{0}$. A
similar random feature model was recently analyzed with diagrammatic
techniques by Maloney et al. (2022). In the high dimensional limit
$M,P,N_{\mathcal{H}}\to\infty$ with $P/M=\alpha,N_{\mathcal{H}}/M=\eta$, our
replica calculation demonstrates that test error is self-averaging (the same
for every random instance of ${\bm{A}}$) which we describe in Appendix D.5 and
Figure 16.

###  4.2 Explaining Feature Learning Benefits and Error Plateaus

(a) $E_{g}^{NN}$ for different $N$ (b) Small $N$ $\approx$ Large $\sigma^{2}$
(c) Variance Limited Transition (d) Feature Scalings (e) Richness $\approx$
Amplified $\bm{\Sigma}_{M}$ (f) $P$-dependent Amplification

Figure 6:  A toy model of noisy features reproduces qualitative dependence of
learning curves on kernel fluctuations and feature learning. (a) The empirical
learning curves for networks of varying width $N$ at large $\alpha$. (b) Noisy
kernel regression learning curve with noise
$\bm{\Sigma}_{\epsilon}=\sigma_{\epsilon}^{2}\bm{\Sigma}_{M}$ and ${\bm{A}}$
is a projection matrix preserving $20$-k top eigenmodes of $\bm{\Sigma}_{M}$,
which was computed from the NTK∞ for a depth $3$ ReLU network. (c) This toy
model reproduces the approximate scaling of the transition sample size
$P_{1/2}\sim N^{1/2}$ if $\sigma^{2}_{\epsilon}\sim N^{-1}$. (d) NNs trained
with varying richness $\alpha$. Small $\alpha$ improves the early learning
curve and asymptotic behavior. (e) Theory curves for a kernel with amplified
eigenvalue $\lambda_{k}\to\lambda_{k}+\Delta\lambda_{k}$ for the target
eigenfunction. This amplification mimics the effect of enhanced kernel
alignment in the low $\alpha$ regime. Large amplification improves
generalization performance. (f) $P$-dependent alignment where
$\Delta\lambda_{k}\sim\sqrt{P}$ gives a better qualitative match to (d).

Using this theory, we can attempt to explain some of the observed phenomena
associated with the onset of the variance limited regime. First, we note that
the kernels exhibit fluctuations over initialization with variance $O(1/N)$,
either in the lazy or rich regime. In Figure 6 (a), we show learning curves
for networks of different widths in the lazy regime. Small width networks
enter the variance limited regime earlier and have higher error. Similarly, if
we alter the scale of the noise
$\bm{\Sigma}_{\epsilon}=\sigma^{2}_{\epsilon}\bm{\Sigma}_{M}$ in our toy
model, the corresponding transition time $P_{1/2}$ is smaller and the
asymptotic error is higher. In Figure 6 (c), we show that our theory also
predicts the onset of the variance limited regime at $P_{1/2}\sim\sqrt{N}$ if
$\sigma^{2}_{\epsilon}\sim N^{-1}$. We stress that this scaling is a
consequence of the structure of the task. Since the target function is an
eigenfunction of the kernel, the infinite width error goes as $1/P^{2}$
(Bordelon et al., 2020). Since variance scales as $1/N$, bias and variance
become comparable at $P\sim\sqrt{N}$. Often, realistic tasks exhibit power law
decays where $E_{g}^{N=\infty}=P^{-\beta}$ with $\beta<2$ (Spigler et al.,
2020; Bahri et al., 2021), where we’d expect a transition around $P_{1/2}\sim
N^{1/\beta}$.

Using our model, we can also approximate the role of feature learning as
enhancement in the signal correlation along task-relevant eigenfunctions. In
Figure 6 (d) we plot the learning curves for networks trained with different
levels of feature learning, controlled by $\alpha$. We see that feature
learning leads to improvements in the learning curve both before and after
onset of variance limits. In Figure 6 (e)-(f), we plot the theoretical
generalization for kernels with enhanced signal eigenvalue for the task
eigenfunction $y({\bm{x}})=\phi_{k}({\bm{x}})$. This enhancement, based on the
intuition of kernel alignment, leads to lower bias and lower asymptotic
variance. However, this model does not capture the fact that feature learning
advantages are small at small $P$ and that the slopes of the learning curves
are different at different $\alpha$. Following the observation of Paccolat et
al. (2021a) that kernel alignment can occur with scale $\sqrt{P}$, we plot the
learning curves for signal enhancements that scale as $\sqrt{P}$. Though this
toy model reproduces the onset of the variance limited regime $P_{1/2}$ and
the reduction in variance due to feature learning, our current result is not
the complete story. A more refined future theory could use the structure of
neural architecture to constrain the structure of the ${\bm{A}}$ distribution.

##  5 Conclusion

We performed an extensive empirical study for deep ReLU NNs learning a fairly
simple polynomial regression problems. For sufficiently large dataset size
$P$, all neural networks under-perform the infinite width limit, and we
demonstrated that this worse performance is driven by initialization variance.
We show that the onset of the variance limited regime can occur early in the
learning curve with $P_{1/2}\sim\sqrt{N}$, but this can be delayed by
enhancing feature learning. Finally, we studied a simple random-feature model
to attempt to explain these effects and qualitatively reproduce the observed
behavior, as well as quantitatively reproducing the relevant scaling
relationship for $P_{1/2}$. This work takes a step towards understanding
scaling laws in regimes where finite-size networks undergo feature learning.
This has implications for how the choice of initialization scale, neural
architecture, and number networks in an ensemble can be tuned to achieve
optimal performance under a fixed compute and data budget.

## References

  * Adlam & Pennington (2020a) Ben Adlam and Jeffrey Pennington.  The neural tangent kernel in high dimensions: Triple descent and a multi-scale theory of generalization.  In _International Conference on Machine Learning_ , pp. 74–84. PMLR, 2020a.
  * Adlam & Pennington (2020b) Ben Adlam and Jeffrey Pennington.  Understanding double descent requires a fine-grained bias-variance decomposition.  _Advances in neural information processing systems_ , 33:11022–11032, 2020b.
  * Atanasov et al. (2021) Alexander Atanasov, Blake Bordelon, and Cengiz Pehlevan.  Neural networks as kernel learners: The silent alignment effect.  In _International Conference on Learning Representations_ , 2021.
  * Ba et al. (2022) Jimmy Ba, Murat A Erdogdu, Taiji Suzuki, Zhichao Wang, Denny Wu, and Greg Yang.  High-dimensional asymptotics of feature learning: How one gradient step improves the representation.  _arXiv preprint arXiv:2205.01445_ , 2022.
  * Bahri et al. (2021) Yasaman Bahri, Ethan Dyer, Jared Kaplan, Jaehoon Lee, and Utkarsh Sharma.  Explaining neural scaling laws, 2021.  URL https://arxiv.org/abs/2102.06701.
  * Baratin et al. (2021) Aristide Baratin, Thomas George, César Laurent, R Devon Hjelm, Guillaume Lajoie, Pascal Vincent, and Simon Lacoste-Julien.  Implicit regularization via neural feature alignment.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 2269–2277. PMLR, 2021.
  * Bordelon & Pehlevan (2022) Blake Bordelon and Cengiz Pehlevan.  Self-consistent dynamical field theory of kernel evolution in wide neural networks.  _arXiv preprint arXiv:2205.09653_ , 2022.
  * Bordelon et al. (2020) Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan.  Spectrum dependent learning curves in kernel regression and wide neural networks.  In Hal Daumé III and Aarti Singh (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 1024–1034. PMLR, 13–18 Jul 2020.  URL https://proceedings.mlr.press/v119/bordelon20a.html.
  * Bradbury et al. (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, et al.  Jax: composable transformations of python+ numpy programs.  _Version 0.2_ , 5:14–24, 2018.
  * Canatar et al. (2021) Abdulkadir Canatar, Blake Bordelon, and Cengiz Pehlevan.  Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks.  _Nature Communications_ , 12, 2021.
  * Chizat et al. (2019) Lénaïc Chizat, Edouard Oyallon, and Francis R. Bach.  On lazy training in differentiable programming.  In _NeurIPS_ , 2019.
  * Cortes et al. (2012) Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh.  Algorithms for learning kernels based on centered alignment.  _The Journal of Machine Learning Research_ , 13(1):795–828, 2012.
  * d’Ascoli et al. (2020) Stéphane d’Ascoli, Levent Sagun, and Giulio Biroli.  Triple descent and the two kinds of overfitting: Where & why do they appear?  _Advances in Neural Information Processing Systems_ , 33:3058–3069, 2020.
  * Dhifallah & Lu (2020) Oussama Dhifallah and Yue M Lu.  A precise performance analysis of learning with random features.  _arXiv preprint arXiv:2008.11904_ , 2020.
  * Dyer & Gur-Ari (2020) Ethan Dyer and Guy Gur-Ari.  Asymptotics of wide networks from feynman diagrams.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=S1gFvANKDS.
  * d’Ascoli et al. (2020) Stéphane d’Ascoli, Maria Refinetti, Giulio Biroli, and Florent Krzakala.  Double trouble in double descent: Bias and variance (s) in the lazy regime.  In _International Conference on Machine Learning_ , pp. 2280–2290. PMLR, 2020.
  * Fort et al. (2020) Stanislav Fort, Gintare Karolina Dziugaite, Mansheej Paul, Sepideh Kharaghani, Daniel M Roy, and Surya Ganguli.  Deep learning versus kernel learning: an empirical study of loss landscape geometry and the time evolution of the neural tangent kernel.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 5850–5861. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/405075699f065e43581f27d67bb68478-Paper.pdf.
  * Geiger et al. (2020a) Mario Geiger, Arthur Jacot, Stefano Spigler, Franck Gabriel, Levent Sagun, Stéphane d’Ascoli, Giulio Biroli, Clément Hongler, and Matthieu Wyart.  Scaling description of generalization with number of parameters in deep learning.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(2):023401, 2020a.
  * Geiger et al. (2020b) Mario Geiger, Stefano Spigler, Arthur Jacot, and Matthieu Wyart.  Disentangling feature and lazy training in deep neural networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(11):113301, 2020b.
  * Gerace et al. (2020) Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Generalisation error in learning with random features and the hidden manifold model.  In _International Conference on Machine Learning_ , pp. 3452–3462. PMLR, 2020.
  * Ghorbani et al. (2020) Behrooz Ghorbani, Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  When do neural networks outperform kernel methods?  In _NeurIPS_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/a9df2255ad642b923d95503b9a7958d8-Abstract.html.
  * Hanin & Nica (2019) Boris Hanin and Mihai Nica.  Finite depth and width corrections to the neural tangent kernel.  In _International Conference on Learning Representations_ , 2019.
  * Hoffmann et al. (2022) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_ , 2022.
  * Hu & Lu (2020) Hong Hu and Yue M Lu.  Universality laws for high-dimensional learning with random features.  _arXiv preprint arXiv:2009.07669_ , 2020.
  * Jacot et al. (2018) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: convergence and generalization in neural networks (invited paper).  _Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing_ , 2018.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kornblith et al. (2019) Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey Hinton.  Similarity of neural network representations revisited.  In _International Conference on Machine Learning_ , pp. 3519–3529. PMLR, 2019.
  * Lee et al. (2019) Jaehoon Lee, Lechao Xiao, Samuel S. Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jascha Sohl-Dickstein.  Wide neural networks of any depth evolve as linear models under gradient descent.  _ArXiv_ , abs/1902.06720, 2019.
  * Long (2021) Philip M. Long.  Properties of the after kernel.  _CoRR_ , abs/2105.10585, 2021.  URL https://arxiv.org/abs/2105.10585.
  * Loureiro et al. (2021) Bruno Loureiro, Cédric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Capturing the learning curves of generic features maps for realistic data sets with a teacher-student model.  _CoRR_ , abs/2102.08127, 2021.  URL https://arxiv.org/abs/2102.08127.
  * Maloney et al. (2022) Alexander Maloney, Daniel A. Roberts, and James Sully.  A solvable model of neural scaling laws, 2022.  URL https://arxiv.org/abs/2210.16859.
  * Mei & Montanari (2022) Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and the double descent curve.  _Communications on Pure and Applied Mathematics_ , 75(4):667–766, 2022.
  * Mei et al. (2018) Song Mei, Andrea Montanari, and Phan-Minh Nguyen.  A mean field view of the landscape of two-layer neural networks.  _Proceedings of the National Academy of Sciences_ , 115(33):E7665–E7671, 2018.
  * Novak et al. (2020) Roman Novak, Lechao Xiao, Jiri Hron, Jaehoon Lee, Alexander A. Alemi, Jascha Sohl-Dickstein, and Samuel S. Schoenholz.  Neural tangents: Fast and easy infinite neural networks in python.  In _International Conference on Learning Representations_ , 2020.  URL https://github.com/google/neural-tangents.
  * Ortiz-Jiménez et al. (2021) Guillermo Ortiz-Jiménez, Seyed-Mohsen Moosavi-Dezfooli, and Pascal Frossard.  What can linearized neural networks actually say about generalization?  _Advances in Neural Information Processing Systems_ , 34:8998–9010, 2021.
  * Paccolat et al. (2021a) Jonas Paccolat, Leonardo Petrini, Mario Geiger, Kevin Tyloo, and Matthieu Wyart.  Geometric compression of invariant manifolds in neural networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2021(4):044001, apr 2021a.  doi: 10.1088/1742-5468/abf1f3.  URL https://doi.org/10.1088/1742-5468/abf1f3.
  * Paccolat et al. (2021b) Jonas Paccolat, Leonardo Petrini, Mario Geiger, Kevin Tyloo, and Matthieu Wyart.  Geometric compression of invariant manifolds in neural networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2021(4):044001, 2021b.
  * Roberts et al. (2021) Daniel A. Roberts, Sho Yaida, and Boris Hanin.  The principles of deep learning theory, 2021.
  * Simon et al. (2021) James B. Simon, Madeline Dickens, and Michael R. DeWeese.  Neural tangent kernel eigenvalues accurately predict generalization, 2021\.
  * Spigler et al. (2020) Stefano Spigler, Mario Geiger, and Matthieu Wyart.  Asymptotic learning curves of kernel methods: empirical data versus teacher–student paradigm.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(12):124001, 2020.
  * Tan & Le (2019) Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In _International conference on machine learning_ , pp. 6105–6114. PMLR, 2019.
  * Vyas et al. (2022) Nikhil Vyas, Yamini Bansal, and Preetum Nakkiran.  Limitations of the ntk for understanding generalization in deep learning.  _arXiv preprint arXiv:2206.10012_ , 2022.
  * Wei et al. (2022) Alexander Wei, Wei Hu, and Jacob Steinhardt.  More than a toy: Random matrix models predict how real-world neural representations generalize.  _arXiv preprint arXiv:2203.06176_ , 2022.
  * Yang & Hu (2020) Greg Yang and Edward J. Hu.  Feature learning in infinite-width neural networks.  _ArXiv_ , abs/2011.14522, 2020.
  * Zagoruyko & Komodakis (2017) Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks, 2017.
  * Zavatone-Veth et al. (2022) Jacob A Zavatone-Veth, William L Tong, and Cengiz Pehlevan.  Contrasting random and learned features in deep bayesian linear regression.  _arXiv preprint arXiv:2203.00573_ , 2022.

##  Appendix A Details on Experiments

We generated the dataset $\mathcal{D}=\\{\bm{x}^{\mu},y^{\mu}\\}_{\mu=1}^{P}$
by sampling $\bm{x}^{\mu}$ uniformly on $\mathbb{S}^{D-1}$, the unit sphere in
$\mathbb{R}^{D}$. $\tilde{y}$ was then generated as a Gegenbauer polynomial of
degree $k$ of a 1D projection of $\bm{x}$,
$\tilde{y}=Q_{k}(\bm{\beta}\cdot\bm{x})$. Because the scale of the output of
the neural network relative to the target is a central quantity in this work,
it is especially important to make sure the target is appropriately scaled to
unit norm. We did this by defining the target to be $y=\tilde{y}/\sqrt{\langle
Q_{k}(\bm{\beta}\cdot\bm{x})^{2}\rangle_{\bm{x}\sim\mathbb{S}^{D-1}}}$. The
denominator can be easily and accurately approximated by Monte Carlo sampling.

We used JAX (Bradbury et al., 2018) for all neural network training. We built
multi-layer perceptrons (MLPs) of depth 2 and 3. Most of the results are
reported for depth 3 perceptrons, where there is a separation between the
width of the network $N$ and the number of parameters $N^{2}$. Sweeping over
more depths and architectures is possible, but because of the extensive
dimensionality of the hyperparameter search space, we have not yet
experimented with deeper networks.

We considered MLPs with no bias terms. Since the Gegenbauer polynomials are
mean zero, we do not need biases to fit the training set and generalize well.
We have also verified that adding trainable biases does not change the final
results in any substantial way.

As mentioned in the main text, we consider the final output function to be the
initial network output minus the output at initialization:

$f_{\theta}(\bm{x})=\tilde{f}_{\theta}(\bm{x})-\tilde{f}_{\theta_{0}}(\bm{x}).$
(6)

Here, only $\theta$ is differentiated through, while $\theta_{0}$ is held
fixed. The rationale for this choice is that without this subtraction, in the
lazy limit the trained neural network output can be written as

$\tilde{f}_{\theta}^{*}(\bm{x})=\tilde{f}_{\theta_{0}}(\bm{x})+\sum_{\mu\nu}\bm{k}_{\mu}(\bm{x})[\bm{K}^{-1}]_{\mu\nu}(y^{\nu}-\tilde{f}_{\theta_{0}}(\bm{x})).$
(7)

This is the same as doing eNTK0 regression on the shifted targets
$y^{\mu}-\tilde{f}_{\theta_{0}}(\bm{x})$. At large initialization the shift
$\tilde{f}_{\theta_{0}}(\bm{x})$ amounts to adding random, initialization-
dependent noise to the targets. By instead performing the subtraction, the
lazy limit can be interpreted as a kernel regression on the targets
themselves, which is preferable.

We trained this network with full batch gradient descent with a learning rate
$\eta$ so that

$\displaystyle\Delta\theta$
$\displaystyle=-\eta\nabla_{\theta}\mathcal{L}(\mathcal{D},\theta),$ (8)
$\displaystyle\mathcal{L}(\mathcal{D},\theta)$
$\displaystyle:=\frac{1}{P}\sum_{\mu=1}^{P}|f_{\theta}(\bm{x}^{\mu})-y^{\mu}|^{2}.$

Each network was trained to an interpolation threshold of $10^{-6}$. If a
network could not reach this threshold in under 30k steps, we checked if the
training error was less than $10$ times the generalization error. If this was
not satisfied, then that run of the network was discarded.

For each fixed $P,k$, we generated 20 independent datasets. For each fixed
$N,\alpha$ we generated 20 independent neural network initializations. This
$20\times 20$ table yields a total of 400 neural networks trained on every
combination of initialization and dataset choice.

The infinite width network predictions were calculated using the Neural
Tangents package (Novak et al., 2020). The finite width eNTK0 s were also
calculated using the empirical methods in Neural Tangents. They were trained
to interpolation using the gradient_descent_mse method. This is substantially
faster than training the linearized model using standard full-batch gradient
descent, which we have found to take a very long time for most networks. We
use the same strategy for the eNTKf s.

For the experiments in the main text, we have taken the input dimension to be
$D=10$ and sweep over $k=1,2,3,4$. We swept over 15 values $P$ in logspace
from size $30$ to size 10k, and over 6 values of $N$ in logspace from size
$30$ to size $2150$. We then swept over alpha values $0.1,0.5,1.0,10.0,20.0$.
Depending on $\alpha,N$, we tuned the learning rate $\eta$ of the network
small enough to stay close to the gradient flow limit, but allow for the
interpolation threshold to be feasibly reached.

For each of the 1800 settings of $P,N,\alpha,k$ and each of the 400 networks,
400 eNTK0 s, 400 eNTKf s, and 20 NTK∞ s, the generalization error was
saved, as well as a vector of $\hat{y}$ predictions on a test set of 2000
points. In addition, for the neural networks we saved both initial and final
parameters. All are saved as lists of numpy arrays in a directory of about
1TB. We plan to make the results of our experiments publicly accessible,
alongside the code to generate them.

###  A.1 CIFAR Experiments

We apply the same methodology of centering the network and allowing $\alpha$
to control the degree of laziness by redefining

$f_{\theta}(\bm{x})=\alpha(\tilde{f}_{\theta}(\bm{x})-\tilde{f}_{\theta_{0}}(\bm{x})).$
(9)

We consider the task of binary classification for CIFAR-10. In order to allow
$P$ to become large we divide the data into two classes: animate and inanimate
objects. We choose to subsample eight classes and superclass them into two:
(cat, deer, dog, horse) vs (airplane, automobile, ship, truck). Each
superclass consists of 20,000 training examples and 4,000 test examples
retrieved from the CIFAR-10 dataset.

On subsets of this dataset, we train wide residual networks (ResNets)
Zagoruyko & Komodakis (2017) of width $64$ and block size $1$ with the NTK
parameterization Jacot et al. (2018) on this task using mini-batch gradient
descent with batch size of 256 and MSE loss. Step sizes are governed by the
Adam optimizer Kingma & Ba (2014) with initial learning rate
$\eta_{0}=10^{-3}.$ Every network is trained for 24,000 steps, such that under
nearly all settings of $\alpha$ and dataset size the network has attained
infinitesimal train loss.

We sweep $\alpha$ from $10^{-3}$ to $10^{0}$ and $P$ from $2^{9}$ to $2^{15}$.
For each value of $P$, we randomly sample five training datasets of size $P$
and compute ensembles of size 20. For each network in an ensemble the
initialization and the order of the training data is randomly chosen
independently of those for the other networks.

##  Appendix B Fine-grained bias-variance decomposition

###  B.1 Fine Grained Decomposition of Generalization Error

Let $\mathcal{D}$ be a dataset of $(\bm{x}^{\mu},y^{\mu})_{\mu=1}^{P}\sim
p(\bm{x},y)$ viewed as a random variable. Let $\theta_{0}$ represent the
initial parameters of a neural network, viewed as a random variable. In the
case of no label noise, as in section 2.2.1 of Adlam & Pennington (2020b), we
derive the symmetric decomposition of the generalization error in terms of the
variance due to initialization and the variance due to the dataset. We have

$\displaystyle E_{g}(f^{*}_{\theta_{0},\mathcal{D}})$
$\displaystyle=\langle(f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})-y)^{2}\rangle_{\bm{x},y}=\langle(\langle
f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})\rangle_{\theta_{0},\mathcal{D}}-y)^{2}\rangle_{\bm{x},y}+\mathbb{E}_{\bm{x}}\mathrm{Var}_{\theta_{0},\mathcal{D}}f^{*}_{\theta_{0},\mathcal{D}}(y)$
(10)
$\displaystyle=\mathrm{Bias}^{2}+V_{\mathcal{D}}+V_{\theta_{0}}+V_{\mathcal{D},\theta_{0}}.$

Here we have defined

$\displaystyle\mathrm{Bias}^{2}$ $\displaystyle=\langle(\langle
f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})\rangle_{\theta_{0},\mathcal{D}}-y)^{2}\rangle_{\bm{x},y},$
(11) $\displaystyle V_{\mathcal{D}}$
$\displaystyle=\mathbb{E}_{\bm{x}}\,\mathrm{Var}_{\mathcal{D}}\,\mathbb{E}_{\theta_{0}}[f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})|\mathcal{D}]=\mathbb{E}_{\bm{x}}\,\mathrm{Var}_{\mathcal{D}}\,\bar{f}_{\mathcal{D}}^{*}(\bm{x}),$
(12) $\displaystyle V_{\theta_{0}}$
$\displaystyle=\mathbb{E}_{\bm{x}}\,\mathrm{Var}_{\theta_{0}}\,\mathbb{E}_{\mathcal{D}}[f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})|\theta_{0}],$
(13) $\displaystyle V_{\mathcal{D},\theta_{0}}$
$\displaystyle=\mathbb{E}_{\bm{x}}\mathrm{Var}_{\theta_{0},\mathcal{D}}f^{*}_{\theta_{0},\mathcal{D}}(y)-V_{\theta_{0}}-V_{\mathcal{D}}.$
(14)

$V_{\mathcal{D}}$ and $V_{\theta_{0}}$ give the components of the variance
explained by variance in $\mathcal{D},\theta_{0}$ respectively.
$V_{\mathcal{D},\theta_{0}}$ is the remaining part of the variance not
explained by either of these two sources. As in the main text,
$\bar{f}^{*}_{\mathcal{D}}(\bm{x})$ is the ensemble average of the trained
predictors over initializations.
$E_{\mathcal{D}}[f^{*}_{\theta_{0},\mathcal{D}}(\bm{x})|\theta_{0}]$ is
commonly referred to as the bagged predictor. In the next subsection we study
these terms empirically.

###  B.2 Empirical Study of Dataset Variance

Using the network simulations, one can show that the bagged predictor does not
have substantially lower generalization error in the regimes that we are
interested in. This implies that most of the variance driving higher
generalization error is due to variance over initializations. In figure 7, we
make phase plots of the fraction of $E_{g}$ that arises from variance due to
initialization, variance over datasets, and total variance for width 1000.
This can be obtained by computing the ensembled predictor, the bagged
predictor, and the ensembled-bagged predictor respectively.

(a) Initialization variance $k=2$ (b) Initialization variance $k=3$ (c)
Initialization variance $k=4$ (d) Dataset variance $k=2$ (e) Dataset variance
$k=3$ (f) Dataset variance $k=4$ (g) Total variance $k=2$ (h) Total variance
$k=3$ (i) Total variance $k=4$

Figure 7: Phase plots of the fraction of the generalization error due to the
initialization variance, the dataset variance, and their combined
contribution. The columns correspond to the tasks of polynomial regression for
degree $2,3$ and $4$ polynomials. Neural network has width 1000 and depth 3.
Notice that the initialization variance dominates in the large $P$ large
$\alpha$ regime

###  B.3 Relating Ensembled Network Generalization to Infinite Width
Generalization

Making use of the fact that at leading order, the eNTKf (either in the rich or
lazy regime) of a trained network has $\theta_{0}$-dependent fluctuations with
variance $1/N$, one can write the kernel Gram matrices as

$\displaystyle\,[\bm{K}_{\theta_{0}}]_{\mu\nu}$
$\displaystyle=[\bm{K}_{\infty}]_{\mu\nu}+\frac{1}{\sqrt{N}}[\bm{\delta
K}_{\theta_{0}}]_{\mu\nu}+O(1/N)$ (15)
$\displaystyle[\bm{k}_{\theta_{0}}(\bm{x})]_{\mu}$
$\displaystyle=[\bm{k}_{\infty}(\bm{x})]_{\mu}+\frac{1}{\sqrt{N}}[\bm{\delta
k}_{\theta_{0}}(\bm{x})]_{\mu}+O(1/N).$

Here, $\bm{\delta K}_{\theta_{0}},\bm{\delta k}_{\theta_{0}}$ are the leading
order fluctuations around the infinite width network. Because of how we have
written them, their variance is $O(1)$ with respect to $N$. Using perturbation
theory (Dyer & Gur-Ari, 2020), one can demonstrate that these leading order
terms have mean zero around their infinite-width limit.

The predictor for the eNTK0 (or for a sufficiently large $\alpha$ neural
network) for a training set with target labels $\bm{y}$ is given by:

$\displaystyle f^{*}(\bm{x})_{\theta_{0}}$
$\displaystyle=\bm{k}_{\theta_{0}}(\bm{x})_{\mu}^{\top}\bm{K}_{\theta_{0}}^{-1}\cdot\bm{y}$
(16) $\displaystyle=f^{\infty}(\bm{x})+\frac{1}{\sqrt{N}}\bm{\delta
k}_{\theta_{0}}(\bm{x})^{\top}\bm{K}_{\infty}^{-1}\cdot\bm{y}-\frac{1}{\sqrt{N}}\bm{k}_{\infty}(\bm{x})^{\top}\bm{K}_{\infty}^{-1}\bm{\delta}K_{\theta_{0}}\bm{K}_{\infty}^{-1}\cdot\bm{y}+O\left(N^{-1}\right).$

This implies (Geiger et al., 2020a):

$\langle(f^{*}_{\theta_{0}}(\bm{x})-f^{\infty}(\bm{x}))^{2}\rangle_{\bm{x}}=O\left(N^{-1}\right).$
(17)

Upon taking the ensemble, because of the mean zero property of the deviations,
we get that

$\displaystyle\langle f^{*}(\bm{x})_{\theta_{0}}\rangle_{\theta_{0}}$
$\displaystyle=f^{\infty}(\bm{x})+O\left(N^{-1}\right)$ (18)
$\displaystyle\Rightarrow\langle(\langle
f^{*}(\bm{x})_{\theta_{0}}\rangle_{\theta_{0}}-f^{\infty}(\bm{x}))^{2}\rangle=O\left(N^{-2}\right).$

We can now bound the generalization error of the ensemble of networks in terms
of the infinite-width generalization:

$\displaystyle\langle(\langle
f^{*}_{\theta_{0}}(\bm{x})\rangle_{\theta_{0}}-y)^{2}\rangle_{\bm{x},y}$
$\displaystyle=\langle(f^{\infty}(\bm{x})-y)^{2}\rangle_{\bm{x},y}+\langle(\langle
f^{*}_{\theta_{0}}(\bm{x})\rangle_{\theta_{0}}-f^{\infty}(\bm{x}))^{2}\rangle_{\bm{x}}$
(19)
$\displaystyle\quad-2\langle(f^{\infty}(\bm{x})-y)(f^{\infty}(\bm{x})-\langle
f^{*}_{\theta_{0}}(\bm{x})\rangle_{\theta_{0}})\rangle_{\bm{x},y}.$

By equation 18, the second term yields a positive contribution going as
$O(N^{-2})$. The last term can be bounded by Cauchy-Schwarz:

$\displaystyle|\langle(f^{\infty}(\bm{x})-y)(f^{\infty}(\bm{x})-\langle
f^{*}_{\theta_{0}}(\bm{x})\rangle_{\theta_{0}})\rangle_{\bm{x},y}|$
$\displaystyle\leq\sqrt{\langle(f^{\infty}(\bm{x})-y)^{2}\rangle\langle(f^{\infty}(\bm{x})-\langle
f^{*}_{\theta_{0}}(\bm{x})\rangle_{\theta_{0}})^{2}\rangle_{\bm{x},y}}$ (20)
$\displaystyle=\sqrt{E_{g}^{\infty}(P)c_{1}}/N.$

After we enter the variance limited regime by taking $P>P_{1/2}$ we get
$E_{g}^{\infty}\leq O(1/N)$ so this last term is bounded by $N^{-3/2}$.
Consequently, the difference in generalization error between the infinite
width NTK and an ensemble of lazy network or eNTK0 predictors is subleading in
$1/N$ compared to the generalization gap, which goes as $N^{-1}$.

The same argument can be extended to any predictor that differs from some
infinite width limit. In particular Bordelon & Pehlevan (2022) show that the
fluctuations of the eNTKf in any mean field network are asymptotically mean
zero with variance $N^{-1}$. The above argument then applies to the predictor
obtained by ensembling networks that have learned features. This implies that
in the variance limited regime, ensemble averages of feature learning networks
have the same generalization as the infinite-width mean field solutions up to
a term that decays faster than $N^{-3/2}$.

##  Appendix C Feature Learning

###  C.1 Controlling Feature Learning Through Initialization Scale

Given the feed-forward network defined in equation 1, one can see that the
components of the activations satisfy $h_{i}^{(\ell)}=O(\sigma
h_{i})^{(\ell-1)}$ and consequently that the output
$h^{(L)}_{1}=O(\sigma^{L})$. Because of the way the network is parameterized,
the changes in the output $\frac{\partial f}{\partial\theta}$ also scale as
$O(\sigma^{L})$. This implies that the eNTK at any given time scales as

$K_{\theta}(\bm{x},\bm{x}^{\prime})=\sum_{\theta}\frac{\partial
f(\bm{x})}{\partial\theta}\frac{\partial
f(\bm{x}^{\prime})}{\partial\theta}=O(\sigma^{2L}).$ (21)

After appropriately rescaling learning rate to $\eta=\sigma^{-2L}$ we get

$\frac{df(\bm{x})}{dt}=-\eta\sum_{\mu}K_{\theta}(\bm{x},\bm{x}^{\mu})(f(\bm{x}^{\mu})-y^{\mu}).$
(22)

Under the assumption that $\sigma^{L}\ll 1$ and $y^{\mu}=O(1)$ so that the
error term is $O(1)$ we get that the output changes in time as $df/dt=O(1)$.

On the other hand, using the chain rule one can show that the features change
as a product of the gradient update and the features in the prior layer,
yielding the scaling

$\frac{dh^{(\ell)}}{dt}=\eta\frac{\sigma^{L}}{\sqrt{N}}=\frac{1}{\sigma^{L}\sqrt{N}}=(\alpha\sqrt{N})^{-1}.$
(23)

This gives us that the change in the features scales as
$(\alpha\sqrt{N})^{-1}$ while the change in the output scales as $O(1)$. Thus,
for $\alpha\sqrt{N}$ sufficiently small, the features can move dramatically.

###  C.2 Output Rescaling without Rescaling Weights

In the main text, we use the scale $\sigma$ at every layer to change the scale
of the output function. This relies on the homogeneity of the activation
function so that $W^{\ell}\to\sigma W^{\ell}$ for all $\ell$ leads to a
rescaling $f\to f\sigma^{L}$. This would not work for nonhomogenous
activations like $\phi(h)=\tanh(h)$. However, following Chizat et al. (2019);
Geiger et al. (2020a), we note that we can set all weights to be
$O_{\alpha}(1)$ and introduce the $\alpha$ only in the definition of the
neural network function

$\displaystyle
f=\frac{\alpha}{\sqrt{N}}\sum_{i=1}^{N}w^{L+1}_{i}\varphi(h^{L}_{i})\ ,\
h^{\ell}_{i}=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}W_{ij}^{\ell}\varphi(h^{\ell-1}_{j})\
,\ h_{i}^{1}=\frac{1}{\sqrt{D}}W_{ij}^{1}x_{j}.$ (24)

We note that all preactivations $h^{\ell}$ have scale $O_{\alpha}(1)$ for any
choice of nonlinearity, but that $f=\Theta_{\alpha}(\alpha)$. Several works
have established that the $\alpha\sim\frac{1}{\sqrt{N}}$ allows feature
learning even as the network approaches infinite width Mei et al. (2018); Yang
& Hu (2020); Bordelon & Pehlevan (2022). This is known as the mean field or
$\mu$-limit.

###  C.3 Kernel Alignment

In this section we comment on our choice of kernel alignment metric

$A(\bm{K}):=\frac{\bm{y}^{\top}\bm{K}\bm{y}}{\operatorname{Tr}\bm{K}|\bm{y}|^{2}}.$
(25)

For kernels that are diagonally dominant, such as those encountered in the
experiments, this metric is related to another alignment metric

$A_{F}(\bm{K}):=\frac{\bm{y}^{\top}\bm{K}\bm{y}}{|\bm{K}|_{F}|\bm{y}|^{2}}.$
(26)

Here $|\bm{K}|_{F}$ is the Frobenius norm of the Gram matrix of the kernel.
This metric was extensively used in Baratin et al. (2021). The advantage of
the first metric over the second is that one can quickly estimate the
denominator of $A(\bm{K})$ via Monte Carlo estimation of
$\langle\bm{u}^{\top}\bm{K}\bm{u}\rangle_{\bm{u}\sim\mathcal{N}(0,\bm{1})}$.

We use $A(\bm{K}_{f})$ as a measure of feature learning, as we have found that
this more finely captures elements of feature learning than other related
metrics. We list several metrics we tried that did not work.

One option for a representation-learning metric involves measuring the
magnitude of the change between the initial and final kernels,
$\bm{K}_{i},\bm{K}_{f}$:

$\Delta\bm{K}:=|\bm{K}_{f}-\bm{K}_{i}|_{F}.$ (27)

However, this is more sensitive to the raw parameter change than any task-
relevant data. If one instead were to normalize the kernels to be unit norm at
the beginning and the end, the modified metric

$\Delta\bm{K}:=\left|\frac{\bm{K}_{f}}{|\bm{K}_{f}|_{F}}-\frac{\bm{K}_{i}}{|\bm{K}_{i}|_{F}}\right|_{F}.$
(28)

This metric however remains remarkably flat over the whole range of
$\alpha,P$, as does the centered kernel alignment (CKA) of Cortes et al.
(2012)

$\mathrm{CKA}(\bm{K}_{i},\bm{K}_{f})=\frac{\operatorname{Tr}[\bm{K}_{i}^{c}\bm{K}_{f}^{c}]}{\sqrt{|\bm{K}_{i}^{c}|_{F}|\bm{K}_{f}^{c}|_{F}}},\quad\bm{K}^{c}=\bm{C}\bm{K}\bm{C},\quad\bm{C}=\bm{1}-\frac{1}{P}\vec{1}\,\vec{1}^{\;T}.$
(29)

Here $C$ is the centering matrix that subtracts off the mean components of the
kernel for a $P\times P$ kernel. This alignment metric has been shown to be
useful in comparing neural representations (Kornblith et al., 2019). For our
task, however, because the signal is low-dimensional, only a small set of
eigenspaces of the kernel align to this task. As a result, the CKA, which
counts all eigenspaces equally, appears to be too coarse to capture the low-
dimensional feature learning that is happening.

On the other hand, we find that $A(\bm{K}_{f})$ (with $\bm{K}_{f}$ given by
the eNTKf evaluated on a test set) can very finely detect alignment along the
task relevant directions. This produces a clear signal of feature learning at
small $\alpha$ and large $P$ as shown in Figure 2c.

$A(\bm{K}_{f})$ can be related to the centered kernel alignment between the
eNTKf and the (mean zero) task kernel $\bm{y}\bm{y}^{\top}$, where $\bm{y}$ is
a vector of draws from the population distribution $p(\bm{x},y)$.

###  C.4 Relationship Between Trained Network and Final Kernel

In general, the learned function contains contributions from the instantaneous
NTKs at every point in the training. Concretely, following Atanasov et al.
(2021) we have the following formula for the final network predictor $f(x)$

$f(x)=\int_{0}^{\infty}dt\,\bm{k}(x,t)\cdot\exp\left(-\int_{0}^{t}ds\bm{K}(s)\right)\bm{y},$
(30)

where $[\bm{k}(x,t)]_{\mu}=K(x,x_{\mu},t)$ and
$[\bm{K}(s)]_{\mu\nu}=K(x_{\mu},x_{\nu},s)$ and $[\bm{y}]_{\mu}=y_{\mu}$. In
general there are contributions from earlier kernels $\bm{k}(x,t)$ for
$t<\infty$ and so the function $f$ cannot always be written as a linear
combination of the final NTK $K_{f}$ on training data:
$f=\sum_{\mu}\alpha_{\mu}K_{f}(x,x_{\mu})$. However, as Vyas et al. (2022);
Atanasov et al. (2021) have shown, the final predictions of the network are
often well modeled by regression with the final NTK. We verify this for our
task in section 3.2.

##  Appendix D Generic Random Feature Model

###  D.1 Setting up the Problem: Feature Definitions

For a random kernel, $K({\bm{x}},{\bm{x}}^{\prime};\theta)$, we first compute
its Mercer decomposition

$\displaystyle\int d{\bm{x}}\
p({\bm{x}})K({\bm{x}},{\bm{x}}^{\prime};\theta)\phi_{k}({\bm{x}})=\lambda_{k}\phi_{k}({\bm{x}}^{\prime}).$
(31)

From the eigenvalues $\lambda_{k}$ and eigenfunctions $\phi_{k}$, we can
construct the square root

$\displaystyle
K^{1/2}({\bm{x}},{\bm{x}}^{\prime};\theta)=\sum_{k}\sqrt{\lambda_{k}}\phi_{k}({\bm{x}})\phi_{k}({\bm{x}}^{\prime}).$
(32)

Lastly, using $K^{1/2}$, we can get a feature map by projecting against a
static basis $\\{b_{k}\\}$ giving

$\displaystyle\psi_{k}({\bm{x}})=\int
d{\bm{x}}^{\prime}p({\bm{x}}^{\prime})K^{1/2}({\bm{x}},{\bm{x}}^{\prime};\theta)b_{k}({\bm{x}}^{\prime}).$
(33)

These features reproduce the kernel so that
$K({\bm{x}},{\bm{x}}^{\prime};\theta)=\sum_{k}\psi_{k}({\bm{x}})\psi_{k}({\bm{x}}^{\prime})$.
This can be observed from the following observation

$\displaystyle\psi_{k}({\bm{x}})$
$\displaystyle=\sum_{\ell}\sqrt{\lambda_{\ell}}\phi_{\ell}({\bm{x}})U_{\ell
k}\ ,\ U_{\ell k}=\left<\phi_{\ell}({\bm{x}})b_{k}({\bm{x}})\right>$ (34)
$\displaystyle\Rightarrow\sum_{k}\psi_{k}({\bm{x}})\psi_{k}({\bm{x}}^{\prime})$
$\displaystyle=\sum_{\ell,m}\sqrt{\lambda_{\ell}\lambda_{m}}\phi_{\ell}({\bm{x}})\phi_{m}({\bm{x}}^{\prime})\sum_{k}U_{\ell
k}U_{mk}=\sum_{\ell}\lambda_{\ell}\phi_{\ell}({\bm{x}}),\phi_{\ell}({\bm{x}}^{\prime})$
(35)

where the last line follows from the orthogonality of $U_{km}$ and recovers
$K(\bm{x},\bm{x}^{\prime};\theta)$.

###  D.2 Decomposition of Finite Width Features

We now attempt to characterize the variance in the features over the sample
distribution. We will first consider the case of a fixed realization of
$\bm{\theta}_{0}$ before providing a typical case analysis over random
$\bm{\theta}_{0}$. For a fixed initialization $\bm{\theta}_{0}$ we define the
following covariance matrices

$\displaystyle\bm{\Sigma}_{M}=\left<\bm{\psi}_{M}({\bm{x}})\bm{\psi}_{M}({\bm{x}})^{\top}\right>\in\mathbb{R}^{M\times
M}.$ (36)

where $\bm{\psi}_{M}$ are the truncated (but deterministic) features induced
by the deterministic infinite width kernel. We will mainly be interested in
the case where $M\to\infty$ and where the target function can be expressed as
the linear combination $y({\bm{x}})={\bm{w}}^{*}\cdot\bm{\psi}_{M}({\bm{x}})$
of these features. For example, in the case of our experiments on the sphere,
$\bm{\psi}_{M}$ could be the spherical harmonic functions. Further, in the
$M\to\infty$ limit, we will be able to express the target features $\bm{\psi}$
as linear combinations of the features $\bm{\psi}_{M}$

$\displaystyle\bm{\psi}({\bm{x}},\bm{\theta}_{0})={\bm{A}}(\bm{\theta}_{0})\bm{\psi}_{M}({\bm{x}})\
,\ {\bm{A}}(\bm{\theta})\in\mathbb{R}^{N_{\mathcal{H}}\times M}.$ (37)

The matrix $\bm{A}(\bm{\theta}_{0})$ are the coefficients of the decomposition
which can vary over initializations. Crucially $\bm{A}(\bm{\theta}_{0})$
projects to the subspace of dimension $N_{\mathcal{H}}$ where the finite width
features have variance over ${\bm{x}}$. The population risk for this
$\bm{\theta}_{0}$ has an irreducible component

$\displaystyle E_{g}(\bm{\theta}_{0})$
$\displaystyle=\left<\left({\bm{w}}^{*}\cdot\bm{\psi}_{M}-{\bm{w}}\cdot\bm{\psi}\right)^{2}\right>$
$\displaystyle\geq{\bm{w}}^{*\top}\left[\bm{\Sigma}_{M}-\bm{\Sigma}_{M}{\bm{A}}(\bm{\theta})^{\top}\left({\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}{\bm{A}}(\bm{\theta}_{0})^{\top}\right)^{-1}{\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}\right]{\bm{w}}^{*}.$
(38)

where the bound is tight for the optimal weights
${\bm{w}}=\left({\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}{\bm{A}}(\bm{\theta}_{0})^{\top}\right)^{-1}{\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}{\bm{w}}^{*}$.
The irreducible error is determined by a projection matrix which preserves the
subspace where the features $\bm{\psi}({\bm{x}},\bm{\theta}_{0})$ have
variance:
${\bm{I}}-{\bm{A}}(\bm{\theta})^{\top}\left({\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}{\bm{A}}(\bm{\theta}_{0})^{\top}\right)^{-1}{\bm{A}}(\bm{\theta}_{0})\bm{\Sigma}_{M}$.
In general, this will preserve some fraction of the variance in the target
function, but some variance in the target function will not be expressible by
linear combinations of the features $\bm{\psi}({\bm{x}},\bm{\theta})$. We
expect that random finite width $N$ neural networks will have unexplained
variance in the target function on the order $\sim 1/N$.

###  D.3 Gaussian Covariate Model

Following prior works on learning curves for kernel regression (Bordelon et
al., 2020; Canatar et al., 2021; Loureiro et al., 2021), we will approximate
the learning problem with a Gaussian covariates model with matching second
moments.

The features $\bm{\psi}_{M}({\bm{x}})$ will be treated as Gaussian over random
draws of datapoints. We will assume centered features. We decompose the
features in the orthonormal basis ${\bm{b}}({\bm{x}})$, which we approximate
as a Gaussian vector ${\bm{b}}\sim\mathcal{N}(0,{\bm{I}})$.

$\displaystyle f$ $\displaystyle=\bm{\psi}(\bm{\theta}_{0})\cdot{\bm{w}}\ ,\
y=\bar{\bm{\psi}}_{M}\cdot{\bm{w}}^{*}$ $\displaystyle{\psi}_{M}$
$\displaystyle=\bm{\Sigma}_{s}^{1/2}{\bm{b}}\ ,\
\bm{\psi}(\bm{\theta}_{0})={\bm{A}}(\bm{\theta}_{0})^{\top}\bm{\psi}_{M}+\bm{\Sigma}_{\epsilon}^{1/2}\bm{\epsilon}$
$\displaystyle{\bm{b}}$ $\displaystyle\sim\mathcal{N}(0,{\bm{I}})\ ,\
\bm{\epsilon}\sim\mathcal{N}(0,{\bm{I}})$ (39)

This is a special case of the Gaussian covariate model introduced by Loureiro
et al. (2021) and subsumes the popular two-layer random feature models (Mei &
Montanari, 2022; Adlam & Pennington, 2020b) as a special case. In a subsequent
section, we go beyond Loureiro et al. (2021) by computing typical case
learning curves over Gaussian ${\bm{A}}(\bm{\theta}_{0})$ matrices. In
particular, we have for the two layer random feature model in the proportional
asymptotic limit $P,N,D\to\infty$ with $P/D=O(1)$ and $P/N=O(1)$ with
$\bm{\psi}({\bm{x}})=\phi(\bm{F}{\bm{x}}_{\mu})$ for fixed feature matrix
${\bm{F}}\in\mathbb{R}^{N\times D}$ nonlinearity $\phi$ and
${\bm{x}}={\bm{b}}\sim\mathcal{N}(0,D^{-1}\bm{I})$

$\displaystyle\bm{\Sigma}_{M}$ $\displaystyle={\bm{I}}\ ,\
\bm{\Sigma}_{\epsilon}=c_{*}^{2}{\bm{I}}\ ,\ {\bm{A}}^{\top}=c_{1}{\bm{F}}$
$\displaystyle c_{1}$
$\displaystyle=\left<z\phi(z)\right>_{z\sim\mathcal{N}(0,1)}\ ,\
c_{*}^{2}=\left<\phi(z)^{2}\right>_{z\sim\mathcal{N}(0,1)}-c_{1}^{2}.$ (40)

We refer readers to Hu & Lu (2020) for a discussion of this equivalence
between random feature regression and this Gaussian covariate model.

###  D.4 Replica Calculation of the Learning Curve

To analyze the typical case performance of kernel regression, we define the
following partition function which is dominated

$\displaystyle Z[\mathcal{D},\bm{\theta}_{0}]$ $\displaystyle=\int
d{\bm{w}}\exp\left(-\frac{\beta}{2\lambda}\sum_{\mu=1}^{P}[{\bm{w}}\cdot\bm{\psi}_{\mu}-{\bm{w}}^{*}\cdot{\bm{\psi}}_{M,\mu}]^{2}-\frac{\beta}{2}|{\bm{w}}|^{2}-\frac{J\beta
M}{2}E_{g}({\bm{w}})\right)$ $\displaystyle
E_{g}({\bm{w}})=\frac{1}{M}|\bm{\Sigma}_{M}^{1/2}{\bm{w}}^{*}-\bm{\Sigma}_{M}^{1/2}{\bm{A}}(\bm{\theta}_{0}){\bm{w}}|^{2}+\frac{1}{M}{\bm{w}}^{\top}\bm{\Sigma}_{\epsilon}{\bm{w}}.$
(41)

For proper normalization, we assume that
$\left<\bm{\psi}_{M}\bm{\psi}_{M}^{\top}\right>=\frac{1}{M}\bm{\Sigma}_{M}$
and
$\left<\bm{\epsilon}\bm{\epsilon}^{\top}\right>=\frac{1}{M}\bm{\Sigma}_{\epsilon}$.
We note that in the $\beta\to\infty$ limit, the partition function is
dominated by the unique minimizer of the regularized least squares objective
(Canatar et al., 2021; Loureiro et al., 2021). Further, for a fixed
realization of $\bm{\theta}_{0}$ the average generalization error over
datasets $\mathcal{D}$ can be computed by differentiation of the source term
$J$

$\displaystyle\frac{2}{\beta M}\frac{\partial}{\partial J}|_{J=0}\left<\ln
Z[\mathcal{D},\bm{\theta}_{0}]\right>_{\mathcal{D}}$ (42)
$\displaystyle=\left<\frac{1}{Z}\int
d{\bm{w}}\exp\left(-\frac{\beta}{2\lambda}\sum_{\mu=1}^{P}[{\bm{w}}\cdot\bm{\psi}_{\mu}-{\bm{w}}^{*}\cdot{\bm{\psi}}_{M,\mu}]^{2}-\frac{\beta}{2}|{\bm{w}}|^{2}\right)E_{g}({\bm{w}})\right>_{\mathcal{D}}.$

Thus the $\beta\to\infty$ limit of the above quantity will give the expected
generalization error of the risk minimizer. We see the need to average the
quantity $\ln Z$ over realizations of datasets $\mathcal{D}$. For this, we
resort to the replica trick $\left<\ln Z\right>=\lim_{n\to
0}\frac{1}{n}\ln\left<Z^{n}\right>$. We will compute the integer moments
$\left<Z^{n}\right>$ for integer $n$ and then analytically continue the
resulting expressions to $n\to 0$ under a symmetry ansatz. The replicated
partition function thus has the form

$\displaystyle\left<Z^{n}\right>=\int\prod_{a=1}^{n}d{\bm{w}}^{a}\mathbb{E}_{\\{{\bm{b}}_{\mu},\bm{\epsilon}_{\mu}\\}}$
$\displaystyle\exp\left(-\frac{\beta}{2\lambda}\sum_{\mu=1}^{P}\sum_{a=1}^{n}[{\bm{w}}^{a}\cdot\bm{\psi}_{\mu}-{\bm{w}}^{*}\cdot{\bm{\psi}}_{M,\mu}]^{2}-\frac{\beta}{2}\sum_{a=1}^{n}|{\bm{w}}^{a}|^{2}\right)$
(43) $\displaystyle\times\exp\left(-\frac{J\beta
M}{2}\sum_{a=1}^{n}E_{g}({\bm{w}}^{a})\right).$

We now need to perform the necessary average over the random realizations of
data points $\mathcal{D}=\\{{\bm{b}}_{\mu},\bm{\epsilon}_{\mu}\\}$. We note
that the scalar quantities
$h^{a}_{\mu}={\bm{w}}^{a}\cdot\bm{\psi}_{\mu}-{\bm{w}}^{*}\cdot{\bm{\psi}}_{M,\mu}$
are Gaussian with mean zero and covariance

$\displaystyle\left<h^{a}_{\mu}h^{b}_{\nu}\right>$
$\displaystyle=\delta_{\mu\nu}Q_{ab}$ $\displaystyle Q_{ab}$
$\displaystyle=\frac{1}{M}\left({\bm{A}}(\bm{\theta}_{0}){\bm{w}}^{a}-{\bm{w}}^{*}\right)\bm{\Sigma}_{M}\left({\bm{A}}(\bm{\theta}_{0}){\bm{w}}^{a}-{\bm{w}}^{*}\right)+\frac{1}{M}{\bm{w}}^{a}\bm{\Sigma}_{\epsilon}{\bm{w}}^{b}.$
(44)

We further see that the generalization error in replica $a$ is
$E_{g}({\bm{w}}^{a})=Q_{aa}$. Performing the Gaussian integral over
$\\{h^{a}_{\mu}\\}$ gives

$\displaystyle\left<Z^{n}\right>\propto\int\prod_{a}$ $\displaystyle
d{\bm{w}}^{a}\prod_{ab}dQ_{ab}d\hat{Q}_{ab}\exp\left(-\frac{P}{2}\ln\det\left(\lambda{\bm{I}}+\beta{\bm{Q}}\right)-\frac{J\beta
M}{2}\text{Tr}{\bm{Q}}-\frac{\beta}{2}\sum_{a}|{\bm{w}}^{a}|^{2}\right)$
$\displaystyle\exp\left(\frac{1}{2}\sum_{ab}\hat{Q}_{ab}\left(MQ_{ab}-[{\bm{A}}(\bm{\theta}_{0}){\bm{w}}^{a}-{\bm{w}}^{*}]^{\top}\bm{\Sigma}_{M}[{\bm{A}}(\bm{\theta}_{0}){\bm{w}}^{a}-{\bm{w}}^{*}]+{\bm{w}}^{a}\bm{\Sigma}_{\epsilon}{\bm{w}}^{b}\right)\right).$

We introduced the Lagrange multipliers $\hat{{\bm{Q}}}$ which enforce the
definition of order parameters ${\bm{Q}}$. We now integrate over
${\bm{W}}=\text{Vec}\\{{\bm{w}}^{a}\\}_{a=1}^{n}$. We let
$\bm{\tilde{\Sigma}}_{s}={\bm{A}}^{\top}\bm{\Sigma}_{M}{\bm{A}}$

$\displaystyle\int d{\bm{W}}$
$\displaystyle\exp\left(-\frac{1}{2}{\bm{W}}\left[\beta{\bm{I}}+\hat{{\bm{Q}}}\otimes[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}]\right]{\bm{W}}\right)$
$\displaystyle\exp\left({\bm{W}}^{\top}\left[\hat{{\bm{Q}}}\otimes{\bm{I}}\right]\left(\bm{1}\otimes{\bm{A}}^{\top}\bm{\Sigma}_{s}{\bm{w}}^{*}\right)\right)$
$\displaystyle=$
$\displaystyle\exp\left(\frac{1}{2}\left(\bm{1}\otimes{\bm{A}}^{\top}\bm{\Sigma}_{s}{\bm{w}}^{*}\right)^{\top}\left[\hat{{\bm{Q}}}\otimes{\bm{I}}\right]\left[\beta{\bm{I}}+\hat{{\bm{Q}}}\otimes[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}]\right]^{-1}\left[\hat{{\bm{Q}}}\otimes{\bm{I}}\right]\left(\bm{1}\otimes{\bm{A}}^{\top}\bm{\Sigma}_{s}{\bm{w}}^{*}\right)\right)$
$\displaystyle\exp\left(-\frac{1}{2}\ln\det\left[\beta{\bm{I}}+\hat{{\bm{Q}}}\otimes[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}]\right]\right).$
(45)

To take the $n\to 0$ limit, we make the replica symmetry ansatz

$\displaystyle\beta{\bm{Q}}=q{\bm{I}}+q_{0}\bm{1}\bm{1}^{\top}\ ,\
\beta^{-1}\hat{{\bm{Q}}}=\hat{q}{\bm{I}}+\hat{q}_{0}\bm{1}\bm{1}^{\top},$ (46)

which is well motivated since this is a convex optimization problem. Letting
$\alpha=P/N$, we find that under the RS ansatz the replicated partition
function has the form

$\displaystyle\left<Z^{n}\right>$ $\displaystyle=\int
dqdq_{0}d\hat{q}d\hat{q}_{0}\exp\left(\frac{nM}{2}S[q,q_{0},\hat{q},\hat{q}_{0}]\right)$
(47) $\displaystyle S$
$\displaystyle=q\hat{q}+q_{0}\hat{q}+q\hat{q}_{0}-\alpha\left[\ln(\lambda+q)+\frac{q_{0}}{\lambda+q}\right]$
$\displaystyle-\frac{\beta}{M}{\bm{w}}^{*}[\hat{q}\bm{\Sigma}_{M}]{\bm{w}}^{*}+\frac{\beta}{M}{\bm{w}}^{*}[\hat{q}\tilde{\bm{\Sigma}}_{s}]{\bm{A}}^{\top}{\bm{G}}{\bm{A}}[\hat{q}\tilde{\bm{\Sigma}}_{s}]{\bm{w}}^{*}$
$\displaystyle-\frac{1}{M}\ln\det{\bm{G}}-\frac{1}{M}\hat{q}_{0}\text{Tr}{\bm{G}}[\tilde{\bm{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}]-J(q+q_{0})$
$\displaystyle{\bm{G}}$
$\displaystyle=\left({\bm{I}}+\hat{q}[\tilde{\bm{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}]\right)^{-1}.$
(48)

In a limit where $\alpha=P/M$ is $O(1)$, then this $S$ is intensive
$S=O_{M}(1)$. We can thus appeal to saddle point integration (method of
steepest descent) to compute the set of order parameters which have dominant
contribution to the free energy.

$\displaystyle\left<Z^{n}\right>$ $\displaystyle=\int
dqdq_{0}d\hat{q}d\hat{q}_{0}\exp\left(\frac{nM}{2}S[q,q_{0},\hat{q},\hat{q}_{0}]\right)\sim\exp\left(\frac{nM}{2}S[q^{*},q_{0}^{*},\hat{q}^{*},\hat{q}_{0}^{*}]\right)$
$\displaystyle\implies\left<\ln Z\right>$
$\displaystyle=\frac{M}{2}S[q^{*},q_{0}^{*},\hat{q}^{*},\hat{q}_{0}^{*}].$
(49)

The order parameters $q^{*},q_{0}^{*},\hat{q}^{*},\hat{q}_{0}^{*}$ are defined
via the saddle point equations $\frac{\partial S}{\partial q}=\frac{\partial
S}{\partial q_{0}}=\frac{\partial S}{\partial\hat{q}}=\frac{\partial
S}{\partial\hat{q}_{0}}=0$. For our purposes, it suffices to analyze two of
these equations

$\displaystyle\frac{\partial S}{\partial q_{0}}$
$\displaystyle=\hat{q}-\frac{\alpha}{\lambda+q}-J=0,$
$\displaystyle\frac{\partial S}{\partial\hat{q}_{0}}$
$\displaystyle=q-\frac{1}{M}\text{Tr}{\bm{G}}\left[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right]=0.$
(50)

We can now take the zero temperature ($\beta\to\infty$) limit to solve for the
generalization error

$\displaystyle E_{g}$ $\displaystyle=-\frac{\partial}{\partial
J}|_{J=0}\lim_{\beta\to\infty}\frac{1}{\beta}F$
$\displaystyle=\frac{1}{M}\partial_{J}{\bm{w}}^{*}\left[\hat{q}\bm{\Sigma}_{M}-\hat{q}^{2}\bm{\Sigma}_{M}{\bm{A}}^{\top}{\bm{G}}{\bm{A}}\bm{\Sigma}_{M}\right]{\bm{w}}^{*}.$
(51)

We see that we need to compute the $J$ derivatives on $\hat{q}$. We let
$\kappa=\lambda+q$ and note

$\displaystyle\partial_{J}\hat{q}$
$\displaystyle=-\alpha\kappa^{-2}\partial_{J}\kappa+1$
$\displaystyle\partial_{J}\kappa$ $\displaystyle=-\partial_{J}\hat{q}\
\frac{1}{M}\text{Tr}{\bm{G}}^{2}\left[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right]^{2}=-\left(-\alpha\kappa^{-2}\partial_{J}\kappa+1\right)\frac{1}{M}\text{Tr}{\bm{G}}^{2}\left[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right]^{2}.$
(52)

We solve the equation for $\partial_{J}\kappa$ which gives
$\partial_{J}\kappa=-\frac{\kappa^{2}}{\alpha}\frac{\gamma}{1-\gamma}$ where
$\gamma=\frac{\alpha}{\kappa^{2}}\frac{1}{M}\text{Tr}{\bm{G}}^{2}\left[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right]^{2}$.
With this definition we have
$\partial_{J}\hat{q}=1+\frac{\gamma}{1-\gamma}=\frac{1}{1-\gamma}$.

$\displaystyle E_{g}$
$\displaystyle=\frac{1}{1-\gamma}\frac{1}{M}{\bm{w}}^{*}\bm{\Sigma}_{M}^{1/2}\left[{\bm{I}}-2\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}+\hat{q}^{2}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}\left[\bm{\tilde{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right]{\bm{G}}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}\right]\bm{\Sigma}_{M}^{1/2}{\bm{w}}^{*}$
$\displaystyle=\frac{1}{1-\gamma}\frac{1}{M}{\bm{w}}^{*}\bm{\Sigma}_{M}^{1/2}\left[{\bm{I}}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}^{2}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}\right]\bm{\Sigma}_{M}^{1/2}{\bm{w}}^{*}.$
(53)

This reproduces the derived expression from Loureiro et al. (2021). The
matching covariance $\tilde{\bm{\Sigma}}_{s}=\bm{\Sigma}_{M}$ and zero
feature-noise limit $\bm{\Sigma}_{\epsilon}=0$ recovers the prior results of
Bordelon et al. (2020); Canatar et al. (2021); Simon et al. (2021). In
general, this error will asymptote to the the irreducible error

$\displaystyle\lim_{P\to\infty}E_{g}=\frac{1}{M}{\bm{w}}^{*}\left[\bm{\Sigma}_{M}-{\bm{\Sigma}}_{s}{\bm{A}}^{\top}\left(\tilde{\bm{\Sigma}}_{s}+\bm{\Sigma}_{\epsilon}\right)^{-1}{\bm{A}}{\bm{\Sigma}}_{s}\right]{\bm{w}}^{*}.$
(54)

We see that this recovers the minimal possible error in the $P\to\infty$
limit. The derived learning curves depend on the instance of random initial
condition $\bm{\theta}_{0}$. To get the average case performance, we take an
additional average of this expression over $\bm{\theta}_{0}$

$\displaystyle\mathbb{E}_{\bm{\theta}_{0}}E_{g}(\bm{\theta}_{0})=\mathbb{E}_{\theta_{0}}\frac{1}{1-\gamma}\frac{1}{M}{\bm{w}}^{*}\bm{\Sigma}_{M}^{1/2}\left[{\bm{I}}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}-\hat{q}{\bm{\Sigma}}_{s}^{1/2}{\bm{A}}^{\top}{\bm{G}}^{2}{\bm{A}}{\bm{\Sigma}}_{s}^{1/2}\right]\bm{\Sigma}_{M}^{1/2}{\bm{w}}^{*}.$
(55)

This average is complicated since $\gamma,\hat{q},{\bm{G}}$ all depend on
$\bm{\theta}_{0}$. In the next section we go beyond this analysis to try
average case analysis for random Gaussian ${\bm{A}}$.

###  D.5 Quenched Average over Gaussian A

In this section we will define a distribution of features which allows an
exact asymptotic prediction over random realizations of disorder
$\bm{\theta}_{0}$ and datasets $\mathcal{D}$. This is a nontrivial extension
of the result of Loureiro et al. (2021) since the number of necessary saddle
point equations to be solved doubles from two to four. However, this more
complicated theory allows us to exactly compute the expectation in equation 55
under an ansatz for the random matrix ${\bm{A}}$. We construct our features
with

$\displaystyle\bm{\psi}|{\bm{A}}$
$\displaystyle=\frac{1}{\sqrt{N}}{\bm{A}}^{\top}{\bm{\psi}}_{M}+\bm{\Sigma}_{\epsilon}^{1/2}\bm{\epsilon}\
,\ A_{ij}\sim\mathcal{N}(0,\sigma^{2}).$

We will now perform an approximate average over both datasets $\mathcal{D}$
and realizations of $\bm{A}$

$\displaystyle\left<Z^{n}\right>=\int\prod_{a=1}^{n}d{\bm{w}}^{a}\mathbb{E}_{\\{{\bm{b}}_{\mu},\bm{\epsilon}_{\mu},{\bm{A}}\\}}$
$\displaystyle\exp\left(-\frac{\beta}{2\lambda}\sum_{\mu=1}^{P}\sum_{a=1}^{n}[{\bm{w}}^{a}\cdot\bm{\psi}_{\mu}-{\bm{w}}^{*}\cdot{\bm{\psi}}_{M,\mu}]^{2}-\frac{\beta}{2}\sum_{a=1}^{n}|{\bm{w}}^{a}|^{2}\right)$
(56)
$\displaystyle\times\exp\left(-\frac{JM\beta}{2}\sum_{a=1}^{n}E_{g}({\bm{w}}^{a})\right).$

As before, we first average over ${\bm{b}}_{\mu},\bm{\epsilon}_{\mu}|A$ and
define order parameters $Q_{ab}$ as before.

$\displaystyle\left<Z^{n}\right>=\int\prod_{a}$ $\displaystyle
d{\bm{w}}^{a}\prod_{ab}dQ_{ab}d\hat{Q}_{ab}\exp\left(-\frac{P}{2}\ln\det\left(\lambda{\bm{I}}+\beta{\bm{Q}}\right)-\frac{J\beta
M}{2}\text{Tr}{\bm{Q}}-\frac{\beta}{2}\sum_{a}|{\bm{w}}^{a}|^{2}\right)$
$\displaystyle\mathbb{E}_{\\{{\bm{g}}^{a}\\}}$
$\displaystyle\exp\left(\frac{1}{2}\sum_{ab}\hat{Q}_{ab}\left(MQ_{ab}-[{\bm{g}}^{a}-{\bm{w}}^{*}]^{\top}\bm{\Sigma}_{M}[{\bm{g}}^{b}-{\bm{w}}^{*}]+{\bm{w}}^{a}\bm{\Sigma}_{\epsilon}{\bm{w}}^{b}\right)\right).$

where we defined the fields
${\bm{g}}^{a}=\frac{1}{\sqrt{N}}{\bm{A}}{\bm{w}}^{a}$ which are mean zero
Gaussian with covariance
$\left<{\bm{g}}^{a}{\bm{g}}^{b\top}\right>=V_{ab}{\bm{I}}$ where
$V_{ab}=\frac{\sigma^{2}}{N}{\bm{w}}^{a}\cdot{\bm{w}}^{b}$. Performing the
Gaussian integral over ${\bm{G}}=\text{Vec}\\{{\bm{g}}^{a}\\}$, we find

$\displaystyle\int\prod_{a}d{\bm{g}}^{a}\exp\left(-\frac{1}{2}{\bm{G}}\left[{\bm{I}}\otimes{\bm{V}}^{-1}+\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}\right]{\bm{G}}+(\bm{\Sigma}_{M}{\bm{w}}^{*}\otimes\hat{{\bm{Q}}}\bm{1}){\bm{G}}-\frac{1}{2}\ln\det\left({\bm{I}}\otimes{\bm{V}}\right)\right)$
$\displaystyle=\exp\left(\frac{1}{2}(\bm{\Sigma}_{M}{\bm{w}}^{*}\otimes\hat{{\bm{Q}}}\bm{1})\left[{\bm{I}}\otimes{\bm{V}}^{-1}+\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}\right]^{-1}(\bm{\Sigma}_{M}{\bm{w}}^{*}\otimes\hat{{\bm{Q}}}\bm{1})-\frac{1}{2}\ln\det\left({\bm{I}}+\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}{\bm{V}}\right)\right).$
(57)

Next, we need to integrate over ${\bm{W}}=\text{Vec}\\{{\bm{w}}^{a}\\}$ which
gives

$\displaystyle\int
d{\bm{W}}\exp\left(-\frac{1}{2}{\bm{W}}\left[\beta{\bm{I}}+\sigma^{2}{\bm{I}}\otimes\hat{{\bm{V}}}+\bm{\Sigma}_{\epsilon}\otimes\hat{{\bm{Q}}}\right]{\bm{W}}\right)=\exp\left(-\frac{1}{2}\ln\det\left[\beta{\bm{I}}+\sigma^{2}{\bm{I}}\otimes\hat{{\bm{V}}}+\bm{\Sigma}_{\epsilon}\otimes\hat{{\bm{Q}}}\right]\right).$
(58)

Now the replicated partition function has the form

$\displaystyle\left<Z^{n}\right>=$ $\displaystyle\int
d{\bm{Q}}d\hat{{\bm{Q}}}d{\bm{V}}d\hat{{\bm{V}}}\exp\left(\frac{M}{2}\text{Tr}[{\bm{Q}}\hat{{\bm{Q}}}+\eta{\bm{V}}\hat{{\bm{V}}}]-\frac{J\beta
M}{2}\text{Tr}{\bm{Q}}-\frac{P}{2}\ln\det\left[\lambda{\bm{I}}+\beta{\bm{Q}}\right]\right)$
(59)
$\displaystyle\times\exp\left(-\frac{1}{2}({\bm{w}}^{*}\otimes\bm{1})^{\top}[\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}]({\bm{w}}^{*}\otimes\bm{1})\right)$
$\displaystyle\times\exp\left(\frac{1}{2}(\bm{\Sigma}_{M}{\bm{w}}^{*}\otimes\hat{{\bm{Q}}}\bm{1})^{\top}\left[{\bm{I}}\otimes{\bm{V}}^{-1}+\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}\right]^{-1}(\bm{\Sigma}_{M}{\bm{w}}^{*}\otimes\hat{{\bm{Q}}}\bm{1})\right)$
$\displaystyle\times\exp\left(-\frac{1}{2}\ln\det\left[{\bm{I}}+\bm{\Sigma}_{M}\otimes\hat{{\bm{Q}}}{\bm{V}}\right]-\frac{1}{2}\ln\det\left[\beta{\bm{I}}+\sigma^{2}{\bm{I}}\otimes\hat{{\bm{V}}}+\bm{\Sigma}_{\epsilon}\otimes\hat{{\bm{Q}}}\right]\right).$

Now we make a replica symmetry ansatz on the order parameters
${\bm{Q}},\hat{{\bm{Q}}},{\bm{V}},\hat{{\bm{V}}}$

$\displaystyle\beta{\bm{Q}}=q{\bm{I}}+q_{0}\bm{1}\bm{1}^{\top}\ ,\
\beta{\bm{V}}=v{\bm{I}}+v_{0}\bm{1}\bm{1}^{\top}$
$\displaystyle\beta^{-1}\hat{{\bm{Q}}}=\hat{q}{\bm{I}}+\hat{q}_{0}\bm{1}\bm{1}^{\top}\
,\ \beta^{-1}\hat{{\bm{V}}}=\hat{v}{\bm{I}}+\hat{v}_{0}\bm{1}\bm{1}^{\top}.$
(60)

We introduce the shorthand for normalized trace of a matrix ${\bm{G}}$ as
$\text{tr}\ {\bm{G}}=\frac{1}{M}\text{Tr}{\bm{G}}$. Under the replica symmetry
ansatz, we find the following free energy

$\displaystyle\frac{2}{M}\left<\ln Z\right>$
$\displaystyle=q\hat{q}+q_{0}\hat{q}+q\hat{q}_{0}+\eta(v\hat{v}+v_{0}\hat{v}+v\hat{v}_{0})-J(q+q_{0})-\alpha\left[\ln(\lambda+q)+\frac{q_{0}}{\lambda+q}\right]$
$\displaystyle-\frac{\beta}{M}{\bm{w}}^{*}[\hat{q}\bm{\Sigma}_{M}]{\bm{w}}^{*}+\frac{\beta}{M}{\bm{w}}^{*}[\hat{q}\bm{\Sigma}_{M}][v^{-1}{\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-1}[\hat{q}\bm{\Sigma}_{M}]$
$\displaystyle-\text{tr}\log\left[{\bm{I}}+\hat{q}v\bm{\Sigma}_{M}\right]-(\hat{q}_{0}v+\hat{q}v_{0})\
\text{tr}[{\bm{I}}+\hat{q}v\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}$
$\displaystyle-\text{tr}\log\left[{\bm{I}}+\sigma^{2}\hat{v}{\bm{I}}+\bm{\Sigma}_{\epsilon}\hat{q}\right]-\text{tr}\left[{\bm{I}}+\sigma^{2}\hat{v}{\bm{I}}+\bm{\Sigma}_{\epsilon}\hat{q}\right]^{-1}\left[\hat{v}_{0}\sigma^{2}{\bm{I}}+\hat{q}_{0}\bm{\Sigma}_{\epsilon}\right].$
(61)

Letting $F=2M^{-1}\left<\ln Z\right>$, the saddle point equations read

$\displaystyle\frac{\partial F}{\partial q_{0}}$
$\displaystyle=\hat{q}-\frac{\alpha}{\lambda+q}-J=0,$
$\displaystyle\frac{\partial F}{\partial\hat{q}_{0}}$ $\displaystyle=q-v\
\text{tr}[{\bm{I}}+\hat{q}v\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}-\text{tr}[{\bm{I}}+\sigma^{2}\hat{v}{\bm{I}}+\bm{\Sigma}_{\epsilon}\hat{q}]^{-1}\bm{\Sigma}_{\epsilon}=0,$
$\displaystyle\frac{\partial F}{\partial v_{0}}$
$\displaystyle=\eta\hat{v}-\hat{q}\
\text{tr}[{\bm{I}}+\hat{q}v\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}=0,$
$\displaystyle\frac{\partial F}{\partial\hat{v}_{0}}$ $\displaystyle=\eta
v-\sigma^{2}\
\text{tr}[{\bm{I}}+\sigma^{2}\hat{v}{\bm{I}}+\bm{\Sigma}_{\epsilon}\hat{q}]^{-1}=0.$
(62)

Now the generalization error can be determined from

$\displaystyle E_{g}$ $\displaystyle=-\frac{\partial}{\partial
J}\lim_{\beta\to\infty}\frac{2}{\beta M}\left<\ln
Z\right>=\partial_{J}\frac{1}{M}{\bm{w}}^{*}\left[\hat{q}\bm{\Sigma}_{M}-(\hat{q}\bm{\Sigma}_{M})[v^{-1}{\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-1}(\hat{q}\bm{\Sigma}_{M})\right]{\bm{w}}^{*}.$
(63)

We see that it is necessary to compute $\partial_{J}\hat{q}$ and
$\partial_{J}v$ in order to obtain the final result. For simplicity, we set
$\sigma^{2}=1$. The equations for the source derivatives are

$\displaystyle\partial_{J}\hat{q}=-\frac{\alpha}{(\lambda+q)^{2}}\partial_{J}q+1,$
$\displaystyle\partial_{J}q=-\text{tr}[{\bm{I}}+v\hat{q}\bm{\Sigma}_{M}]^{-2}[-\partial_{J}v{\bm{I}}+v^{2}\partial_{J}\hat{q}\bm{\Sigma}_{M}]\bm{\Sigma}_{M}-\text{tr}[{\bm{I}}+\hat{v}{\bm{I}}+\hat{q}\bm{\Sigma}_{\epsilon}]^{-2}\bm{\Sigma}_{\epsilon}[\partial_{J}\hat{v}{\bm{I}}+\partial_{J}\hat{q}\bm{\Sigma}_{\epsilon}],$
$\displaystyle\eta\partial_{J}\hat{v}=-\text{tr}[{\bm{I}}+v\hat{q}\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}[-\partial_{J}\hat{q}{\bm{I}}+\hat{q}^{2}\partial_{J}v\bm{\Sigma}_{M}],$
$\displaystyle\eta\partial_{J}v=-\text{tr}[{\bm{I}}+\hat{v}{\bm{I}}+\bm{\Sigma}_{\epsilon}\hat{q}]^{-2}[\partial_{J}\hat{v}{\bm{I}}+\partial_{J}\hat{q}\bm{\Sigma}_{\epsilon}].$
(64)

Once the value of the order parameters $(q,\hat{q},v,\hat{v})$ have been
determined, these source derivatives can be obtained by solving a $4\times 4$
linear system. Examples of these solutions are provided in Figure 16.

####  D.5.1 Asymptotics in Underparameterized Regime

We can compute the asymptotic ($\alpha\to\infty$) generalization error due to
the random projection ${\bm{A}}$ in the limit of $\bm{\Sigma}_{\epsilon}=0$.
First, note that if $\hat{v}\to O_{\alpha}(1)$, then the asymptotic error
would be zero. Therefore, we will assume that $\hat{v}\sim a\alpha^{c}$ for
some $a,c>0$. The saddle point equations give the following asymptotic
conditions

$\displaystyle\hat{q}\sim\frac{\alpha}{\lambda}\ ,\ \eta\sim\hat{q}\
\text{tr}[\hat{v}{\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}$
$\displaystyle\implies\eta=\text{tr}[\lambda
a\alpha^{c-1}\bm{I}+\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}.$ (65)

For $0<\eta<1$, this equation can only be satisfied as $\alpha\to\infty$ if
$c=1$ so that $\hat{v}$ has the same scaling with $\alpha$ as $\hat{q}$. If
$c<1$ then we could get the equation $\eta=1$. If $c>1$, then the equation
would give $\eta=0$. The constant $a$ solves the equation

$\displaystyle\eta=\text{tr}[\lambda
a{\bm{I}}+\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}.$ (66)

Using this fact, our order parameters satisfy the following large $\alpha$
scalings

$\displaystyle\hat{q}\sim\frac{\alpha}{\lambda}\ ,\ q\sim 0\ ,\ \hat{v}\sim
a\alpha\ ,\ v\sim 0.$ (67)

The source derivative equations simplify to $\partial_{J}\hat{q}\sim 1\ ,\
\partial_{J}q\sim 0$ and

$\displaystyle\eta\partial_{J}\hat{v}$
$\displaystyle\sim(\hat{v}\partial_{J}\hat{q}+\hat{q}\partial_{J}\hat{v})\
\text{tr}[\hat{v}\bm{I}+\hat{q}\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}-\hat{q}\hat{v}\text{tr}[\hat{v}{\bm{I}}+\hat{q}\bm{\Sigma}]^{-2}[\partial_{J}\hat{v}\bm{\Sigma}_{M}+\partial_{J}\hat{q}\bm{\Sigma}_{M}^{2}]$
$\displaystyle\sim\hat{q}^{2}\text{tr}[\hat{v}{\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}^{2}\partial_{J}\hat{v}+\hat{v}^{2}\text{tr}[\hat{v}{\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}$
$\displaystyle\Rightarrow\partial_{J}\hat{v}$
$\displaystyle\sim\frac{\text{tr}[{\bm{I}}+a^{-1}\lambda^{-1}\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}}{\eta-\text{tr}[a\lambda{\bm{I}}+\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}^{2}}.$
(68)

We note that $\partial_{J}\hat{v}$ only depends on the product $a\lambda$
which is an implicit function of $\eta$ and $\bm{\Sigma}_{M}$. The
generalization error is
$E_{g}=\frac{1}{M}\partial_{J}\hat{q}(1+\hat{v}){\bm{w}}^{*}[(1+\hat{v}){\bm{I}}+\hat{q}\bm{\Sigma}]^{-1}\bm{\Sigma}_{M}{\bm{w}}^{*}$

$\displaystyle E_{g}$
$\displaystyle\sim\frac{1}{M}{\bm{w}}^{*}\left[{\bm{I}}+a^{-1}\lambda^{-1}\bm{\Sigma}_{M}\right]^{-2}\bm{\Sigma}_{M}{\bm{w}}^{*}$
(69) $\displaystyle\quad+\frac{1}{M}{\bm{w}}^{*}[\lambda
a{\bm{I}}+\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}^{2}{\bm{w}}^{*}\times\frac{\text{tr}[{\bm{I}}+a^{-1}\lambda^{-1}\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}}{\eta-\text{tr}[a\lambda{\bm{I}}+\bm{\Sigma}_{M}]^{-2}\bm{\Sigma}_{M}^{2}}.$

We see that in the generic case, the asymptotic error has a nontrivial
dependence on the task ${\bm{w}}^{*}$ and the correlation structure
$\bm{\Sigma}_{M}$. To gain more intuition, we will now consider the special
case of isotropic features $\bm{\Sigma}_{M}={\bm{I}}$. In this case, we have
$\eta=\frac{1}{1+\lambda a}$ so that $\lambda a=\frac{1-\eta}{\eta}$. This
results in the following generalization error

$\displaystyle E_{g}\sim$
$\displaystyle\frac{1}{M}|{\bm{w}}^{*}|^{2}\left[(1-\eta)^{2}+\eta^{2}\frac{(1-\eta)^{2}}{\eta-\eta^{2}}\right]\sim\frac{1}{M}|{\bm{w}}^{*}|^{2}(1-\eta).$
(70)

We see that as $\eta=\frac{N_{\mathcal{H}}}{M}\to 1$, the asymptotic error
converges to zero since all information in the original features is preserved.

####  D.5.2 Simplified Isotropic Feature Noise

We can simplify the above expressions somewhat in the case where
$\sigma^{2}=1$ and $\bm{\Sigma}_{\epsilon}=\sigma^{2}_{\epsilon}{\bm{I}}$. In
this case, the order parameters become

$\displaystyle\eta v$
$\displaystyle=\eta(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{-1}\implies
v=\frac{1}{1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}}$
$\displaystyle\implies\eta\hat{v}$ $\displaystyle=\hat{q}\
\text{tr}\left[{\bm{I}}+\frac{\hat{q}}{1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}}\bm{\Sigma}_{M}\right]^{-1}\bm{\Sigma}_{M}=\hat{q}(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})\
\text{tr}[(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}){\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-1}\bm{\Sigma}_{M}$
$\displaystyle q$
$\displaystyle=\text{tr}\left[(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}){\bm{I}}+\hat{q}\bm{\Sigma}_{M}\right]^{-1}\bm{\Sigma}_{M}+\frac{\eta\sigma_{\epsilon}^{2}}{1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}}.$
(71)

Letting
${\bm{G}}=[(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q}){\bm{I}}+\hat{q}\bm{\Sigma}_{M}]^{-1}$,
the source derivatives have the form

$\displaystyle\partial_{J}\hat{q}$
$\displaystyle=1-\frac{\alpha}{(\lambda+q)^{2}}\partial_{J}q$ (72)
$\displaystyle=1+\frac{\alpha}{(\lambda+q)^{2}}\left[\text{tr}{\bm{G}}^{2}\bm{\Sigma}[\partial_{J}\hat{v}{\bm{I}}+\partial_{J}\hat{q}\bm{\Sigma}_{M}]+\frac{\eta\sigma^{2}_{\epsilon}}{(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{2}}(\partial_{J}\hat{v}+\sigma^{2}_{\epsilon}\partial_{J}\hat{q})\right],$
$\displaystyle\eta\partial_{J}\hat{v}$
$\displaystyle=((1+\hat{v}+2\sigma^{2}_{\epsilon}\hat{q})\partial_{J}\hat{q}+\hat{q}\partial_{J}\hat{v})\text{tr}{\bm{G}}\bm{\Sigma}_{M}$
(73)
$\displaystyle\quad-\hat{q}(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})\text{tr}{\bm{G}}^{2}\bm{\Sigma}_{M}[(\partial_{J}\hat{v}+\sigma^{2}_{\epsilon}\partial_{J}\hat{q}){\bm{I}}+\partial_{J}\hat{q}\bm{\Sigma}_{M}]$
$\displaystyle=(\partial_{J}\hat{q})(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{2}\text{tr}{\bm{G}}^{2}\bm{\Sigma}+(\partial_{J}\hat{v}+\sigma^{2}_{\epsilon}\partial_{J}\hat{q})\hat{q}^{2}\text{tr}{\bm{G}}^{2}\bm{\Sigma}^{2}.$
(74)

This is a $2\times 2$ linear system

$\displaystyle\begin{bmatrix}1-\frac{\alpha}{(\lambda+q)^{2}}[\text{tr}{\bm{G}}^{2}\bm{\Sigma}^{2}+\frac{\eta\sigma^{4}_{\epsilon}}{(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{2}}]&-\frac{\alpha}{(\lambda+q)^{2}}[\text{tr}{\bm{G}}^{2}\bm{\Sigma}+\frac{\eta\sigma^{2}_{\epsilon}}{(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{2}}]\\\
-(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})^{2}\text{tr}{\bm{G}}^{2}\bm{\Sigma}_{M}-\sigma^{2}_{\epsilon}\hat{q}^{2}\text{tr}{\bm{G}}^{2}\bm{\Sigma}_{M}^{2}&\eta-\hat{q}^{2}\text{tr}{\bm{G}}^{2}\bm{\Sigma}^{2}_{s}\end{bmatrix}\begin{bmatrix}\partial_{J}\hat{q}\\\
\partial_{J}\hat{v}\end{bmatrix}=\begin{bmatrix}1\\\ 0\end{bmatrix}.$

For each $\alpha$, we can solve for $\partial_{J}\hat{q}$ and
$\partial_{J}\hat{v}$ to get the final generalization error with the formula

$\displaystyle E_{g}$
$\displaystyle=\partial_{J}\frac{1}{M}{\bm{w}}^{*}\left[(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})\hat{q}\bm{\Sigma}{\bm{G}}\right]{\bm{w}}^{*}$
$\displaystyle=\frac{1}{M}{\bm{w}}^{*}[\partial_{J}(\hat{q}+\hat{q}\hat{v}+\sigma^{2}_{\epsilon}\hat{q}^{2})\bm{\Sigma}{\bm{G}}-(1+\hat{v}+\sigma^{2}_{\epsilon}\hat{q})\hat{q}\bm{\Sigma}{\bm{G}}^{2}(\partial\hat{v}{\bm{I}}+\sigma^{2}_{\epsilon}\partial\hat{q}{\bm{I}}+\partial\hat{q}\bm{\Sigma}_{M})]{\bm{w}}^{*}.$
(75)

An example of these solutions can be found in Figure 16, where we show good
agreement between theory and experiment.

##  Appendix E ResNet on CIFAR Experiments

(a) Generalization MSE (b) Accuracy

(c) var/$E_{g}$ (d) var/$E_{g}$

Figure 8: A Wide ResNet Zagoruyko & Komodakis (2017) trained on a superclassed
CIFAR task comparing animate vs inanimate objects. Each learning curve is
averaged over $5$ different samples of the train set, yielding the means and
error bars shown in the figures. a) Generalization error $E_{g}$. The dashed
lines are the error of a 20-fold ensemble over different values of $\alpha$.
Across all $P$, lazy networks attain worse generalization error. As with the
MLP task, the best performing networks are ensembles of rich networks. b) The
accuracy also has the same trend: richer networks perform better and
ensembling lazy networks helps them more. c) Once $P$ is large enough, lazier
networks tend to benefit more from ensembling. d) Very lazy networks
transition to variance limited behavior earlier. For ResNets on this task, we
see that rich, feature learning networks eventually begin reducing their
variance on this task. Further details of the experiment are given in section
A.1.

##  Appendix F Additional Experiments

(a) $E_{g}$ for eNTK0 $k=2$ (b) $E_{g}$ for eNTK0 $k=3$ (c) $E_{g}$ for eNTK0
$k=4$ (d) $E_{g}$ for NN $k=2$ (e) $E_{g}$ for NN $k=3$ (f) $E_{g}$ for NN
$k=4$

Figure 9: Phase plots of $\log_{10}E_{g}$ for initial eNTKs (top) and neural
networks (bottom). The large $\alpha$ behavior of the neural network
generalization matches the generalization of the corresponding eNTK0 . As a
sanity check, the eNTK0 generalization error is independent of re-scaling the
network initialization because of the homogeneity of the ReLU network output.

Figure 10: A fine-grained view of the generalization error across different
datasets and ensembles. Solid curves are depth $3$ neural networks, dashed
curves are the infinite width NTK (which only has variance over datasets).
Each color is a set of networks trained on the same dataset but different
initializations. Different colors correspond to different datasets indexed by
$d\in\\{0,\dots,9\\}$.

(a) $E_{g}$ for $k=1,P=600$ (b) $E_{g}$ for $k=1,P=10000$

(c) $E_{g}$ for $k=3,P=600$ (d) $E_{g}$ for $k=3,P=10000$

Figure 11: Phase plots of $\log_{10}E_{g}$ for neural networks in the
$N$-$\alpha$ plane. We plot these at different train set sizes $P$ and
different tasks $k$. The colors are fixed to match across networks trained on
the same task. Figure 12: Empirical plot of of the scaling of the variance of
the eNTK0 with $N$ with variance taken over 10 initializations and averaged 10
different datasets.

(a) $E_{g}$ for $L=2$ (b) $E_{g}$ for $L=3$ (c) $E_{g}$ for $L=4$

(d) Ensembled $E_{g}$ for $L=2$ (e) Ensembled $E_{g}$ for $L=3$ (f) Ensembled
$E_{g}$ for $L=4$

(g) $\mathrm{var}/E_{g}$ for $L=2$ (h) $\mathrm{var}/E_{g}$ for $L=3$ (i)
$\mathrm{var}/E_{g}$ for $L=4$

Figure 13: Sweep over depth $L=\\{2,3,4\\}$. Deeper networks in the rich
regime can more easily outperform the infinite width network for a larger
range of $P$. Also, for larger $L$ it is easier to deviate from the lazy
regime at a given $\alpha$. By contrast, on this task the shallower NTK∞
outperforms deeper NTK∞ s. As before, ensembled lazy networks approach NTK∞
and the variance rises with $P$.

(a) $E_{g}$ for $L=2$ (b) $E_{g}$ for $L=3$ (c) $E_{g}$ for $L=4$

(d) Ensembled $E_{g}$ for $L=2$ (e) Ensembled $E_{g}$ for $L=3$ (f) Ensembled
$E_{g}$ for $L=4$

(g) $\mathrm{var}/E_{g}$ for $L=2$ (h) $\mathrm{var}/E_{g}$ for $L=3$ (i)
$\mathrm{var}/E_{g}$ for $L=4$

Figure 14: Sweep over input dimension $D=\\{5,25,50\\}$. At larger input
dimensions rich networks can more easily outperform NTK∞ . This is a
consequence of the task depending on the low-dimensional projection
$\bm{\beta}\cdot\bm{x}$.

(a) $E_{g}$ for centered vs uncentered (b) Effect of ensembling (c) $E_{g}$
uncentered color plot

Figure 15: a) $E_{g}$ for the centered predictor
$\tilde{f}_{\theta}(\bm{x})-\tilde{f}_{\theta_{0}}(\bm{x})$ (solid) compared
to the generalization of the uncentered predictor $\tilde{f}_{\theta}(\bm{x})$
(dashed). At small $\alpha$, the difference is negligible, while at large
$\alpha$ the uncentered predictor does worse and does not approach eNTK0 . The
worse generalization can be understood as $\tilde{f}_{\theta_{0}}(\bm{x})$
effectively adding an initialization-dependent noise to the target $\bm{y}$.
b) The effect of ensembling becomes less beneficial for uncentered lazy
networks. c) Color plot of $E_{g}$. The lazy regime is different from the
eNTK0 generalization (c.f. Figure 9).

(a) Gaussian ${\bm{A}}$, vary $N_{\mathcal{H}}$ (b) $N_{\mathcal{H}}=750$,
vary $\sigma^{2}_{\epsilon}$ (c) $N_{\mathcal{H}}=1000$, Vary $k$ (d)
$N_{\mathcal{H}}=1000$, Vary $\lambda$

Figure 16: Verification of Gaussian ${\bm{A}}$ model. Solid lines are theory
and dots are experiments. (a) The effect of changing the student’s RKHS
dimension $N_{\mathcal{H}}$. Double descent overfitting peaks occur at
$P=N_{\mathcal{H}}$ (b) The effect of additive noise in the student features
$\bm{\Sigma}_{\epsilon}=\sigma_{\epsilon}^{2}\bm{\Sigma}_{M}$. (c) Learning
curves for fitting the $k$-th eigenfunction. All mode errors exhibit a double
descent peak at $P={N}_{\mathcal{H}}$ regardless of the task. (d)
Regularization can prevent the overfitting peak.

(a) $E_{g}$ for mixed mode task (b) $\mathrm{Var}\,\hat{y}$ for mixed mode
task

Figure 17: Width 500 depth 3 MLP learning a $D=25$ mixture of a linear and
cubic polynomials. a) Generalization error of NTK∞ (solid black) and MLP
(solid colored lines) on mixed mode task. The dashed lines are convex
combinations of the generalization curves for the pure mode $k=1,k=3$ tasks.
For the NTK∞ , the generalization curves sum to give the mixed mode curve, as
observed in Bordelon et al. (2020). We see that this also holds for the eNTK0
for sufficiently lazy networks, as predicted by the simple renadom feature
model considred in section 4 of this paper. b) The variance curves for the
same task. Again, for sufficiently lazy networks the variance is a sum of the
variances of the individual pure mode tasks, as predicted by our random
feature model.

(a) $k=3,N=177,\alpha=0.1$ (b) $k=3,N=177,\alpha=20$ (c)
$k=3,N=1000,\alpha=0.1$ (d) $k=3,N=1000,\alpha=20$

Figure 18: Comparison of the neural network predictor ($x$-axis) to eNTK0 ,
eNTKf (blue, orange respectively, $y$-axis) across several training dataset
sizes. a) $N=177,\alpha=0.1$, b) $N=177,\alpha=20$, c) $N=1000,\alpha=0.1$ d)
$N=1000,\alpha=20$. All networks are depth 3. Note how in the rich regime
there is a much stronger distinction between the eNTK0 and the neural network.
In all regimes, eNTKf matches the NN output.

(a) $k=1,N=177,\alpha=0.1$ (b) $k=1,N=177,\alpha=20$ (c)
$k=1,N=1000,\alpha=0.1$ (d) $k=1,N=1000,\alpha=20$

Figure 19: The same as figure 18 but for fitting a linear $k=1$ mode. Because
the task is simpler, it doesn’t require as large of an $\alpha$ to enter the
lazy regime.