\section{Introduction}

Neural networks perform better with more  resources (data, model size, training time), but different tasks exhibit qualitatively different dependencies of performance on resources.
In particular, while many learning tasks exhibit continuous improvement in performance with increasing resources, other cases show \emph{discontinuous} improvement, where a capability emerges at a certain threshold.
Through a statistical lens, it is well-understood that larger models, trained with more data, can fit more complex and expressive functions. However, far less is known about the analogous \emph{computational} question: how does the scaling of these resources influence the success of gradient-based optimization for training these models?

These emergent \emph{phase transitions} cannot be explained via statistical capacity alone. In many cases we see a phase transition even when the amount of data remains fixed, with only model size or training time increasing.
A timely example is the emergence of reasoning and few-shot learning capabilities when scaling up language models \citep{radford2019language,brown2020language,chowdhery2022palm,hoffmann2022training}. \citet{power2022grokking} give examples exhibiting discontinuous improvements in population accuracy (``grokking'') when running time increases, while data and model size remain fixed.

In this work, we analyze the computational aspect of scaling in deep learning in a simple synthetic setting which already exhibits discontinuous improvements. Specifically, we consider the \emph{sparse parity problem}--- where the label is the parity (XOR) of $k \ll n$ bits in a random length-$n$ binary string. This is a canonical problem which is computationally difficult for a range of algorithms, including gradient-based \citep{kearns1998efficient} and streaming \citep{kol2017time} algorithms. We focus on analyzing the resource measure of \emph{training time}, and demonstrate that the loss curves for sparse parities display a phase transition across a variety of architectures and hyperparameters (see Figure~\ref{fig:intro}, left). A striking observation is that SGD robustly finds the sparse subset (and hence, effectively, reaches 0 error) with a variety of activation functions and initialization schemes, even with \emph{no} over-parameterization. This robust convergence is the starting point for our investigation.

\begin{figure}
    \centering
    \includegraphics[width=0.62\linewidth]{figures/train_curves.pdf}
    \includegraphics[width=0.373\linewidth]{figures/scaling_curves.pdf}
    \caption{Main empirical findings at a glance. A variety of neural networks, with standard training and initialization, can solve the $(n,k)$-parity learning problem, with a number of iterations scaling as $n^{O(k)}$. \emph{Left:} Training curves under various algorithmic choices (architecture, batch size, learning rate) on the $(n=50, k=3)$-parity problem. \emph{Right:} Median convergence times for small $(n,k)$.}
    \label{fig:intro}
\end{figure}

Perhaps the most natural hypothesis to explain SGD's success in learning parities would be that it simply ``stumbles in the dark'', essentially performing random search for the unknown target (e.g. via stochastic gradient Langevin dynamics). If that were the case, we might expect to observe a convergence time of $2^{\Omega(n)}$, like a naive search over parameters or subsets of indices.
However, Figure~\ref{fig:intro} \emph{(right)}, already provides some evidence against this ``random search'' hypothesis: the convergence time is closer to an $n^{O(k)}$ scaling. Notably, such a convergence rate implies that SGD is closer to achieving the \emph{optimal computation time} among a natural class of algorithms, namely SQ (statistical query) algorithms.

In this work, through an extensive empirical analysis of the scaling behavior of a variety of models, as well as theoretical analysis, we give strong evidence against the ``stumbling in the dark'' viewpoint. Instead, there is a \emph{hidden progress measure} under which SGD is steadily improving. Furthermore, and perhaps surprisingly, we show that SGD achieves a computational runtime much closer to the optimal SQ lower bound than simply doing (non-sparse) parameter search. More generally, our theoretical and empirical investigations reveal a number of notable phenomena regarding the dependence of SGD's performance on resources, and we provide further theoretical and empirical evidence of phase transitions with data, model size, and training time.

\subsection{Our contributions}

\paragraph{SGD efficiently learns sparse parities, in theory and practice.}
It is known from SQ lower bounds that (noisy) gradient descent on \emph{any} architecture requires at least $n^{\Omega(k)}$ computational steps to learn $k$-sparse $n$-dimensional parities (for background, see Appendix~\ref{app:preliminaries}). However, with standard architectures and initialization, which do not explicitly encode a sparsity prior, one may expect SGD's performance to be much worse, on the order of $2^{\Omega(n)}$. We give extensive empirical evidence that this is not the case.

\begin{empfind}
For all small instances ($n \leq 30, k \leq 4$) of the sparse parity problem, architectures $\mathcal{A} \in \{$2-layer MLPs, Transformers\footnote{With a smaller range of hyperparameters.}, sinusoidal (and other non-standard) neurons, PolyNets\footnote{A non-standard architecture introduced in this work; see Section~\ref{sec:empirics} for the definition.}$\}$, initializations in $\{$uniform, Gaussian, Bernoulli$\}$, and batch sizes $1 \leq B \leq 1024$, SGD on $\mathcal{A}$ solves the $(n,k)$-sparse parity problem (w.p. $\geq 0.2$) within at most $c \cdot n^{\alpha k}$ steps, for small architecture-dependent constants $c, \alpha$.
\end{empfind}

For a subset of these architectures, we performed more extensive training, to show scaling behaviors of the computation time;\footnote{While our focus is on the performance as a function of \emph{training time}, we also performed some experiments on performance as a function of \emph{model size}, see Section~\ref{sec:discussion}.} see Figure~\ref{fig:intro} \emph{(right)} and Figure~\ref{fig:extra-scaling-curves} in the appendix.

The above results indicate that SGD succeeds at solving parities much faster than the $2^{\Omega(n)}$ steps that would be required by random search without a sparse prior. This suggests that despite flat loss and accuracy curves before the phase transition, SGD makes progress ``under the hood''. We show that this is indeed the case by coming up with a \emph{progress measure} that continuously improves throughout training. We can explain the loss and accuracy curves using this progress measure, by showing that they typically rise sharply once the measure passes a certain threshold.

\paragraph{Theoretical Analysis.} Our empirical results suggest that, in a number of computational steps matching the SQ limit, SGD is able to solve the parity problem and identify the influential coordinates, without an explicit sparse prior. We give a theoretical analysis which validates this claim.
\begin{inftheorem}
On 2-layer MLPs of width $2^{\Theta(k)}$, and with batch size $n^{O(k)}$, SGD converges with high probability to a solution with at most $\epsilon$ error on the $(n,k)$-parity problem in at most $2^{O(k)} \cdot \mathrm{poly}(1/\epsilon)$ iterations.
\end{inftheorem}

We also present a stronger analysis for an idealized architecture (which we call the \emph{disjoint-PolyNet}), which allows for any batch size, and captures the phase transitions observed in the error curves.

\begin{inftheorem}
On disjoint-PolyNets, SGD (with any batch size $B \geq 1$) converges with high probability to a solution with at most $\epsilon$ error on the $(n,k)$-parity problem in at most $n^{O(k)} \cdot \log(1/\epsilon)$ iterations. Continuous-time gradient flow exhibits a phase transition: it spends a $1 - o(1)$ fraction of its time before convergence with error $\geq 49\%$.
\end{inftheorem}

In addition to refuting the alternative ``random search'' hypothesis, our work also poses a counterexample to other models for the computational mechanisms of deep learning; for instance, it provides a setting where deep neural nets successfully learn a concept to 100\% accuracy while the corresponding Neural Tangent Kernel (NTK) only achieves trivial performance, hence showing the importance of \emph{feature learning}.
We also construct a counterexample to the \emph{``deep only works if shallow is good''} principle of \cite{malach2019deeper}, demonstrating a case where a deep network can get near-perfect accuracy even when greedy layerwise training (e.g. \citet{belilovsky2019greedy}) cannot beat trivial performance.
By providing positive theory and empirics which elude these simplified explanations of SGD, we hope to point the way to a more complete understanding of learning dynamics in the challenging cases where no apparent progress is made for extended periods of time.

\subsection{Related work} \label{sec:related}
We present the most directly related work below, on feature learning, learning parities with neural nets, and the statistical mechanics of machine learning. A broader discussion of related work can be found in Appendix \ref{app:relwork}.

\paragraph{SGD and feature learning.}
Theoretical analysis of gradient descent over neural networks is notoriously hard, due to the non-convex nature of the optimization problem. That said, it has been established that in some settings, the dynamics of GD keep the weights close to their initialization, thus behaving like convex optimization over the Neural Tangent Kernel (see, for example, \cite{jacot2018neural, allen2019convergence, du2018gradient}). In contrast, it has been shown that in various tasks, moving away from the fixed features of the NTK is essential for the success of neural networks trained with GD (for example \cite{yehudai2019power, allen2019can, wei2019regularization} and the review in \cite{malach2021quantifying}). These results demonstrate that feature learning is an important part of the GD optimization process. Our work also focuses on a setting where feature learning is essential for the target task. In our theoretical analysis, we show that the initial population gradient encodes the relevant features for the problem. The importance of the first gradient step for feature learning has been recently studied in \citep{ba2022high}.

\paragraph{Learning parities with neural networks.}
The problem of learning parities using neural networks has been investigated in prior works from various perspectives. It has been demonstrated that parities are hard for gradient-based algorithms, using similar arguments as in the SQ analysis \citep{shalev2017failures, abbe2020poly}. One possible approach for overcoming the computational hardness is to make favorable assumptions on the input distribution. Indeed, recent works show that under various assumptions on the input distribution, neural networks can be efficiently trained to learn parities (XORs) \citep{daniely2020learning,shi2021theoretical,frei2022random,malach2021quantifying}. In contrast to these results, this work takes the approach of intentionally focusing on a hard benchmark task, without assuming that the distribution has some favorable (namely, non-uniform) structure. This setting allows us to probe the performance of deep learning at a known computational limit. Notably, the work of \cite{andoni2014learning} provides analysis for learning polynomials (and in particular, parities) under the uniform distribution. However, their main results require a network of size $n^{O(k)}$ (i.e., extremely overparameterized network), and provides only partial theoretical and empirical evidence for the success of smaller networks.
Studying a related subject, some works have shown that neural networks display a spectral bias, learning to fit low-frequency coefficients before high-frequency ones \citep{rahaman2019spectral,cao2019towards}.

\paragraph{Statistical mechanics of machine learning.}
An extensive body of work originating in the statistical physics community has studied phase transitions in the learning curves of neural networks \citep{gardner1989three, watkin1993statistical, engel2001statistical}. These works typically focus on student-teacher learning in the ``thermodynamic limit’’, in which the number of training examples is $\alpha$ times larger than the input dimension and both are taken to infinity. One of the classic toy architectures analyzed in this literature is the \emph{parity machine} \citep{mitchison1989bounds, hansel1992memorization, opper1994learning, kabashima1994perfect, simonetti1996line}. In our work, we introduce PolyNets, a variant of parity machines in which the output is real-valued rather than $\pm 1$; and we theoretically analyze disjoint-Polynets, which are the real-output analogue of the oft-considered parity machines with tree architecture. While much of the statistical mechanics of ML literature focuses on an idealized training limit in which the weights reach a Gibbs distribution equilibrium, there is a strand of the literature that aims to characterize the trajectory of SGD training in the high-dimensional limit with constant-sized sets of ordinary differential equations \citep{saad1995line, saad1995dynamics, goldt2019dynamics}. These papers discuss cases, including problems that share aspects with 2-sparse parities \cite{refinetti2021classifying}, where the network gets stuck in (and then escapes from) a plateau of suboptimal generalization error. Recently, \citet{arous2021online} studied (for rank-one parameter estimation problems) the relative amount of time spent by SGD in an initial high-error ``search'' phase versus a final ``descent'' phase, which is reminiscent of the framing of Theorem~\ref{thm:fraction}. However, to our knowledge prior work has not shown $k$-sparse parities can be learned with a number of iterations that nearly matches known lower bounds, nor has it specifically studied phase transitions in $k$-sparse parity learning during gradient descent.
\section{Preliminaries}

We define necessary notation here; see Appendix~\ref{app:preliminaries} for more background and technical ingredients.

\paragraph{Sparse parities.} For integer $n\geq 1$ and non-empty $S \subseteq [n]$, the $(n,S)$-\emph{parity function} $\chi_S:\{ \pm 1 \}^n \rightarrow \{ \pm 1 \}$ is defined as $\chi_S(x) = \prod_{i\in S}x_i$. We define the $(n,S)$-\emph{parity distribution} $\mathcal{D}_S$ as the distribution over $(x,y=\chi_S(x))$\footnote{Our theoretical analyses and experiments can tolerate noisy parities, that is, random flipping of the label. For ease of presentation, we state the non-noisy setting.} where $x$ is drawn from the uniform distribution over $\{\pm 1\}^n$, which we denote by $\mathrm{Unif}(\{\pm 1\}^n)$. We define the $(n,k)$-\emph{parity learning problem} as the task of recovering $S$ using samples from $\mathcal{D}_S$, where $S$ is chosen at random in $\binom{[n]}{k}$. Statistically, it is possible to do so using $\Theta(\log \binom{n}{k}) \approx k \log n$ samples.
However, in the \emph{statistical query} (SQ) model \citep{kearns1998efficient}, (which has been shown to encapsulate gradient-based methods such as GD or SGD~\citep{abbe2021power}), this task requires $\Omega(n^k)$ queries (assuming constant level of noise).
While learning noiseless parities can be solved by Gaussian elimination using $O(n)$ samples, learning sparse \emph{noisy} parities, even at a very small noise level (i.e., $o(1)$ or $n^{-\delta}$) is believed to inherently require $n^{\Omega(k)}$ computational steps and samples, or exponential computation with $n^{o(k)}$ samples. This was first explicitly conjectured by \citet{alekhnovich2003more}, and has been the basis for several cryptographic schemes (e.g., \citep{ishai2008cryptography, applebaum2009fast, applebaum2010public, bogdanov2019xor}).

\paragraph{Notation for neural networks and training.} Our main results are presented in the online learning setting, with a stream of i.i.d. batches of examples. At each iteration $t = 1, \ldots, T$, a learning algorithm receives a batch of $B$ examples $\{(x_{t,i},y_{t,i})\}_{i = 1}^{B}$ drawn i.i.d. from $\mathcal{D}_S$, then outputs a classifier $\hat y_t : \{\pm 1\}^n \rightarrow \{\pm 1\}$. We say that the algorithm solves the parity task in $t$ steps with $\epsilon$ error, if with probability at least $1-\epsilon$ over both training and internal randomness, as well as $(x,y) \sim \mathcal{D}_S$, $\hat y_t(x)=y$. We will focus on the case that $\hat y_t = \mathop{\mathrm{sign}}(f(x;\theta_t))$ for some parameters $\theta_t$ in a continuous domain $\Theta$ and for a continuous function   $f : \{\pm 1\}^n \times \Theta \rightarrow \RR$.\footnote{When $f(x;\theta) = 0$ in practice (e.g. with sign initialization), we break the tie arbitrarily. We ensure in the theoretical analysis that this does not happen.}
A ubiquitous online learning algorithm is gradient descent (GD). For a choice of loss function $\ell : \{\pm 1\} \times \RR \rightarrow \RR$, initialization $\theta_0$ (that are chosen from some distribution), learning rate schedule $\{\eta_t\}_{t=1}^T \subseteq \RR$ and weight-decay schedule $\{\lambda_t\}_{t=1}^T \subseteq \RR$, GD refers the standard iterative update rule using the regularized, empirical loss function, which is a function of \emph{architecture}  $f$.
The learning rate $\eta_t$ can also be a vector (e.g., allowing different rate schedules for different layers).

\section{Empirical findings}
\label{sec:empirics}

\begin{figure}
    \centering
    \includegraphics[width=0.34\linewidth]{figures/hists.pdf}
    \hspace{-5mm}
    \includegraphics[width=0.34\linewidth]{figures/clusters.pdf}
    \includegraphics[width=0.31\linewidth]{figures/mlp_scaling_large.pdf}
    \caption{Black-box observations about the training dynamics. \emph{Left:} Histograms of convergence times over $10^6$ random trials, with heavy upper tails but no observed successes near $t = 0$ (unlike random search). \emph{Center:} Optimization path (and convergence time) depend heavily on initialization, not the randomness of SGD; $B = 128, \eta = 0.01$ shown here. \emph{Right:} The power-law exponent ($\alpha$ such that $t_c \propto n^{\alpha}$) eventually degrades on larger problem instances. }
    \label{fig:scaling-laws}
\end{figure}

\subsection{SGD on neural networks learns sparse parities}
\label{subsec:empirics-main}

The central phenomenon of study in this work is the empirical observation that neural networks, with standard initialization and training, can solve the $(n,k)$-parity problem in a number of iterations scaling as $n^{O(k)}$ on small instances. We observed robust positive results for randomly-initialized SGD on the following architectures, indexed by Roman numerals:
\begin{itemize}[leftmargin=1em]
    \item \emph{2-layer MLPs:} ReLU $(\sigma(z) = (z)_+)$ or polynomial ($\sigma(z) = z^k$) activation, in a wide variety of width regimes $r \geq k$. Settings (i), (ii), (iii) (resp. (iv), (v), (vi)) use $r = \{10, 100, 1000\}$ ReLU (resp. polynomial) activations. We also consider $r = k$ (exceptional settings (*i), (*ii) ), the minimum width for representing a $k$-wise parity for both activations.
    \item \emph{Single neurons:} Next, we consider non-standard activation functions $\sigma$ which allow a one-neuron architecture $f(x; w) = \sigma(w^\top x)$ to realize $k$-wise parities. The constructions stem from letting $w^* = \sum_{i \in S} e_i$, and constructing $\sigma(\cdot)$ to interpolate (the appropriate scaling of) $\frac{k - {w^*}^\top x}{2} \text{ mod } 2$ with a piecewise linear \emph{$k$-zigzag} activation (vii), or a degree-$k$ polynomial (viii). Going a step further, a single \emph{$\infty$-zigzag} (ix) or \emph{sinusoidal} (x) neuron can represent \emph{all} $k$-wise parities. We also removed the second trainable layer (setting $u = 1$), obtaining settings (xi), (xii), (xiii), (xiv). We found that wider architectures with these activations also trained successfully.
    \item \emph{Transformers:} Motivated by recent theoretical and empirical work on the ability of self-attention to learn sparse functions and parities \citep{lu2021pretrained,edelman2021inductive,hahn2020theoretical}, we consider a simplified specialization of the Transformer architecture to this sequence classification problem.
    This is the less-robust setting (*iii); the architecture and optimizer are described in Appendix~\ref{subsubsec:parity-transformer}.
    \item \emph{PolyNets:} Our final setting (xv) is the PolyNet, a slightly modified version of the parity machine architecture. Parity machines have been studied extensively in the statistical mechanics of ML literature (see the related work section) as well as in a line of work on `neural cryptography' \citep{rosen2002mutual}. A parity machine outputs the sign of the product of $k$ linear functions of the input. A PolyNet simply outputs the product itself. Both architectures can clearly realize $k$-sparse parities. The PolyNet architecture was originally motivated by the search for an idealized setting where an end-to-end optimization trajectory analysis is tractable (see Section~\ref{subsec:theory-theorems}); we found in these experiments that this architecture trains very stably and sample-efficiently.
\end{itemize}

\paragraph{Robust space of positive results.}
All of the networks listed above were observed to successfully learn sparse parities in a variety of settings. We summarize our findings as follows:
for all combinations of $n \in \{10, 20, 30\}$, $k \in \{2, 3, 4\}$, batch sizes $B \in \{1, 2, 4, \ldots, 1024\}$, initializations $\{$uniform, Gaussian, Bernoulli$\}$, loss functions $\{$hinge, square, cross entropy$\}$, and architecture configurations $\{ \text{(i)}, \text{(ii)}, \ldots, \text{(xv)} \}$, SGD solved the parity problem (with $100\%$ accuracy, validated on a batch of $2^{13}$ samples) in at least $20\%$ of $25$ random trials, for at least one choice of learning rate $\eta \in \{0.001, 0.01, 0.1, 1\}$. The models converged in $t_c \leq c \cdot n^{\alpha k} \leq 10^5$ steps, for small architecture-dependent constants $c, \alpha$, shown in Appendix~\ref{sec:appendix-experiments}. Figure~\ref{fig:intro} \emph{(left)} shows some representative training curves.

\paragraph{Less robust configurations.} Settings (*i) and (*ii), where the MLP just barely represents a $k$-sparse parity, and the Transformer setting (*iii), are less robust to small batch sizes. In these settings, the same positive results as above only held for sufficiently large batch sizes: $B \geq 16$. Also, setting (*iii) used the Adam optimizer; see Appendix~\ref{subsubsec:parity-transformer} for details.

\paragraph{Phase transitions in training curves.} For almost all of the architectures, we find that that the training curves exhibit phase transitions in terms of running time (and thus, in the online learning setting, dataset size as well): long durations of seemingly no progress, followed by periods of rapid decrease in the validation error. Strikingly, for architectures (v) and (vi), this plateau is absent: the error in the initial phase appears to decrease with a linear slope. See Appendix~\ref{subsec:appendix-poly-no-plateau} for more plots.

\begin{figure}
    \centering
    \includegraphics[width=0.4\linewidth]{figures/sin_weights.pdf}
    \hspace{-4mm}
    \includegraphics[width=0.37\linewidth]{figures/sin_weights_largebatch.pdf}
    \includegraphics[width=0.232\linewidth]{figures/sin_progress.pdf}
    \caption{Hidden progress when learning parities with neural networks. \emph{Left, center:} Black-box losses and accuracies exhibit a long plateau and sharp phase transition (top), hiding gradual progress in the SGD iterates (bottom). \emph{Right:} A hidden progress measure which distinguishes gradual feature amplification (top) from training on noise (bottom). }
    \label{fig:hidden-progress}
\end{figure}

\subsection{Random search or hidden progress?}
\label{subsec:empirics-random-or-progress}

The remainder of this paper seeks to answer the question:
\emph{``By what mechanism does deep learning solve these emblematic computationally-hard optimization problems?''}

A natural hypothesis would be that SGD somehow implicitly performs Monte Carlo random search, ``bouncing around'' the loss landscape in the absence of a useful gradient signal. Upon closer inspection, several empirical observations clash with this hypothesis:

\begin{itemize}[leftmargin=*]
    \item \emph{Scaling of convergence times:} Without an explicit sparsity prior in the architecture or initialization, it is unclear how to account for the $n^{\Theta(k)}$ runtimes observed in experiments. The initializations, which certainly do not prefer sparse functions\footnote{Indeed, under all standard architectures and initialization, the probability that a random network is $\Omega(1)$-correlated with a sparse parity would be $2^{-\Omega(n)}$, since with that probability $1-o(1)$ of the total influence would be accounted by the $n-k$ irrelevant features.}, are close to the correct solutions with probability $2^{-\Omega(n)} \ll n^{-k}$. On larger instances $(n,k)$, the power-law exponents worsen; see Figure~\ref{fig:scaling-laws} \emph{(right)}, and the discussion in Appendix~\ref{subsec:building-blocks}.
    \item \emph{No early convergence:} Over a large number of random trials, no copies of this randomized algorithm get ``lucky'' (i.e. solve the problem in significantly fewer than the median number of iterations); see Figure~\ref{fig:scaling-laws} \emph{(left)}. The success times of random exhaustive search would be distributed as $\mathrm{Geom}(1/\binom{n}{k})$, whose probability mass is highest at $t=0$ and decreases monotonically with $t$.
    \item \emph{Sensitivity to initialization, not stochastic batches:} Running these training setups over multiple stochastic batches from a common initialization, we find that loss curves and convergence times are highly correlated with the architecture's random initialization; see Figure~\ref{fig:scaling-laws} \emph{(center)}.
\end{itemize}

Even these %
observations, which do not probe the internal state of the algorithm, suggest that
exhaustive search is an insufficient picture of the training dynamics, and a different mechanism is at play.
\section{Theoretical analyses}
\label{sec:theory}

\subsection{Provable emergence of the parity indices in high-precision gradients}
\label{subsec:theory-theorems}

We now provide a theoretical account for the success of SGD in solving the $(n,k)$-parity problem. Our main theoretical observation is that, in many cases, the \emph{population} gradient of the weights at initialization contains enough ``information'' for solving the parity problem. That is, given an accurate enough estimate of the initial gradient (by e.g. computing the gradient over a large enough batch size), the relevant subset $S$ can be found.

As a warm-up example, consider training a single ReLU neuron $f(x;w) = \sigma(w^\top x)$ w.r.t. the correlation loss $\ell(y,\hat{y}) = -y\hat{y}$, from the initialization $w = [1,\dots, 1]$. While a single neuron cannot \emph{express} the parity, we observe that the population gradient can indicate what the correct subset is:
$\mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S}\left[\nabla_{w_i} \ell(y,f(x;w))\right]
= \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S}\left[-y x_i \, \sigma'(w^\top x)\right]$ which corresponds to either the order-($k-1$) Fourier coefficient $S\setminus \{ i \}$ of the function $x \mapsto \sigma'(w^\top x)$  (if $i\in S$ is a relevant coordinate) or the order-$(k+1)$ coefficient $S \cup \{i\}$ (if $i\not\in S$ is irrelevant). When $\sigma$ is the ReLU function and $w = [1,\dots, 1]$, $\sigma'(w^\top x) = \frac{\mathop{\mathrm{sign}}\left(\sum_i x_i\right)+1}{2}$ is just a shifted majority function of $x$. The Fourier spectrum of majority is well-understood: for even $k$, there is a gap between these Fourier coefficients that is detectable using $n^{O(k)}$ samples.

This analysis can be further extended to a ReLU neuron initialized with weights $w \in \{\pm 1\}^n$ and a small bias. In fact, we can show that taking a single gradient step with large enough batch size on a ReLU network trained with the hinge-loss $\ell(y,\hat{y}) = \max \{1 - y\hat{y},0\}$, already finds features that can solve the parity problem. Using this, we show that training a ReLU MLP with SGD solves the parity problem:

\begin{theorem}
\label{thm:mlp-main}
Fix some $\epsilon \in (0,1)$, let $k$ be some even number, and assume that $n$ is some odd number satisfying $n \ge \Omega(k^4 \log(nk/\epsilon))$. There exist a symmetric random initialization scheme and learning rate and weight decay schedules s.t. for every $S \subseteq [n]$ of size $k$, training a ReLU MLP of size $r = O(2^k k \log(k/\epsilon))$ with batch size $B = n^{O(k)}\log^2(r/\epsilon)$ for $T = \mathrm{poly}(k,r,1/\epsilon)$ iterations on $\mathcal{D}_S$ w.r.t. the hinge loss, finds a network $f(x;\theta_t)$ with expected loss:
$
\mathop{\mathbb{E}} \left[\ell(f(x;\theta_t),y)\right] \le \epsilon
$,
where the expectation is over the randomness of initialization, training and sampling $(x,y) \sim \mathcal{D}_S$.
\end{theorem}

The analysis presented above shows that when the \emph{batch size} scales with $n^{O(k)}$, SGD over MLP solves the parity problem. However, in our experimental setting, the \emph{number of steps}, and not the batch size, scales with $n^{O(k)}$.
While we believe that running SGD with small batch size and small learning rate essentially amplifies the signal in the population gradient, behaving similarly to performing a large step over a large batch size, we do not have a complete analysis for training MLPs with SGD in the small batch size regime.
To complement the above result, in Section \ref{subsec:polynet} we analyze the trajectory for a variant of the PolyNet architecture trained with gradient flow.

We note that, while the above analysis applies for ReLU MLPs with a specific initialization scheme, a similar feature emergence phenomenon can be observed in a broader set of architectures and setting. Indeed, for feature emergence to occur, we only require that there is a ``gap'' between the relevant and irrelevant Fourier coefficients. Formally, denote $\hat{f}(S) := \mathop{\mathbb{E}} \left[f(x)\chi_S(x)\right]$ the Fourier coefficient of $f$ at $S$, and observe the following definition:

\begin{definition}[Fourier gap]
\label{def:fourier-gap}
For some function $f : \{\pm 1\}^n \to \RR$ and some subset $S$ of size $k$, we say that $f$ has a $\gamma$-\emph{Fourier gap} at $S$ if (1) for every $(k-1)$-element subset $S' \subseteq S$, it holds that $|\hat{f}(S')| \ge \gamma$, and (2) for every subset $S' \subseteq [n]$ of size $k+1$ it holds that $|\hat{f}(S')| \le \gamma/2$.
\end{definition}

Now, given a network architecture where some neuron has a $\gamma$-Fourier gap with respect to the target subset $S$, we can generalize the result of the ReLU neuron. That is, we show that the subset $S$ can be determined by observing an estimate of the population gradient at initialization:
\begin{proposition}
\label{prop:mlp-from-gap}
Let $\sigma$ be some activation function and let $\ell$ be some loss function. Denote $f(x;w) = \sigma(w^{\top} x)$. Fix some subset $S \subseteq [n]$.
Let $g \in \RR^n$ be an estimate of the population gradient such that
$
\left\|g-\mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S}\left[\nabla_w \ell(y,f(x;w))\right]\right\|_\infty \le \gamma/4.
$
Then, for every $w$ s.t. $\sigma'(w^\top x)$ has a $\gamma$-Fourier gap w.r.t. to $S$ and $\ell'(f(x;w),y) = -y$ for all $x$,
the target subset $S$ is detected by $g$, namely $S = \{i \in [n] ~:~|g_i| \ge 3\gamma/4\}$.
\end{proposition}

\paragraph{Comparison with NTK analysis.} In recent years, many theoretical works have studied the behavior of SGD on neural networks through the lens of the neural tangent kernel (NTK) \citep{jacot2018neural}. It is therefore important to highlight the fact that the NTK (or, in fact, any kernel) cannot solve the sparse parity problem. The following result shows that no kernel can achieve small error on the sparse parity problem, unless the size of the kernel is $n^{\Omega(k)}$:
\begin{theorem}
Let $\Psi : \{\pm 1\}^n \to \RR^D$ be some $D$-dimensional embedding with $\sup_x \|\Psi(x)\|_2 \le 1$, and let $R > 0$ be some number. If $D R^2 < \epsilon^2 \cdot \binom{n}{k}$, then there exists some $(n,k)$-parity distribution $\mathcal{D}_S$ s.t.
$$\inf_{\|w\| \le R} \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S} \left[\ell(\Psi(x)^\top w, y)\right] > 1-\epsilon.$$
\end{theorem}

The above is a corollary of the results in \citep{kamath2020approximate,malach2022hardness}.

\subsection{Disjoint-PolyNet: an idealized architecture for trajectory analyses}
\label{subsec:polynet}
In this section, we present an idealized architecture (a version of PolyNets (xv)) that exhibits similar behavior to MLPs (experimentally) and is technically easier to analyze. More specifically, we consider the \emph{disjoint-PolyNet} which takes a product over $k$ linear functions where the linear functions are computed on $k$ disjoint partitions of the input $P_1, \ldots, P_k$ with each $P_i = \{(i-1)n' + 1, \ldots, in'\}$ with $n' = n/k$\footnote{We assume for simplicity that $n$ is divisible by $k$.}, that is, $f(x;w_{1:k}) := \prod_{i=1}^k w_i^\top x_{P_i}.$ As noted in the related work section, this is equivalent to a tree parity machine but with real-valued rather than $\pm 1$ output.

In order for the disjoint-PolyNet to be able to express the class of parity problems, we assume that the set $S$ of size $k$ in the $(n,k)$-parity problem is selected such that exactly one index belongs to each disjoint partition, that is, for all $i \in [k]$, $S\cap P_i = 1$. We refer to this problem as the $(n,k)$-disjoint parity problem. Note that there are still $(n')^k = (n/k)^k$ different possibilities for set $S$ under this restriction. For fixed $k$, these represent a constant portion of the $\binom{n}{k} \approx (ne/k)^k$ (by Stirling's approximation) possibilities for $S$ in the general non-disjoint case.

Consider training a disjoint-PolyNet w.r.t. the correlation loss. WLOG, let $S=\{1,n'+1,\ldots,(k-1)n'+1\}$ and $e_1=(1,0,..0)$. The population gradient is non-zero at $i$ iff $i\in S$: $g_i(w_{1:k})= \mathop{\mathbb{E}}\left[\nabla_{w_i} \ell(f(x;w_{1:k}), y)\right] = - \mathop{\mathbb{E}} \left[ y \left(\prod_{j \ne i} w_j^\top x_{P_j}\right) x_{P_i} \right]=  -\left(\prod_{j \ne i} w_{j,1} \right)e_1$.%

Now we consider the gradient flow dynamics of disjoint-PolyNet, which provide a mathematically tractable case study for the trajectory of sparse parity learning. For each $i \in [k]$, in this section we treat $w_i$ as a function from $\RR_{\geq 0} \to \RR^{n'}$ which satisfies the following differential equation: $\dot{w_i} = -g_i(w_{1:k}(t))$. For clarity of exposition, assume all-ones initialization.\footnote{Results for $\pm 1$ initialization and Gaussian initialization are qualitatively similar and can be found in the Appendix.
} Then all of the relevant weights $\{w_{i,1}: i \in [k]\}$ follow the same trajectory, which we denote by $v: \RR_{\geq 0} \to \RR$.

By analyzing the resulting differential equations, we can formally exhibit ``phase transition''-like behavior in the fully deterministic gradient flow setting.

\begin{theorem}[Gradient flow on disjoint-PolyNets]\label{thm:fraction}
Suppose $k \geq 3$. Let $T(.49)$ be the time it takes for error to fall below $.49$, and let $T(0)$ be the time it takes to reach zero error. Then
\[
\frac{T(.49)}{T(0)} = 1 - O\left((n')^{1-k/2}\right).
\]
\end{theorem}

In other words, the network takes much longer to reach slightly better than trivial accuracy than it takes to go from slightly better than trivial to perfect accuracy.

We can also analyze the trajectory of disjoint-PolyNets trained with online SGD, confirming that a neural network trained with batch size 1 SGD can learn $k$-sparse parities within $n^{O(k)}$ iterations.

\begin{theorem}[SGD on disjoint-PolyNets]
Suppose we train the disjoint-PolyNet, initialized as above, with online SGD with any batch size $B \geq 1$. Then there exists an adaptive learning rate schedule such that for any $\epsilon > 0$, with probability .99, the error falls below $\epsilon$ within
\[\tilde{O}\left((n')^{(2k-1)} \log(1/\epsilon)\right)\]
steps.
\end{theorem}

Extended versions of these theorems, along with proofs, can be found in Appendix~\ref{app:polynet}.
\section{Hidden progress: discussion and additional experiments}
\label{sec:discussion}

\begin{figure}
    \centering
    \includegraphics[width=0.36\linewidth]{figures/width.pdf}
    \includegraphics[width=0.63\linewidth]{figures/grok.pdf}
    \caption{Parity as a sandbox for understanding the effects of model size and dataset size. \emph{Left:} Success times vs. network width $r$ on a fixed $(40,3)$-parity task: in accordance with the theory, parallelization experiences diminishing returns (unlike expected success times for random search, shown in green). Underparameterized models ($r = 1, 2$) were considered successful upon reaching $55\%$ accuracy. \emph{Right:} Training curves for an identical setup ($(50,3)$-parity task, architecture, and training algorithm), varying only the sample size $m$. The two center panels display \emph{``grokking''}: a large gap between the time to zero train error vs. zero test error. }
    \label{fig:extra-experiments}
\end{figure}

In this section, we advocate for sparse parities as an idealized testbed for understanding algorithms and phenomena in modern deep learning.
These are accompanied by experimental vignettes which are auxiliary to the core results from Section~\ref{sec:empirics}, with more systematic studies deferred to future work.
Details are given in Appendix~\ref{sec:experiment-details}.

\paragraph{A progress measure for parity.} To begin, the black-box experiments in Section~\ref{sec:empirics} suggest that random search is the incorrect model of SGD's behavior in this setting. Using the theoretical insight that amplifying a precise initial population gradient suffices for learning parities, we construct one possible progress measure, which is a function of the sequence of weights $w_0, \ldots w_t \in \RR^n$ so far: $\rho(w_{0:t}) := \left\lVert w_t - w_0 \right\rVert_\infty$, using the fact that $w_t - w_0$ is an estimate for the initial population gradient in the linearized setting. Figure~\ref{fig:hidden-progress} shows how gradual weight movement (and thus, progress) can be hidden behind plateauing losses.

\paragraph{Roles of overparameterization vs. \textit{oversampling}.}
An interesting consequence of our analysis is that it illuminates scaling behaviors with respect to a third fundamental resource parameter: \emph{model size}, which we study in terms of network width $r$.
If SGD operated by a ``random search'' mechanism, one would expect width to provide a parallel speedup. Instead, we SGD sequentially amplifies progress. The sharp lower tails in Figure~\ref{fig:scaling-laws} \emph{(left)} imply that running $r$ identical copies of SGD does \emph{not} give $(1/r)\times$ speedups; more directly, in Figure~\ref{fig:extra-experiments} \emph{(left)}, convergence times for sparse parities empirically plateau at large model sizes.\looseness=-1

It is a significant challenge to understand the interactions between network \emph{depth} and computation, and largely outside the scope of this work. However, in Appendix~\ref{subsec:deep-quad}, we provide a brief note on using parities and polynomial-activation MLPs to construct a simple counterexample to the \emph{``deep only works if shallow is good''} principle of \cite{malach2019deeper}, demonstrating a case where a deep network can get near-perfect accuracy even when greedy layerwise training (e.g. \citep{belilovsky2019greedy}) cannot beat trivial performance.

\paragraph{Learning and grokking in the finite-sample multi-pass setting.}
The main theoretical and empirical results in this work consider online learning algorithms which couple the resources of training time and independent samples. However, due to the computational-statistical gap in parity learning, these positive results are suboptimal in terms of sample efficiency.
We find that minibatch SGD (with weight decay) can empirically solve sparse parities, even from a sample of size $m \ll n^k$. For small values of $m$, we reliably observe the \emph{grokking} phenomenon \citep{power2022grokking}: overfitting for a long time, then a \emph{delayed} phase transition for the generalization error; see the two center panels of Figure~\ref{fig:extra-experiments} \emph{(right)}.\looseness=-1
\section{Conclusion}

This work puts forward parity learning as a stylized test case to explore some of the puzzling features of the role of computational (as opposed to statistical) resources in deep learning. These include discontinuous improvements (a.k.a. \emph{emergent capabilities}), feature learning, and universality of architectures. In particular, we show that deep learning on parities exhibits a phase transition behavior, that it is successfully learned by a variety of deep-net architectures, and that this success cannot be explained as a ``random exhaustive search'', nor through frameworks such as the neural tangent kernel or layer-by-layer learning.

However, there are more experimental and theoretical questions, even for this simplified case of parity learning. Our focus in this work was on the online learning case, where training time and samples arise in tandem. However, we believe it would be instructive to investigate parity learning when three resources of samples, time, and model size are scaled separately. Some very preliminary findings along these lines are presented in Section~\ref{sec:empirics}.

Extending our theoretical results to the small-batch setting, as well as to more architectures, is an open problem. Resolving it would require a better understanding of anti-concentration (lower bound on deviation from mean) of Fourier coefficients, a phenomenon that is much less studied than the concentration of these coefficients. We would also want to extend the analysis beyond parities to tasks that are not aligned with the elementary basis such as low-rank tensor recovery.

Another important follow-up direction is understanding the extent that these insights extend from parity learning to real-world problems, as well as the extent into which non-synthetic tasks (in, e.g., natural language processing and program synthesis) embed within them parity-like subtasks of exhaustive combinatorial search.

\paragraph{Broader impact.} This work seeks to contribute to the foundational understanding of computational scaling behaviors in deep learning. Our theoretical and empirical analyses are in a heavily-idealized synthetic problem setting. Hence, we see no direct societal impacts of the results in this study.

\paragraph{Acknowledgements.} We would like to thank Lenka Zdeborová for providing us with references to the statistical physics literature on phase transitions in the learning curves of neural networks.
\section{Proofs}

\subsection{Global convergence for MLPs}
For some even number $r$, consider a ReLU MLP of size $r$:
\[
f(x;\theta) = \sum_{i=1}^r u_i\sigma(w_i^\top x+b_i)
\]
Where $\sigma$ is the ReLU activation $\sigma(x) = \max \{x,0\}$, $w_1, \dots, w_r \in \RR^{r \times n}, b \in \RR^r, u \in \RR^r$ and we denote the set of parameters by $\theta = \{w_1, \dots, w_r, b, u\}$. We denote by $u_i^{(t)}, w_i^{(t)}, b^{(t)}$ and $\theta_t$ the value of the relevant parameters at iteration $t$ of gradient-descent. For brevity, we sometimes denote $w_i = w_i^{(0)}, b_i = b_i^{(0)}, u_i = u_i^{(0)}$.
W.l.o.g., we assume that $S = [k]$, and so $y = \chi_{[k]}(x) = \prod_{i=1}^k x_i$. Indeed, since the weights' initialization we consider is permutation symmetric, this does not limit the generality of the results. We take $\ell$ to be the \textbf{hinge loss}: $\ell(y, \hat{y}) = \max \{1-y\hat{y}, 0\}$. We use the following unbiased initialization:
\begin{itemize}
    \item For all $1 \le i \le r/2$, randomly initialize
    \[
    w^{(0)}_i \sim \mathrm{Unif} \left(\{\pm 1\}^n\right), u_i^{(0)} \sim \mathrm{Unif} \left(\{\pm 1\}\right), b_i^{(0)} \sim \mathrm{Unif}(\left\lbrace-1+1/k, \dots, 1-1/k\right\rbrace)
    \]
    \item For all $r/2 < i \le r$, initialize
    \[
    w^{(0)}_i = w^{(0)}_{i-r/2}, b_i^{(0)} = b_{i-r/2}^{(0)}, u^{(0)}_i = -u^{(0)}_{i-r/2}
    \]
\end{itemize}

We start by computing the population gradient at initialization. Using the fact that $\ell'(0,y) = -y$ we get the following:
\begin{equation}
\mathop{\mathbb{E}}\left[\nabla_{w_{i,j}} \ell(f(x;\theta_0), y)\right] =
\mathop{\mathbb{E}} \left[ -y\nabla_{w_{i,j}} f(x;\theta_0) \right] \label{eq:pop-grad}
\end{equation}
\[= \mathop{\mathbb{E}} \left[ -y u_i 1\left\lbrace{w_{i} \cdot x+b_i > 0}\right\rbrace x_j\right].\]

For $j \in [k]$ we have:
\[
\mathop{\mathbb{E}}\left[\nabla_{w_{i,j}} \ell(f(x;\theta_0), y)\right] = -u_i\mathop{\mathbb{E}} \left[\left(\prod_{j' \in [k] \setminus \{j\}} x_{j'} \right) 1\left\lbrace{w_{i} \cdot x+b_i > 0}\right\rbrace\right]
\]
For $j \notin [k]$ we have:
\[
\mathop{\mathbb{E}}\left[\nabla_{w_{i,j}} \ell(f(x;\theta_0), y)\right] = -u_i\mathop{\mathbb{E}} \left[\left(\prod_{j' \in [k] \cup \{j\}} x_{j'}\right) 1\left\lbrace{w_{i} \cdot x+b_i > 0}\right\rbrace\right]
\]
Finally, we have:
\[
\mathop{\mathbb{E}}\left[\nabla_{b_i} \ell(f(x;\theta_0), y)\right] = -u_i\mathop{\mathbb{E}} \left[\left(\prod_{j' \in [k]} x_{j'}\right) 1\left\lbrace{w_{i} \cdot x+b_i > 0}\right\rbrace\right]
\]

Denote $$g_{i,j} = \mathbb{E} \left[\nabla_{w_{i,j}} \ell(f(x;\theta_0), y)\right],~ \gamma_i = \mathbb{E} \left[\nabla_{b_i} \ell(f(x;\theta_0), y)\right]$$

For some function $f$ and some subset $S \subseteq [n]$, denote $\hat{f}(S) = \mathop{\mathbb{E}}[f(x)\chi_S(x)]$.

Denote $\mathrm{LTF}_{w,b}(x) = 1\{w \cdot x+b > 0\}$ and
let $\mathrm{MAJ}(x) = \mathop{\mathrm{sign}} \left(\sum_{i=1}^n x_i\right)$ and observe that, if $|b| < 1$,
$$\mathrm{LTF}_{w,b}(x) = \frac{1}{2}+\frac{1}{2}\mathrm{MAJ}(w \odot x)$$
where $w \odot x = (w_1 x_1, \dots, w_n x_n) \in \{\pm 1\}$.
Since $\{\chi_S\}_{S \subseteq [n]}$ is a Fourier Basis, we can write $\mathrm{MAJ} = \sum_{S \subseteq [n]}\hat{\mathrm{MAJ}}(S) \chi_S$ and therefore:
\begin{align*}
\mathrm{LTF}_{w,b}(x) &= \frac{1}{2}+\frac{1}{2}\mathrm{MAJ}(w \odot x)
= \frac{1}{2}+\frac{1}{2}\sum_{S \subseteq [n]} \hat{\mathrm{MAJ}}(S) \chi_S(w \odot x)\\
&= \frac{1}{2}+\frac{1}{2}\sum_{S \subseteq [n]} \hat{\mathrm{MAJ}}(S) \chi_S(w) \chi_S(x)
\end{align*}
So, for every $S \subseteq [n]$ with $|S| \ge 1$ we have $\hat{\mathrm{LTF}}_{w,b}(S) = \frac{1}{2}\hat{\mathrm{MAJ}}(S) \chi_S(w)$ and so $|\hat{\mathrm{LTF}}_{w,b}(S)| = \frac{1}{2}|\hat{\mathrm{MAJ}}(S)|$.

\begin{lemma}
\label{lem:maj_bound}
Fix some $k$, and assume that $n \ge 2k^2$. Then, for every $S \subseteq [n]$ s.t. $|S| = k$ it holds that:
\begin{itemize}
    \item If $k$ is even: $\hat{\mathrm{MAJ}}(S) = 0$.
    \item If $k$ is odd:
    \[
    c_1^2 k^{-2/3} \le \binom{n}{k} \hat{\mathrm{MAJ}}(S)^2 \le c^2_2 k^{-2/3}
    \]
    for some universal constants $c_1, c_2$.
\end{itemize}
\end{lemma}

\begin{proof}
See \cite{o2014analysis}.
\end{proof}

Observe that for all $S \subseteq [n]$ s.t. $|S| = k$ it holds that $\hat{\mathrm{MAJ}}(S) = \hat{\mathrm{MAJ}}([k])$ and denote $\xi_k := \hat{\mathrm{MAJ}}([k])$.

Therefore, by the previous lemma we get that
for every even $k$ the following holds:
\begin{equation}
\label{eq:fourier_ratio}
\frac{\xi_{k+1}^2}{\xi_{k-1}^2} \le \frac{c^2_2 (k-1)^{2/3}\binom{n}{k-1}}{c^2_1(k+1)^{2/3}\binom{n}{k+1}} \le \frac{c_2}{c_1} \frac{(k+1)k}{(n-k)(n-k+1)} \le \frac{8c_2}{c_1} \frac{k^2}{n^2}
\end{equation}

\begin{lemma}
\label{lem:grad}
Assume that $k$ is even and that $n \ge 2(k+1)^2$.
Then, the following hold:
\begin{enumerate}
    \item If $j \in \{1, \dots, k\}$ then:
    \[
    g_{i,j} = -\frac{1}{2} u_i \xi_{k-1} \cdot \chi_{[k]\setminus \{j\}}(w_i)
    \]
    \item If $j \in \{k+1, \dots, n\}$ then:
    \[
    g_{i,j} = -\frac{1}{2} u_i \xi_{k+1} \cdot \chi_{[k] \cup \{j\}}(w_i)
    \]
    \item $\gamma_i = 0$.
\end{enumerate}
\end{lemma}

\begin{proof}
If $j \in [k]$ then:
\begin{align*}
g_{i,j} &= -u_i \mathop{\mathbb{E}} \left[\chi_{[k] \setminus \{j\}}(x) \mathrm{LTF}_{w_i,b_i}(x)\right] = -u_i \hat{\mathrm{LTF}}_{w_i,b_i}([k] \setminus \{j\}) \\
&= - \frac{1}{2} u_i \hat{\mathrm{MAJ}}([k]\setminus \{j\})\chi_{[k]\setminus \{j\}}(w_i) = -\frac{1}{2} u_i \xi_{k-1} \cdot \chi_{[k]\setminus \{j\}}(w_i)
\end{align*}
Similarly, if $j \notin [k]$ we have:
\begin{align*}
g_{i,j} &= -u_i \mathop{\mathbb{E}} \left[\chi_{[k] \cup \{j\}}(x) \mathrm{LTF}_{w_i,b_i}(x)\right] = -u_i \hat{\mathrm{LTF}}_{w_i,b_i}([k] \cup \{j\}) \\
&= - \frac{1}{2} u_i \hat{\mathrm{MAJ}}([k]\cup \{j\})\chi_{[k]\cup \{j\}}(w_i) = -\frac{1}{2} u_i \xi_{k+1} \cdot \chi_{[k]\cup \{j\}}(w_i)
\end{align*}
Finally, we have:
\[
\gamma_i = -u_i \hat{\mathrm{LTF}}_{w_i,b_i}([k]) = -u_i \hat{\mathrm{MAJ}}([k]) = 0
\]
\end{proof}

\begin{lemma}
\label{lem:mlp_partial_step_one}
Let $\tau > 0$ be some tolerance parameter, fix $\epsilon \in (0,1)$ and let $\eta = \frac{1}{k|\xi_{k-1}|}$. Assume that $k$ is an even number.
Fix some $w_1, \dots, w_k \in \{\pm 1\}^n$, $b_1, \dots, b_r \in (-1,1)$ and $u_1, \dots, u_k \in \{\pm 1\}$.
Let $\hat{w}_i = -\eta \hat{g}_i$ and $\hat{b}_i = b_i-\eta \hat{\gamma}_{i}$ s.t. $\|\hat{g}_i - g_i\|_\infty \le \tau$ and $\|\hat{\gamma}_{i}-\gamma_i\| \le \tau$. Assume the following holds:
\begin{itemize}
    \item For all $i,j \in [k]$ it holds that $w_{i,j} = u_i \cdot \mathop{\mathrm{sign}} \xi_{k-1}$.
    \item $b_i = -\frac{1}{2} + \frac{i+1}{k}$
\end{itemize}
Then, if $\tau \le \frac{|\xi_{k-1}|}{16k\sqrt{2n \log(2k/\epsilon)}}$ and $n \ge \frac{2^{11} k^4c_2^2\log (2k/\epsilon)}{c_1^2}$ there exists some $\hat{u} \in \mathbb{R}^k$ with $\|\hat{u}\|_\infty \le 8k$ s.t.
$f(x) = \sum_{i=1}^k \hat{u}_i \sigma(\hat{w}_i \cdot x+\hat{b}_i)$
satisfies $$\mathop{\mathbb{E}}_x[\ell(f(x),\chi_{[k]}(x))] \le 16 \epsilon k^2 n$$
Additionally, for all $i$ and all $x$ it holds that $|\sigma(\hat{w}_i \cdot x + \hat{b}_i)| \le n+1$.
\end{lemma}

\begin{proof}
We start with the following claim.

\textbf{Claim 1}: For all $i$ and for all $j \in [k]$ it holds that $\left|\hat{w}_{i,j} -\frac{1}{2k}\right| \le \frac{\tau}{|\xi_{k-1}|}$.

\textbf{Proof}: First, observe that by the assumption it holds that for all $i,j \in [k]$,
\[u_i \cdot \chi_{[k] \setminus \{j\}}(w_i) = u_i \cdot \prod_{j' \in [k] \setminus \{j\}}^k w_{i,j} = u_i \cdot (u_i \cdot \mathop{\mathrm{sign}} \xi_{k-1})^{k-1} = \mathop{\mathrm{sign}} \xi_{k-1}
\]
Now, from Lemma \ref{lem:grad} we have:
\[
g_{i,j} = -\frac{1}{2} u_i \xi_{k-1} \cdot \chi_{[k]\setminus \{j\}}(w_i) = - \frac{1}{2}|\xi_{k-1}|
\]
and so
\[
\left|w_{i,j}-\frac{1}{2k}\right| = \left|-\eta\hat{g}_{i,j}-\frac{1}{2k}\right| =
\frac{1}{k}\left|-k\eta \hat{g}_{i,j} +\frac{g_{i,j}}{|\xi_{k-1}|} \right| = \frac{1}{k|\xi_{k-1}|} \left|g_{i,j}-\hat{g}_{i,j}\right| \le \frac{\tau}{|\xi_{k-1}|}
\]

\textbf{Claim 2}: For all $i$
 and for all $j > k$ it holds that $|\hat{w}_{i,j}| \le \frac{|\xi_{k+1}|+2\tau}{2k|\xi_{k-1}|}$

\textbf{Proof}: Using Lemma \ref{lem:grad} we have:
\[
|\hat{w}_{i,j}| = \eta |\hat{g}_{i,j}| \le \eta(|g_{i,j}| + |g_{i,j}-\hat{g}_{i,j}|) \le \eta \left( \frac{|\xi_{k+1}|}{2} +\tau\right) = \frac{|\xi_{k+1}|+2\tau}{2k|\xi_{k-1}|}
\]

\textbf{Claim 3}: For all $i$
 it holds that $|\hat{b}_{i}-b_i| \le \frac{\tau}{k|\xi_{k-1}|}$

\textbf{Proof}: Using Lemma \ref{lem:grad} we have:
\[
|\hat{b}_{i}-b_i| = \eta |\hat{g}_{i,0}| = \eta|g_{i,j}-\hat{g}_{i,j}| \le \eta \tau = \frac{\tau}{k|\xi_{k-1}|}
\]

\textbf{Claim 4}: Fix $\delta > 0$. Let $h_i$ be a function s.t. $h_i(x) = \sigma(\frac{1}{2k}\sum_{j=1}^k x_j+b_i)$ and $\hat{h}_i$ a function s.t. $\hat{h}_i(x) = \sigma(\hat{w}_{i,j} \cdot x+\hat{b}_i)$. Then, if $\tau \le \frac{\delta}{2} \frac{k|\xi_{k-1}|}{\sqrt{2n \log(2k/\epsilon)}}$ and $n \ge \frac{32c_2^2\log (2k/\epsilon)}{c_1^2 \delta^2}$, the following holds:
\begin{enumerate}
    \item $\mathbb{P}_{x \sim \{\pm 1\}^n} \left[|h_i(x) - \hat{h}_i(x)| \ge \delta \right] \le \frac{\epsilon}{k}$
    \item $|\hat{h}_i(x)| \le n+1$
\end{enumerate}

\textbf{Proof}: Let $x \sim \{\pm 1\}^n$, and denote $\tilde{x}_j = \hat{w}_{i,j} x_j$. Denote $\Delta = \frac{|\xi_{k+1}|+2\tau}{2k|\xi_{k-1}|}$. So, for $j > k$, $\tilde{x}_j$ is a random variable satisfying $|\tilde{x}_j| \le \Delta$. Furthermore, it holds that $\mathop{\mathbb{E}}\left[\sum_{j > k} \tilde{x}_j\right] = 0$. Therefore, from Hoeffding's inequality:
\[
\mathbb{P}\left[\left|\sum_{j > k} \tilde{x}_j\right| \ge \Delta \sqrt{\frac{n \log (2k/\epsilon)}{2}}\right] \le 2 \exp \left(-\frac{2\left(\Delta \sqrt{\frac{n \log (2k/\epsilon)}{2}}\right)^2}{n\Delta^2}\right) \le \frac{\epsilon}{k}
\]
Now, fix some $x$ s.t. $\left|\sum_{j>k}\tilde{x}_j\right| \le \Delta \sqrt{\frac{n \log (2k/\epsilon)}{2}}$. In this case we have:
\begin{align*}
|h_i(x)-\hat{h}_i(x)|
&\le \left|\frac{1}{2k}\sum_{j=1}^k x_j + b_i - \hat{w}_i \cdot x + \hat{b}_i\right| \\
&\le \sum_{j=1}^k\left| \frac{1}{2k}x_j-\hat{w}_{i,j}x_j\right| + \left|\sum_{j > k} \hat{w}_{i,j} x_j\right| + \left|b_i-\hat{b}_i\right| \\
&\le \frac{k \tau}{|\xi_{k-1}|} + \Delta \sqrt{\frac{n \log (2k/\epsilon)}{2}} + \frac{\tau}{k|\xi_{k-1}|} \\
&= \frac{\tau(k^2+1)}{k|\xi_{k-1}|} + \frac{|\xi_{k+1}|+2\tau}{2k|\xi_{k-1}|} \cdot \sqrt{\frac{n \log(2k/\epsilon)}{2}} \\
&= \frac{\tau\left(\sqrt{2}(k^2+1)+\sqrt{n \log (2k/\epsilon)}\right)}{\sqrt{2}k|\xi_{k-1}|} + \frac{|\xi_{k+1}|\sqrt{n \log (2k/\epsilon)}}{2\sqrt{2}k|\xi_{k-1}|} \\
&\le \frac{\tau \sqrt{2n \log (2k/\epsilon)}}{k|\xi_{k-1}|} + \frac{2\sqrt{2}c_2\sqrt{\log (2k/\epsilon)}}{c_1\sqrt{n}}
\end{align*}
where in the last inequality we use Eq. (\ref{eq:fourier_ratio}).
So, choosing $\tau \le \frac{\delta}{2} \frac{k|\xi_{k-1}|}{\sqrt{2n \log(2k/\epsilon)}}$ and $n \ge \frac{32c_2^2\log (2k/\epsilon)}{c_1^2 \delta^2}$ gives the required.

\textbf{Claim 5:} Let $h_1, \dots, h_k$ be the functions defined in the previous claim. Then, there exists weights $u^*$ with $\|u^*\|_\infty \le 8k$ s.t. for
$f^*(x) = \sum_{i=1}^k u^*_i h_i(x)$
it holds that $f^*(x) = 2\chi_{[k]}(x)$ for all $x \in \{\pm 1\}^n$.

\textbf{Proof:} For $i \le k-2$ define $u^*_i = 8k(-1)^{i+1}$ and $u^*_{k-1} = 6k$, $u^*_{k} = -2k$.

\textbf{Proof of Lemma \ref{lem:mlp_partial_step_one}:}
Choose $\hat{u} = u^*$. Using Claim 4 and the union bound, w.p. $1-\epsilon$ over $x \sim \{\pm1\}^n$ it holds that for all $i \in [k]$, $|h_i(x) -\hat{h}_i(x)| \le \delta$. Therefore, w.p. $\ge 1- \epsilon$
\[
|f(x)-f^*(x)| = \left|\sum_{i=1}^k u_i^*(h_i(x)-\hat{h}_i(x)) \right| \le \sum_{i=1}^{k} \left|u_i^*\right| \left|h_i(x)-\hat{h}_i(x) \right| \le 8k^2\delta
\]
so, choosing $\delta = \frac{1}{8k^2}$ we get that, w.p. at least $1-\epsilon$ over the choice of $x$ it holds that $f(x)\chi_{[k]}(x) \ge 1$.
Additionally, for every $x$ it holds that
\[
|f(x)|  \le \sum_{i=1}^{k} \left|u_i^*\right| \left|\hat{h}_i(x) \right| \le 8k^2 (n+1) \le 16k^2 n
\]
Therefore, we get:
\begin{align*}
\mathop{\mathbb{E}}_x\left[\ell(f(x), \chi_{[k]}(x))\right]
&\le \epsilon \mathop{\mathbb{E}}_x \left[\ell(f(x), \chi_{[k]}(x)) | f(x) \chi_{[k]}(x) < 1\right] \\
&\le \epsilon \mathop{\mathbb{E}}_x \left[|f(x)| | f(x) \chi_{[k]}(x) < 1\right] \le 16 \epsilon k^2 n
\end{align*}

 \end{proof}

\begin{lemma}
\label{lem:empirical_one_step}
Assume we randomly initialize an MLP using the unbiased initialization defined previously. Consider the following update:
\[
w_i^{(1)} = (1-\lambda_0)w_i^{(1)} - \eta_0 \hat{g}_{i}, ~b_i^{(1)} = b_i^{(1)} - \eta_0 \hat{\gamma}_{i}
\]
where
\[
\hat{g}_i = \frac{1}{B} \sum_{l=1}^B\nabla_{w_i} \ell(f(x_{0,l};\theta_0),y_{0,l}),
~\hat{\gamma}_i = \frac{1}{B} \sum_{l=1}^B \nabla_{b_i} \ell(f(x_{0,l};\theta_0), y_{0,l})
\]

Let $k$ be even number. Then, for every $\epsilon,\delta \in (0,1/2)$, denoting $\tau = \frac{|\xi_{k-1}|}{16k\sqrt{2n \log(2k/\epsilon)}}$,  if $\eta = \frac{1}{k |\xi_{k-1}|}$, $\lambda_0 = 1$, $r \ge k \cdot 2^{k} \log(k/\delta)$, $B \ge \frac{2}{\tau^2} \log(4nr/\delta)$  and $n \ge \frac{2^{11} k^4c_2^2\log (2k/\epsilon)}{c_1^2}$, w.p. at least $1-2\delta$ over the initialization and the sample, there exists $\hat{u} \in \mathbb{R}^r$ with $\|\hat{u}\|_\infty \le 8k$ and $\|\hat{u}\|_2 \le 8k \sqrt{k}$ s.t. $f(x) = \sum_{i=1}^r \hat{u}_i \sigma\left(w_i^{(1)} \cdot x+b^{(1)}_i\right)$ satisfies
\[
\mathop{\mathbb{E}}_x[\ell(f(x), \chi_{[n]}(x))] \le 16 \epsilon k^2 n
\]
Additionally, it holds that $\|\sigma(W^{(1)} \cdot x + b^{(1)})\|_\infty \le n+1$.
\end{lemma}

\begin{proof}
Note that by the choice of initialization it holds that $f(x;W^{(0)}) = 0$, and by the assumption on the loss function $\ell'(f(x;W^{(0)}), y) = -y$. Therefore, we get that $\mathop{\mathbb{E}}\left[\nabla_{w_{i}} \ell(f(x;W^{(0)}), y)\right] = g_i$ and $\mathop{\mathbb{E}}\left[\nabla_{b_{i}} \ell(f(x;W^{(0)}), y)\right] = \gamma_i$.

\textbf{Claim:} with probability at least $1-\delta$,
\begin{equation}
\label{eq:emp_grad_bound}
    \mathrm{for~all~} i,j:~\left\lVert\hat{g}_i-\mathop{\mathbb{E}}\left[\nabla_{w_{i}} \ell(f(x;W^{(0)}), y)\right]\right\rVert_\infty \le \tau~\mathrm{and}~ \left\lvert\hat{\gamma}_{i}-\mathop{\mathbb{E}}\left[\nabla_{b_{i}} \ell(f(x;W^{(0)}), y)\right]\right\rvert \le \tau
\end{equation}

\textbf{Proof:} Fix some $i,j$ and note that by Hoeffding's inequality,
\[
\Pr \left[|\hat{g}_{i,j} - \mathop{\mathbb{E}} \hat{g}_{i,j}| \ge \tau\right] \le 2 \exp \left(-B\tau^2/2\right) \le \frac{\delta}{nr + r}
\]
and similarly we get $\Pr \left[|\hat{\gamma}_i - \mathop{\mathbb{E}} \hat{\gamma}_i| \ge \tau \right] \le \frac{\delta}{nr+r}$. The required follows from the union bound.

Now, assume that Eq. (\ref{eq:emp_grad_bound}) holds.
For some random $w_i \sim \{\pm 1\}^n$, the probability that $w_{i,j} = w_{i,j'}$ for all $j,j' \in [k]$ is $2^{-k+1}$. Additionally, for some fixed $i'\in [k]$, the probability that $b_i = -\frac{1}{2} + \frac{i'}{k}$ is $\frac{1}{2k}$. Therefore, for some fixed $i \in [r/2]$ and $i' \in [k]$, with probability $\frac{1}{k\cdot 2^{k-1}}$, $b_i = b_{i+r/2} = -\frac{1}{2} + \frac{i'}{k}$ and either $w_{i,j} = u_i \mathop{\mathrm{sign}} \xi_{k-1}$ or $w_{i+r/2,j} = u_{i+r/2,j} \mathop{\mathrm{sign}} \xi_{k-1}$. Taking $r \ge k \cdot 2^{k} \log(k/\delta)$, we get that the probability that there is no $i \in [r/2]$ that satisfies the above condition (for fixed $i'$) is:
\[
\left(1-\frac{1}{k 2^{k-1}}\right)^{r/2} \le \exp \left(-\frac{r}{2k2^{k-1}}\right) \le \frac{\delta}{k}
\]
Using the union bound, with probability at least $1-\delta$, there exists a set of $k$ neurons satisfying the conditions of Lemma \ref{lem:mlp_partial_step_one}, and therefore the required follows from the Lemma.
\end{proof}

We use the following well-known result on convergence of SGD (see for example \cite{shalev2014understanding}):

\begin{theorem}
\label{thm:sgd}
Let $M, \rho > 0$. Fix $T$ and let $\eta = \frac{M}{\rho\sqrt{T}}$ Let $F$ be a convex function and $u^* \in \arg \min_{\|u\|_2 \le M} f(u)$. Let $u^{(0)} = 0$ and for every $t$, let $v_t$ be some random variable s.t. $\mathop{\mathbb{E}} \left[v_t | u^{(t)}\right] = \nabla_{u^{(t)}} F(u^{(t)})$ and let $u^{(t+1)} = u^{(t)} - \eta v^{(t)}$. Assume that $\|v_t\|_2 \le \rho$ w.p. $1$. Then,
\[
\frac{1}{T} \sum_{t=1}^T F(u^{(t)}) \le F(u^*) + \frac{M \rho}{\sqrt{T}}
\]
\end{theorem}

We prove the following theorem:

\begin{theorem}
Let $k$ be an even number.
Assume we randomly initialize an MLP using the unbiased initialization defined previously. Fix $\epsilon \in (0,1/2)$ and let
let $T \ge \frac{2^9k^3rn^2}{\epsilon^2}, r \ge k \cdot 2^{k} \log(8k/\epsilon), B \ge c_1^{-1}2^8k^{7/6} n \binom{n}{k-1} \log(128k^3n/\epsilon)  \log(32nr/\epsilon), n \ge \frac{2^{11} k^4c_2^2\log (128 k^3 n/\epsilon)}{c_1^2}$. Choose the following learning rate and weight decay schedule:
\begin{itemize}
    \item At the first step, use $\eta_0 = \frac{1}{k |\xi_{k-1}|}$, $\lambda_0 = 1$ for all weights.
    \item After the first step, use $\eta_t = 0$ for the first layers weights and biases and $\eta_t = \frac{4k^{1.5}}{n\sqrt{r(T-1)}}$ for the second layer weights, with $\lambda_t = 0$ for both layers.
    \item Bias terms are never regularized.
\end{itemize}

Then, the following holds, with expectation over the randomness of the initialization and the sampling of the batches:
\[
\mathop{\mathbb{E}} \left[\min_{t \in [T]} \ell(f(x;\theta_t),y) \right] \le \epsilon
\]
\end{theorem}

\begin{proof}
Let $F(u) = \mathop{\mathbb{E}}_x\left[\ell(u^\top \sigma(W^{(1)}x + b^{(1)}),y)\right]$ and notice that $F$ is a convex function. For every $t$, denote
$$v_t = \frac{1}{B} \sum_{l=1}^B \nabla_{u^{(t)}} \ell(f(x_{t,l};\theta_t),y) = \frac{1}{B} \sum_{l=1}^B \nabla_{u^{(t)}}\ell\left(\left(u^{(t)}\right)^\top \sigma(W^{(1)}x_{l,t} + b^{(1)}),y_{l,t}\right)$$
where we use the fact that we don't update the weights of the first layer after the first step.
From the above we get $\mathop{\mathbb{E}}[v_t|u^{(t)}] = \nabla_{u^{(t)}}F(u^{(t)})$.

Now, we will show that w.h.p. there exists $u^*$ with good loss. Let $\epsilon' = \frac{\epsilon}{64 k^2 n}, \delta' = \frac{\epsilon}{8}$.
Denote $\tau =  \frac{|\xi_{k-1}|}{16k\sqrt{2n \log(128k^3n/\epsilon)}} = \frac{|\xi_{k-1}|}{16k\sqrt{2n \log(2k/\epsilon')}}$. Observe that $r \ge k \cdot 2^{k} \log(k/\delta')$, and using the fact that $|\xi_{k-1}| \ge c_1(k-1)^{-1/3}\binom{n}{k-1}^{-1}$ we get

$$B \ge \frac{2^8k^2\cdot  n \log(128k^3n/\epsilon)}{|\xi_{k-1}|} \log(32nr/\delta) = \frac{2}{\tau^2} \log(4nr/\delta')$$
and additionally
$n \ge \frac{2^{11} k^4c_2^2\log (2k/\epsilon')}{c_1^2}$.

From the above, applying Lemma \ref{lem:empirical_one_step} with $\epsilon', \delta'$ we get that w.p. $1-\epsilon/4$ there exists $u^* \in \mathbb{R}^r$ with $\|u^*\|_2 \le 8k \sqrt{k}$ s.t.
$F(u^*) \le \epsilon/4$
and for all $i$ and all $x$ it holds that $\|\sigma(W^{(1)} \cdot x + b^{(1)})\|_\infty \le n+1$. Using this, we get:
\[
\|v_t\|_2 \le \frac{1}{B}\sum_{l=1}^B \|\sigma(W^{(1)}x_{l,t} + b^{(1)}),y_{l,t} \|_2 \le \sqrt{r}(n+1)
\]
So, we can apply Theorem \ref{thm:sgd} with $M = 8k\sqrt{k}$ and $\rho = 2\sqrt{r}n$ and get that, w.p. $1-\epsilon/4$ over the initialization and the first step, it holds that
\begin{align*}
\mathop{\mathbb{E}}_{\mathrm{steps~ 2\dots T}} \left[\min_{t \in \{2, \dots, T\}} \ell(f(x;\theta_t),y)\right]
&\le \mathop{\mathbb{E}} \left[\frac{1}{T-1} \sum_{t=2}^T \ell(f(x;\theta_t),y)\right] \\
&= \mathop{\mathbb{E}}\left[\frac{1}{T-1} \sum_{t=2}^T F(u^{(t)})\right] \le \epsilon/4 + \frac{16k^{1.5}\sqrt{r}n}{\sqrt{T-1}} \le \epsilon/2
\end{align*}
Now, that after the first step we have $u^{(1)} = 0$ and therefore $\ell(f(x;\theta_1),y) = 1$, and so we always have $\min_{t \in [T]} \ell(f(x;\theta_t),y) \le 1$. Therefore, taking expecation over all steps we get:
\[
\mathop{\mathbb{E}}_{\mathrm{steps~1\dots T}} \left[\min_{t \in [T]} \ell(f(x;\theta_t),y)\right] \le \epsilon/2 + \epsilon/2 = \epsilon
\]

\end{proof}

\subsection{Proof of Proposition \ref{prop:mlp-from-gap}}
\begin{proof}
Denote $h(x) = \sigma'(w^\top x)$. Observe the gradient w.r.t. $w_i$:
\begin{align*}
    g^*_i = \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S} \left[\nabla_{w_i} \ell(y,f(x;w))\right]
    &= -\mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S} \left[\sigma'(w^\top x) yx_i\right] \\
    &= \begin{cases}
    \mathop{\mathbb{E}} \left[\sigma'(w^\top x) \prod_{j \in S \setminus \{i\}} x_i \right] & i \in S \\
    \mathop{\mathbb{E}} \left[\sigma'(w^\top x) \prod_{j \in S \cup \{i\}} x_i \right] & i \notin S \\
    \end{cases} \\
    &= \begin{cases}
    \hat{h}(S \setminus \{i\}) & i \in S \\
    \hat{h}(S \cup \{i\}) & i \notin S \\
    \end{cases}
\end{align*}
So, from the Fourier gap property, it holds that for all $i \in S$ we have $|g^*_i| \ge \gamma$ and for all $i \notin S$ we have $|g^*_i| \le \gamma/2$. From this, the required follows.
\end{proof}

\newcommand{\mathsf{dis}}{\mathsf{dis}}
\subsection{Global convergence for disjoint-PolyNets}\label{app:polynet}
In this section we will develop theory for disjoint-PolyNets trained with correlation loss, as illustrated in Figure~\ref{fig:disjoint-plot}. Section~\ref{app:flow} will consider optimization with gradient flow, and section~\ref{app:sgd} will consider optimization with batch size 1 SGD.

\subsubsection{Gradient flow analysis}\label{app:flow}

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/disjoint-polynet-plot.png}
    \vspace{-3mm}
    \caption{Example training curves for the disjoint-PolyNet trained with correlation loss, which is the setting of Appendix~\ref{app:polynet}. The left plot shows validation error under population gradient descent with small step size, approximating the setting of Section~\ref{app:flow}, and the right plot shows a run of SGD with batch size 1, as in Section~\ref{app:sgd}, though with a constant learning rate schedule. Initialization is i.i.d. standard Gaussian.}
    \label{fig:disjoint-plot}
\end{figure}

For $i$, $\mathop{\mathbb{E}}\left[\nabla_{w_i} \ell(f(x;w_{1:k}), y)\right] = 0$, so the irrelevant weights remain fixed at initialization. For each $i \in [k]$, let $v_i^{(t)}$ be the relevant weight in the $k$th partition.

In gradient flow, the relevant weights evolve according to the following differential equations:
\[
\forall i \in [n]: \quad \dot{v_i} = \prod_{j \neq i} v_j
\]

\begin{lemma}\label{lem:flow}
Suppose disjoint-PolyNet for $k>2$ is initialized such that $\prod_i v_i(0) > 0$, and optimized with gradient flow. Let $\bar{v}_a := \frac{1}{k}\sum_{i=1}^k (v_i(0))^2$, and $\bar{v}_g := \left(\prod_{i=1}^k (v_i(0))^2 \right)^{1/k}$. For any $b \geq 0$ and $i \in [k]$, let $T_i(b) := \arg\sup_{t \geq 0} (|v_i(t)| \leq b)$. Then
\[
T_i(b) \geq \frac{1}{k-2}\left(\bar{v}_a^{1-k/2} - (\bar{v}_a + b^2 - v_i(0)^2 )^{1-k/2} \right).
\]
Let $T_i(\infty) := \arg\sup_{t \geq 0} (|v_i(t)| < \infty)$. Then
\[
T_i(\infty) - T_i(b) \leq
\frac{1}{k-2} (\bar{v}_g + b^2 - v_i(0)^2)^{1-k/2}.
\]

\end{lemma}
\begin{proof}
First, observe that the product of the relevant weights is non-decreasing during gradient flow.
\[
\frac{d\left(\prod_{i=1}^k v_i(t)\right)}{d t}
= \sum_{i=1}^k \frac{\partial\left(\prod_{i=1}^k v_i\right)}{\partial v_i} \cdot \dot{v}_i
= \sum_{i=1}^k \left(\prod_{j \neq i} v_i\right)^2
> 0
\]

Thus, $\prod_i v_i(0) > 0$ implies that $\prod_i v_i(t) > 0$ for all $t$.

Observe that
\[
\frac{d v_i^2}{d t} = 2v_i \dot{v_i} = 2 \prod_{j = 1}^k v_j = 2\prod_{j = 1}^k |v_j|.
\]
This implies that for all $i, l \in [t]$,
\begin{equation}\label{eq:vi-equiv}
\frac{d v_i^2}{d t} = \frac{d v_l^2}{d t}.
\end{equation}
In other words, the squares of the relevant weights each follow the same trajectory, shifted according to their initializations. Let $q(t) := ({v_i(t)})^2 - ({v_i(0)})^2$, for any $i$. This quantity evolves as follows:
\[
\dot{q} = 2 \prod_{i=1}^k |v_i| = 2 \left(\prod_{i=1}^k (q(t) + (v_i(0))^2 )\right)^{1/2}
\]

Since $q(t)$ is strictly increasing, its inverse $q^{-1}$ is well-defined, and we can use the inverse function theorem to characterize $q^{-1}$ for all $t \geq 0$:
\[
q^{-1}(c) = \int_0^c \frac{1}{2} \left(\prod_{i=1}^k (\gamma + (v_i(0))^2 )\right)^{-1/2} d\gamma.
\]

We can upper- and lower-bound the integrand by applying Maclaurin's inequality (see page 52 in \citet{hardy1952inequalities}):
\[
\left(\gamma + \bar{v}_g \right)^k \leq \prod_{i=1}^k (\gamma + (v_i(0))^2 ) \leq \left(\gamma + \bar{v}_a \right)^k.
\]

The amount of time it takes for $q$ to reach a value of $c$ is thus:
\[
\int_0^c \frac{1}{2} \left(\prod_{i=1}^k (\gamma + (v_i(0))^2 )\right)^{-1/2} d\gamma \geq \int_0^c \frac{1}{2} (\gamma + \bar{v}_a)^{-k/2} d\gamma
= \frac{1}{k-2} \left(\bar{v}_a^{1-k/2} - (c+\bar{v}_a)^{1-k/2} \right).
\]

Hence, for any $b \geq 0$, for each $i$, $|v_i(t)| \leq b$ as long as
\[
t \leq \frac{1}{k-2}\left(\bar{v}_a^{1-k/2} - (\bar{v}_a + b^2 - v_i(0)^2 )^{1-k/2} \right).
\]

Meanwhile, the amount of time after $q^{-1}(c)$ it takes for $q$ to explode to infinity is
\[
\int_c^\infty \frac{1}{2} \left(\prod_{i=1}^k (\gamma + (v_i(0))^2 )\right)^{-1/2} d\gamma \leq \int_c^\infty \frac{1}{2} (\gamma + \bar{v}_g)^{-k/2} d\gamma
= \frac{1}{k-2} (c+\bar{v}_g)^{1-k/2}.
\]
Substituting $c = b^2 - v_i(0)^2$, we obtain that the amount of time it takes for $|v_i|$ to grow from $b$ to $\infty$ is
\[
\frac{1}{k-2} (\bar{v}_g + b^2 - v_i(0)^2)^{1-k/2}.
\]

We can upper- and lower-bound $\dot{q}$ by applying Maclaurin's inequality (see page 52 in \citet{hardy1952inequalities}):
\[
2\left(q(t) + \bar{v}_g \right)^{k/2} \leq \dot{q} \leq 2\left(q(t) + \bar{v}_a \right)^{k/2}.
\]

When $k=2$, solving the LHS and RHS differential inequalities yields:
\[
\bar{v}_g(e^{2t}-1) \ \leq\ q(t) \ \leq\ \bar{v}_a(e^{2t}-1).
\]

When $k>2$, we obtain:
\begin{equation}\label{eq:q-bounds-3}
(\bar{v}_g^{-(k/2-1)}-(k-2)t)^{-\frac{1}{k/2-1}} - \bar{v}_g \ \leq\ q(t) \ \leq\ (\bar{v}_a^{-(k/2-1)}-(k-2)t)^{-\frac{1}{k/2-1}} - \bar{v}_a
\end{equation}

From the lower bound on $q(t)$, we can infer that the relevant weights all explode to infinity by the following time:
\[
t = \frac{1}{(k-2)}\bar{v}_g^{-(k/2-1)}
= \frac{1}{(k-2)\left(\prod_i v_i(0)\right)^{1-2/k}}
\]

From the upper bound, we can infer that for any $c > 0$, it is the case that $q(t) \leq c$ so long as
\[
t \leq \frac{1}{k-2}\left(\bar{v}_a^{-(k/2-1)} - (\bar{v}_a + c)^{-(k/2-1)} \right).
\]
Hence, for each $i$, $|v_i(t)| \leq b$ for all
\[
t \leq \frac{1}{k-2}\left(\bar{v}_a^{-(k/2-1)} - (\bar{v}_a + b^2 - (v_i(0))^2 )^{-(k/2-1)} \right).
\]
\end{proof}

Now we analyze the relationship between the relevant weights and the accuracy of the disjoint-PolyNet.

For $y \in \RR$, let
\[
\mathop{\mathrm{sign}}(y) =
\begin{cases}
1 & \text{if}\ y>0 \\
0 & \text{if}\ y=0 \\
-1 & \text{if}\ y<0
\end{cases}
\]

Then define the error of $f$ with parameters $w_{1:k}$ as
\[
    \mathsf{err}(w_{1:k}) := \Pr_{x \sim \{\pm 1\}^n} [\mathop{\mathrm{sign}}(f(x;w_{1:k})) \neq \chi_S(x)]
\]

\begin{lemma}\label{lem:err}
Let $w$ be any setting of the weights of a disjoint-PolyNet such that $\prod_i v_i > 0$. For ease of notation, let $u_i := w_{i,2:n'}$ be the irrelevant portion of $w_i$. There is a constant $c$ such that
\[
\frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \left(\erf\left(\frac{|v_i|}{\|u_i\|_2\sqrt{2}}\right) + \frac{c\|u_i\|_\infty}{\|u_i\|_2^3}\right) \leq
\mathsf{err}(w_{1:k})
\leq
2\sum_{i=1}^k \exp\left(-\frac{|v_i|^2}{\|u_i\|_2^2}\right)
\]
where $\erf$ is the Gauss error function $\erf(y) := \frac{2}{\sqrt\pi}\int_0^y e^{-\tau^2}\,d\tau$.
\end{lemma}

\begin{proof}
Let $z_i = x_{(i-1)n'+1}$ be the $i$th relevant coordinate of $x$, and let $z_i^- = x_{(i-1)n'+2:in'}$ be the irrelevant coordinates in $P_i$.

Then we have that the error of $f$ with parameters $w_{1:k}$ is
\begin{align*}
    \mathsf{err}(w_{1:k}) &= \Pr_x \left[\mathop{\mathrm{sign}}\left(\prod_{i=1}^k w_i^\top x_{P_i}\right) \neq \chi_S(x)\right] \\
    &= \Pr_x \left[\mathop{\mathrm{sign}}\left(\prod_{i=1}^k w_i^\top x_{P_i}\right) = -\chi_S(x)\right] + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right]\\
    &= \Pr_x \left[\mathop{\mathrm{sign}}\left(\prod_{i=1}^k w_i^\top x_{P_i}\right) = -\chi_S(x)\right] + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right]\\
    &= \Pr_x \left[\#\{i \in [k]: \mathop{\mathrm{sign}}(v_i z_i + u_i^\top z_i^-) = -z_i \} \text{ is odd}\right] + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right] \\
    &= \Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- > v_i\} \text{ is odd}\right]/2 \\ &+ \Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- < -v_i\} \text{ is odd}\right]/2 + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right] \numberthis \label{line:indep}\\
    &= \Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- > v_i\} \text{ is odd}\right] + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right] \numberthis \label{line:sym}\\
    &= \Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- > |v_i|\} \text{ is odd}\right] + \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right] \numberthis \label{line:pos}
\end{align*}

Line~\ref{line:indep} follows because $u_i^\top z_i^-$ and $z_i$ are independent. In line~\ref{line:sym} we use that $u_i^\top z_i^-$ is symmetric about 0. Finally, line~\ref{line:pos} uses the assumption that $\prod_{i=1}^k v_i > 0$.

We can bound $\Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right]$ using Hoeffding's inequality:
\[
0 \leq \Pr_x \left[\prod_{i=1}^k w_i^\top x_{P_i} = 0\right]
\leq \sum_{i=1}^k \Pr_{z_i} [u_i^\top z_i^- = v_i]
\leq \sum_{i=1}^k \Pr_{z_i} [u_i^\top z_i^- \geq |v_i|]
\leq \sum_{i=1}^k \exp\left(-\frac{|v_i|^2}{\|u_i\|_2^2}\right).
\]

The indicator random variables $\ind[u_i^\top z_i^- > |v_i|]$ are independent of each other, so the first term in line~\ref{line:pos} can be characterized using the distribution of the parity of a sum of independent Bernoulli random variables. Let $X_i \sim \Ber(p_i)$ for $i \in [k]$, and let $X = \sum_{i=1}^k X_i$. The generating function for $X$ is $f(z) = \prod_{i=1}^k ((1-p_i) + p_i z)$. The parity of $X$ then satisfies

\[\Pr[X \text{ is odd}] = \frac{f(1)}{2} - \frac{f(-1)}{2}
= \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k (1 - 2p_i).
\]

First we will prove the upper bound on $\mathsf{err}$. First, observe that
\[
\frac{1}{2} - \frac{1}{2}\prod_{i=1}^k (1 - 2p_i) \leq \sum_{i=1}^k p_i.
\]
Thus,
\begin{align*}
\Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- > |v_i|\} \text{ is odd}\right] &\leq \sum_{i=1}^k \Pr_{z_i}[u_i^\top z_i^- > |v_i|] \\
&\leq \sum_{i=1}^k \exp\left(-\frac{|v_i|^2}{\|u_i\|_2^2}\right)
\end{align*}
by Hoeffding's inequality, and we are done.

Now we will prove the lower bound on $\mathsf{err}$. We have
\begin{align*}
\Pr_{z_1^-, \dots, z_k^-} \left[\#\{i \in [k]: u_i^\top z_i^- > |v_i|\} \text{ is odd}\right] &=
\frac{1}{2} - \frac{1}{2}\prod_{i=1}^k(1-2\Pr_{z_i}[u_i^\top z_i^- > |v_i|]) \\
&= \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \Pr_{z_i}[|u_i^\top z_i^-| \leq |v_i|]. \numberthis \label{eq:errform}
\end{align*}

We can bound this expression using the Berry-Esseen theorem \citep{berry1941accuracy, esseen1942liapunov}. Let $\beta \sim N(0,\|u_i\|_2^2)$. Then, the Berry-Esseen theorem states that there is a constant $c$ (which in practice can be $.56$) such that for any $i \in [k]$,
\[
\left| \Pr_{z_i}[|u_i^\top z_i^-| \leq |v_i|]
- \Pr[|\beta| \leq |v_i|] \right| \leq \frac{c\|u_i\|_\infty}{\|u_i\|_2^3}
\]
and we can use the characterization
\[
\Pr[|\beta| \leq |v_i|] = \erf\left(\frac{|v_i|}{\|u_i\|_2\sqrt{2}}\right).
\]

Plugging this into equation~\ref{eq:errform}, we obtain the lower bound on $\mathsf{err}$.
\end{proof}

First, we'll apply Lemma~\ref{lem:flow} and Lemma~\ref{lem:err} to the situation where the $w_i$'s have $\pm 1$ initialization. This generalizes Theorem~\ref{thm:fraction}.
\begin{corollary}\label{cor:pm1}

Suppose all the weights in a disjoint-Polynet are initialized randomly in $\pm 1$ and $k \geq 3$.

Let $T(\alpha) := \arg \sup_{t \geq 0}(\mathsf{err}(w_{1:k}(t)) \geq \alpha$. Then, for $\gamma \in (0, 1/2)$, if $\prod_i v_i(0) > 0$ (which happens w.p. 1/2),
\[
\frac{T(1/2 - \gamma)}{T(0)} = 1 - O((n')^{1-k/2} \cdot \gamma^{2/k-1}).
\]
Thus, even for $\gamma$ arbitrarily close to $0$, when the input is sufficiently long, the network spends almost all of training with error above $1/2 - \gamma$.

\end{corollary}
\begin{proof}

For $\pm 1$ initialization,
\[
\frac{\|u_i\|_\infty}{\|u_i\|_2^3} = (n')^{-3/2}
\]
so by Lemma~\ref{lem:err},
\begin{align*}
\mathsf{err}(w_{1:k}(t)) &\geq \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \left(\erf\left(\frac{|v_i(t)|}{\sqrt{2n'}}\right) + c(n')^{-3/2}\right) \\
&\geq \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \left(\sqrt{\frac{2}{\pi n'}}\cdot |v_i(t)| + c(n')^{-3/2}\right)
\end{align*}
for some $c > 0$.

By Lemma~\ref{lem:flow}, $|v_i(t)| \leq b$ for all $i$ whenever
\[
t \leq \frac{1}{k-2}\left(1 - b^{2-k} \right).
\]

Setting
\[
b = \sqrt{\pi n'/2}\left(\gamma^{1/k} - c(n')^{-3/2}\right),
\]
we obtain that $T(1/2 - \gamma) = \frac{1}{k-2} (1 - O((n')^{1-k/2} \cdot \gamma^{2/k-1}) )$.

Also by Lemma~\ref{lem:flow}, using the language of that lemma statement, for all $i$
\[
T_i(\infty) - T_i(b) \leq \frac{1}{k-2} \cdot b^{2-k} = \frac{1}{k-2}\cdot O((n')^{1-k/2}\cdot \gamma^{2/k-1}).
\]

Once all the relevant weights have exploded to infinity, the error of the network will have zero error, so the result follows.
\end{proof}

Now let us apply Lemma~\ref{lem:flow} and Lemma~\ref{lem:err} to the situation where the $w_i$'s have standard normal initialization. Again we find that with high probability, the phase of learning with near-trivial accuracy is much longer than the subsequent period until perfect accuracy, as illustrated by the left-hand plot in Figure~\ref{fig:disjoint-plot}.

\begin{corollary}
Suppose all the weights in a disjoint-Polynet are initialized $\sim N(0,1)$, and $k \geq 3$. Then, conditioned on $\prod_i v_i(0) > 0$, with probability $1-1/\mathrm{poly}(n')$ over the randomness of the initialization, for $\gamma \in (0, 1/2)$,
\[
\frac{T(1/2 - \gamma)}{T(0)} = 1 - \tilde{O}((n')^{1-k/2} \cdot \gamma^{2/k-1}).
\]
where $T(\cdot)$ is defined as in Corollary~\ref{cor:pm1}.

\end{corollary}
\begin{proof}
By a standard application of the generic Chernoff bound \citep{wainwright2019high}, for each $i \in [k], j \in [n'-1]$, we have
\[
\Pr[|u_{i,j}| > \tau] \leq 2e^{-\tau^2/2} \quad \text{for all}\ \tau \geq 0.
\]
Applying the union bound, we obtain
\[
\Pr[\exists i,j \text{ s.t. }|u_{i,j}| > \tau] \leq 2n'k e^{-\tau^2/2} \quad \text{for all}\ \tau \geq 0.
\]
This implies that w.p. $\geq 1-\epsilon/2$,
\begin{equation}\label{eq:uinfnorm}
\forall i: \quad \|u_i\|_\infty \leq \sqrt{2\log (4n'k/\epsilon)}.
\end{equation}

For each $i$, $\|u_i\|_2^2$ follows a chi-squared distribution with $n'-1$ degrees of freedom. By \citet{laurent2000adaptive}, for any $\tau \geq 0$,
\[
\|u_i\|_2^2 \in [(n'-1) - 2\sqrt{(n'-1)\tau}, (n'-1) + 2\sqrt{(n'-1)\tau} + 2\tau] \quad \text{w.p. } \geq 1-2e^{-\tau}.
\]
Hence, w.p. $\geq 1-\epsilon/2$,
\begin{equation}\label{eq:u2norm}
\forall i: \quad \|u_i\|_2 = \sqrt{n'} + O(\sqrt{\log(k/\epsilon)})
\end{equation}

With probability $1-\epsilon$ both $\|u_i\|_\infty$ and $\|u_i\|_2$ are bounded as above, in which case we obtain that for $\epsilon = 1/\mathrm{poly}(n',k)$, and for some $c_1, c_2 > 0$,
\[
\frac{\|u_i\|_\infty}{\|u_i\|_2^3} \le \frac{c_1\sqrt{\log (n'k)}}{\left(\sqrt{n'} + c_2\sqrt{\log(n'k)}\right)^{3}} \le \tilde{O}\left((n')^{-3/2}\right).
\]

Plugging this into Lemma~\ref{lem:err} gives us
\begin{align*}
\mathsf{err}(w_{1:k}(t)) &= \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \left(\erf\left(\frac{|v_i(t)|}{\|u_i\|_2\sqrt{2}}\right) + O\left(\frac{\|u_i\|_\infty}{\|u_i\|_2^3}\right)\right)\\
&\ge \frac{1}{2} - \frac{1}{2}\prod_{i=1}^k \left(\tilde{O}\left(\frac{|v_i(t)|}{\sqrt{n'}}\right) + \tilde{O}\left((n')^{-3/2}\right)\right).
\end{align*}

Thus, we can choose $ b = \tilde{\Omega}(\gamma^{1/k}\sqrt{n'})$ such that if $v_i(t) \leq b$ for all $i$, then $\mathsf{err}(w_{1:k}(t)) \geq \frac{1}{2}-\gamma$.

By Lemma~\ref{lem:flow}, for all $i$,
\[
T_i(b) \geq \frac{1}{k-2}\left(\bar{v}_a^{1-k/2} - (\bar{v}_a + b^2 - v_i(0)^2 )^{1-k/2} \right) \geq \frac{1}{k-2}\left(\bar{v}_a^{1-k/2} - (b^2 - O(\log(k)) )^{1-k/2} \right).
\]

By the Chernoff bound, $\max_i |v_i(0)| \leq \sqrt{2\log(4k/\delta)}$ w.p. $\geq 1-\delta/2$. And, by the same bound we used for $\|u_i\|_2^2$, we have that $\bar{v}_a = 1 + O(\sqrt{\log(1/\delta)/k})$ w.p. $1-\delta/2$. Thus, for $\delta = 1/\mathrm{poly}(k)$, we have that with probability $1-\delta$,
\begin{equation}\label{eq:tb}
T_i(b) \geq \frac{1}{k-2}\left(\tilde{\Omega}(1) - (b^2 - O(\log(k)) )^{1-k/2} \right)
= \frac{1}{k-2} \cdot \tilde{\Omega}(1).
\end{equation}

Also by Lemma~\ref{lem:flow},
\[
T_i(\infty) - T_i(b) \leq \frac{1}{k-2}(\bar{v}_g + b^2 - v_i(0)^2 )^{1-k/2}
\leq \frac{1}{k-2}(b^2 - v_i(0)^2 )^{1-k/2}
.
\]

With probability $1-1/\mathrm{poly(k)}$,
\begin{equation}\label{eq:tinf}
T_i(\infty) - T_i(b) \leq \frac{1}{k-2}(b^2 - O(\log(k)))^{1-k/2}
= \frac{1}{k-2} \cdot \tilde{O}\left((n')^{1-k/2} \cdot \gamma^{2/k - 1}\right).
\end{equation}

Combining Equation~\ref{eq:tb} and Equation~\ref{eq:tinf}, we obtain the desired statement.

\end{proof}

\subsubsection{SGD analysis}\label{app:sgd}
In this section we will analyze the training of a disjoint-PolyNet using SGD with online (i.i.d.) batches. We will show a convergence result for the 0-1 error of the learned classifier.
\begin{theorem}
Assume we randomly initialize the disjoint-PolyNet with weights drawn uniformly from $\{\pm 1\}$. Fix $\epsilon \in (0,1/2)$ and run SGD at any batch size $B \geq 1$ for $T\ge 6\log (2nT/\delta)\log(2k/\epsilon) (3n'-2)^{2k - 1}$ iterations. There exists an adaptive learning rate schedule, such that, with probability 1/2 over the randomness of the initialization and $1 - \delta$ over the sampling of SGD, the following holds:
\[
\mathrm{err}\left(w_{1:k}^{(T+1)}\right) \le \epsilon.
\]
\end{theorem}
\begin{proof}
For simplicity of presentation, we will assume $B = 1$.
Let the sample at iteration $t$ be $(x^{(t)}, y^{(t)})$ where $x^{(t)} \sim \mathsf{Unif}(\{\pm1\}^n)$ and $y^{(t)} = \chi_S^{(t)} (x^{(t)})$. Denote the population and stochastic gradient at time $t$ as:
\begin{align*}
    \hat{g}_i^{(t)} &= - y^{(t)}\left(\prod_{j\ne i}{w_j^{(t)}}^\top x_j^{(t)}\right)x_i^{(t)}\\
    g_i^{(t)} &= -\left(\prod_{j \ne i} w_{j,1}^{(t)} \right)e_1
\end{align*}
Observe that $|\hat{g}_{i,j}^{(t)}|, |g_{i,j}^{(t)}| \le \prod_{l \ne i} \|w_l^{(t)}\|_1$; note that this is also true when $B > 1$. Let the learning rate be $\eta^{(t)}_i = \frac{1}{2\sqrt{2T \log (2nT/\delta)}\cdot \prod_{l \ne i} \|w_{l}^{(t)}\|_1}$.
Let $\Delta_{i,j}^{(t)} = \eta^{(t)}_i(g_{i,j}^{(t)} - \hat{g}_{i,j}^{(t)})$ and $s_{i,j}^{(t)} = \sum_{\tau=1}^t \Delta^{(\tau)}_{i,j}$.
Observe that $\mathop{\mathbb{E}} \Delta^{(t)}_{i,j} = 0$ and therefore $s^{(1)}_{i,j}, \dots, s^{(t)}_{i,j}$ form a martingale. Furthermore, observe that
\begin{align*}
|s_{i,j}^{(t)}-s^{(t-1)}_{i,j}| = |\Delta^{(t)}_{i,j}| \le \frac{|g_{i,j}^{(t)}| + |\hat{g}_{i,j}^{(t)}|}{\sqrt{T\log (2nT/\delta)} \cdot \prod_{l \ne i} \|w_{l}^{(t)}\|_1} \le \frac{1}{\sqrt{2T \log (2nT/\delta)}}
\end{align*}
Then, for every $t < T$ and $i,j$, by the Azuma-Hoeffding inequality, with probability $1- \frac{\delta}{nT}$, $|s^{(t)}_{i,j}| \le 1/2$.
By the union bound, w.p. at least $1-\delta$, for every $t < T$ and all $i,j$ it holds that $|s^{(t)}_{i,j}| \le 1/2$. Let us assume this holds.

\begin{claim}\label{claim:upper}
For all $t \le T + 1$, for $j > 1$, $|w_{i,j}^{(t)}| \le 3/2$.
\end{claim}
\begin{proof}
Note that $ w_{i,j}^{(t+1)}= w_{i,j}^{(1)}  + s_{i,j}^{(t)}$, therefore, we have,
\begin{align*}
    | w_{i,j}^{(t+1)}|&\le  |w_{i,j}^{(1)}| + |s_{i,j}^{(t)}| \le  3/2.
\end{align*}
since $|w_{i,j}^{(1)}| = 1$ and $|s_{i,1}^{(t-1)}| \le 1/2$.
\end{proof}

Denote $\xi_i = \mathop{\mathrm{sign}} \left(w_{i,1}^{(1)}\right)$ and assume that $\prod_{j=1}^k \xi_j> 0$, and note that this happens w.p. $1/2$. We will assume this holds for the next lemma.
\begin{claim} \label{claim:lower}
For all $t \le T+1$ and $i \in [k]$ it holds that
\[
\xi_i w_{i,1}^{(t)} = |w_{i,1}^{(t)}| \ge \frac{1}{2} + \frac{1}{4\sqrt{2}}\left(\frac{1}{3n'-2}\right)^{k-1} \sqrt{\frac{t-1}{\log (nT/\delta)}}.
\]
\end{claim}
\begin{proof}
Observe that the claim holds for $t=1$ since $|w_{i,`}^{(1)}| = 1$. By induction on $t$, assume the claim holds for all $\tau \le t$. Now we will prove it holds for $t+1$.

Observe that for all $\tau \le t$, by the assumption:
\[
g_{i,1}^{(\tau)} = - \prod_{j \ne i}w_{j,1}^{(\tau)}= - \xi_i\prod_{j \ne i}\xi_j w_{j,1}^{(\tau)}= - \xi_i \prod_{j \ne i} |w_{j,1}^{(\tau)}|.
\]
Note that $ w_{i,1}^{(t+1)}= w_{i,1}^{(1)} +  \sum_{\tau=1}^{t}\eta^{(\tau)}_i\left(\prod_{j \ne i} w_{j,1}^{(\tau)} \right)  + s_{i,1}^{(t)}$, then we have
\begin{align}
    \xi_i w_{i,1}^{(t+1)}&= \xi_iw_{i,1}^{(1)} +  \xi_i\sum_{\tau=1}^{t}\eta^{(\tau)}_i\left(\prod_{j \ne i} w_{j,1}^{(\tau)} \right)  + \xi_is_{i,1}^{(t)} \nonumber\\
    & = 1 +   \sum_{\tau=1}^{t}\left(\frac{1}{4\sqrt{2\tau \log (nT/\delta)}} \cdot \prod_{j \ne i} \frac{\xi_j w_{j,1}^{(\tau)}}{\|w_j^{(\tau)}\|_1} \right)  + \xi_is_{i,1}^{(t)} \label{eq:1}\\
    & \ge  \frac{1}{2} +  \frac{1}{2\sqrt{2T \log (2nT/\delta)}} \cdot\sum_{\tau=1}^{t}\left( \prod_{j \ne i} \frac{|w_{j,1}^{(\tau)}|}{\|w_j^{(\tau)}\|_1} \right) -
    |s_{i,1}^{(t)}| \label{eq:2}\\
     & \ge \frac{1}{2} + \frac{1}{2\sqrt{2T \log (2nT/\delta)}} \cdot \sum_{\tau=1}^{t}\left(\prod_{j \ne i} \frac{|w_{j,1}^{(\tau)}|}{|w_{j,1}^{(\tau)}| + \frac{3(n'-1)}{2}} \right) \label{eq:3} \\
     & \ge \frac{1}{2} + \frac{1}{2\sqrt{2T \log (2nT/\delta)}} \cdot \sum_{\tau=1}^{t}\left(\prod_{j \ne i} \frac{1}{1 + 3(n'-1)} \right) \label{eq:4}\\
     & \ge \frac{1}{2} \left(1 +  \left(\frac{1}{3n'-2}\right)^{k-1} \frac{t}{\sqrt{2T\log (2nT/\delta)}}\right) \nonumber.
\end{align}
\eqref{eq:1} follows from observing that $\xi_iw_{i,1}^{(1)} = |w_{i,1}^{(1)}| = 1$ and $\prod_{j = 1}^k \xi_j = 1$. \eqref{eq:2} follows from the inductive hypothesis $\xi_jw_{j,1}^{(\tau)} = |\xi_jw_{j,1}^{(\tau)}|$. \eqref{eq:3} follows from our assumption that $|s_{i,1}^{(t)}| \le 1/2$ and Claim \ref{claim:upper}.\eqref{eq:4} follows from the inductive hypothesis $|w_{i,1}^{(\tau)}| \ge 1/2$.
\end{proof}
Setting $T$ such that $T\ge 2\log (2nT/\delta)\alpha^2 (3n'-2)^{2k - 2}$ for $\alpha = 3\sqrt{(n'-1)\log(2k/\epsilon)} - 1$, from the above claims, after iteration $T$, we have
\begin{align*}
w_{i,1}^{(T + 1)} &\ge \frac{\alpha + 1}{2}\\
w_{i,j}^{(T+1)} &\le 3/2 \text{ for }j > 1.
\end{align*}
Using Lemma \ref{lem:err}, we have with probability $1- \delta$,
\[
\mathrm{err}(w_{1:k}^{(T + 1)}) \le 2k \exp\left(-\frac{(\alpha+1)^2}{9(n' -1)}\right) = \epsilon.
\]
\end{proof}

\section{Additional figures, experiments, and discussion}
\label{sec:appendix-experiments}

This section contains our unabridged empirical results, visualizations, and accompanying discussion. Additional example training curves (like the assortment in Figure~\ref{fig:intro} \emph{(left)}) are shown in Figure~\ref{fig:extra-training-curves}; more examples can be found in the subsections below.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/extra-examples.pdf}
    \caption{Additional training curves (on the $(50,3)$-parity problem, except the first row); full details are in Appendix~\ref{subsec:extra-training-curves-details}. \emph{1\textsuperscript{st} row:} The same architecture, initialization, and training algorithm (a width-$100$ ReLU MLP in this case), without an explicit sparse prior, adapts to the computational difficulty parameters $n,k$. \emph{2\textsuperscript{nd} row:} Our positive empirical results hold over a wide range of batch sizes $B$, all the way down to $B=1$. Training is unstable (more outliers) at very large and very small batch sizes. \emph{3\textsuperscript{rd} row:} Even the least overparameterized neural networks, which are barely wide enough to represent parity, converge with reasonable probability (sometimes failing to reach a global minimum). \emph{4\textsuperscript{th} row:} Larger models (width-$1000$ MLPs and $d_\mathrm{emb}=1024$ Transformers) are robust to a wide range of batch sizes. Note the lack of plateaus in setting (vi), which is revisited in Appendix~\ref{subsec:appendix-poly-no-plateau} }
    \label{fig:extra-training-curves}
\end{figure}

\paragraph{Convergence times, success probabilities, and scaling laws.} We first present the full empirical results outlined in Section~\ref{sec:empirics} of the main paper. Figure~\ref{fig:converge-time-table} shows convergence times $t_c$ on small parity instances for all of the architecture configurations enumerated in Section~\ref{subsec:empirics-main}. In some of these settings, $t_c$ exhibits high variance due to unlucky initializations (see Figure~\ref{fig:converge-prob-table}); thus, we report $10^{\text{th}}$ percentile convergence times. Figure~\ref{fig:scaling-table} gives coarse-grained estimates for how $t_c$ scales with $(n, k)$, based on small examples.
For selected architectures, Figure~\ref{fig:extra-scaling-curves} shows how these convergence times scale with $n$ and $k$ more precisely: for small $n$, power law relationships $t_c \propto n^{\Theta(k)}$ are observed for all configurations. Note that for larger $n$, the exponent (i.e. the slope in the log-log plot) increases: with a constant learning rate and standard training, the $n^{\Theta(k)}$ does not continue indefinitely. All additional details are in Appendix~\ref{sec:experiment-details}.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/converge_time_table.pdf}
    \vspace{-3mm}
    \caption{10th percentile convergence times $t_c$ of SGD (with $n = 30$, hinge loss, uniform initialization, and best learning rate $\eta \in \{1, 10^{-1}, 10^{-2}, 10^{-3}\}$), for various architectures, parity degrees $k$, and batch sizes $B$. See Appendix~\ref{sec:experiment-details} for full details.}
    \label{fig:converge-time-table}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/converge_probs.pdf}
    \vspace{-3mm}
    \caption{Percentage of converged runs ($t_c < 10^5$) out of 25 random trials, with the same architectures and training parameters as Figure~\ref{fig:converge-time-table}. With sufficiently large batch sizes, training is extremely robust in settings (i) through (xv). }
    \label{fig:converge-prob-table}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/scaling_table.pdf}
    \caption{Coarse scaling estimates for the 10th percentile convergence time of SGD (with hinge loss, uniform initialization, and best learning rate) on every architecture: $c, \alpha$ such that $t_c \leq c \cdot (n - n_0)^{\alpha}$ on small parity instances $n \in \{10, 20, 30\}, k \in \{2, 3, 4\}$. Missing entries denote cases where $<10\%$ of trials were convergent for any $n$ (see Figure~\ref{fig:converge-prob-table}). Boxes are colored according to $\alpha$. See Appendix~\ref{sec:experiment-details} for full details.}
    \label{fig:scaling-table}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/appendix_scaling_curves.pdf}
    \caption{Finer-grained plots of convergence times for selected architecture configurations, for $n \in \{10, 11, \ldots, 39, 40\}$ and $k \in \{2, 3, \ldots, 8\}$. Medians over $1000$ runs are shown, with $95\%$ bootstrap confidence intervals. \emph{Top row:} standard MLP configurations (i) through (vi). \emph{Bottom row:} Miscellaneous settings ($1$-neuron networks; Transformers; larger $n$). For details on each setting, see Appendix~\ref{subsec:scaling-details}.}
    \label{fig:extra-scaling-curves}
\end{figure}

\paragraph{Guide to this section.} The remainder of Appendix~\ref{sec:appendix-experiments} expands on the various discussions and figures from Sections~\ref{sec:theory} and \ref{sec:discussion}.

\begin{itemize}[leftmargin=*]
    \item Appendix~\ref{subsec:bool} gives experimental evidence that Fourier gaps are present at iterates $w_t$ and initializations $w_0$ other than sign vectors, as well as for activation functions other than ReLU. This suggests that the feature amplification mechanism is robust, and illuminates directions for strengthening the theoretical results.
    \item Appendix~\ref{subsec:building-blocks} discusses how the building blocks of deep learning (activation functions, biases, initializations, learning rates, and batch sizes) play multiple, sometimes conflicting roles in this setting.
    \item Appendix~\ref{subsec:appendix-hidden-progress} provides additional white-box visualizations of hidden progress from Figure~\ref{fig:hidden-progress}.
    \item Appendix~\ref{subsec:appendix-width} explores the implications of the feature amplification mechanism for scaling model size-- namely, unlike random search, large width does not impart parallel speedups.
    \item Appendix~\ref{subsec:appendix-grok} shows that our results hold in the finite-sample setting (allowing for multiple passes over a training set of size $m$). In particular, we show that in low-data regimes, the models exhibit the \emph{grokking} phenomenon.
    \item Appendix~\ref{subsec:noise-experiments} extends our results to \emph{noisy} parities (which comprise the true ``emblematic computationally-hard problem'').
    \item Appendix~\ref{subsec:deep-quad} builds a counterexample for ``layer-by-layer learning'', using parity distributions whose degrees are higher than those of the individual layers' polynomial activations. Preliminary experiments show that standard training works in this setting.
    \item Appendix~\ref{subsec:appendix-poly-no-plateau} presents examples of training curves for wide polynomial-activation MLPs, where, unlike the other settings, there is no initial plateau in the model's error.
\end{itemize}

\subsection{Fourier gaps at initialization and SGD iterates}
\label{subsec:bool}

Proposition~\ref{prop:mlp-from-gap} shows that if the function $x \mapsto \sigma'(w_0^\top x)$ has a Fourier gap at $S$, then $S$ can be identified from a batch gradient at initialization $w_0$ with $B = O(1/\gamma^2)$ samples.
Our end-to-end result (Theorem~\ref{thm:mlp-main}) requires ReLU activations and sign vector initialization, because
the Fourier gap condition (Definition~\ref{def:fourier-gap}) arises from exact formulas for the Fourier coefficients of the majority function.
Stronger end-to-end theoretical guarantees would follow from analogous Fourier gaps in more general population gradients. This requires $x \mapsto \sigma'(w^\top x)$ to satisfy these conditions simultaneously:
\begin{itemize}[leftmargin=*]
    \item \emph{Fourier concentration:} upper bounds on the degree-$(k+1)$ coefficients $\hat f(S \cup \{i\})$, for $i \not \in S$. The term is borrowed from \citet{klivans2004learning}, who use upper bounds on Fourier coefficients of LTFs to approximate them (thus, learn halfspaces) with low-degree polynomials.
    \item \emph{Fourier anti-concentration:} lower bounds on the degree-$(k-1)$ coefficients $\hat f(S \setminus \{i\})$, for $i \in S$.
\end{itemize}

A natural question is:
\emph{which Boolean functions, other than majority, satisfy the $\gamma$-Fourier gap property at $S$, for $\gamma \geq n^{-\Omega(k)}$?}

We present some numerical evidence for large Fourier gaps in functions $x \mapsto \sigma'(w^\top x)$ other than majority, which arise from gradients of architectures other than ReLU MLPs with sign initialization. This shows that the mechanism of feature emergence is empirically robust in settings not fully explained by our current theory.
Establishing corresponding theoretical guarantees would enable stronger end-to-end global convergence guarantees for MLPs and other architectures.

In these experiments, population gradients were computed by brute force integration over all $2^n$ Boolean inputs $x \in \{\pm 1\}^n$. In all cases, for various choices of $\sigma, w$, we measure a slightly weaker notion of Fourier gap $\Gamma$ in the population gradient:
\[ \Gamma_{\sigma,k}(w) := \max_{i \in [k]} | g_i | - \max_{i \not\in [k]} | g_i | ,
\qquad g := \mathop{\mathbb{E}}[y x \, \sigma'(w^\top x)]. \]
If $\Gamma > 0$, then \emph{one}\footnote{Replacing the first $\max$ in the definition of $\Gamma$ with $\min$ would give us a notion closer to Definition~\ref{def:fourier-gap}: if \emph{all} the relevant coordinates are larger than \emph{all} of the irrelevant ones, estimating the population gradient allows us to recover the relevant coordinates.} coordinate from the parity can be identified from $O(1/\Gamma^2)$ samples of the gradient at $w$.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/random_ltf.pdf}
    \caption{Distributions of exact Fourier gaps $\Gamma_{\sigma,k}(w)$ when $\sigma = \mathrm{ReLU}$, and $w$ is an i.i.d. \{uniform, Gaussian\} random vector. These are derived from the Fourier coefficients of the corresponding random LTFs, computed here in exponential time for $n=15$.}
    \label{fig:random-fourier-gaps}
\end{figure}

\paragraph{Random LTFs.}
For ReLU activations and symmetric Bernoulli (i.e. random sign) initialization $w_i \sim \mathrm{Unif}(\{\pm c\}^n)$, the Fourier coefficients are the same as those of majority; thus, there is a Fourier gap of $\gamma \geq n^{-\Omega(k)}$ at every set $S$ (and the same is true of $\Gamma$). We probe the Fourier gaps of linear threshold functions (LTFs) under other ubiquitous initializations: i.i.d. uniform and Gaussian. These are shown in Figure~\ref{fig:random-fourier-gaps}, which indicates (at least for small $n,k$) that the Fourier gap is comparable to that of majority with non-negligible probability.

\paragraph{Random non-LTFs.} The successful convergence of architectures with smoother activations (in the parity setting and beyond) motivates the question of whether large Fourier gaps are present in population gradients corresponding to functions other than LTFs. Figure~\ref{fig:other-fourier-gaps} \emph{(left)} shows that this is the case for sinusoidal activations.

\begin{figure}
    \centering
    \includegraphics[width=0.48\linewidth]{figures/random_act.pdf}
    \includegraphics[width=0.48\linewidth]{figures/sgd_gaps.pdf}
    \caption{Numerical evidence for Fourier gaps beyond random LTFs. \emph{Left:} Fourier gaps of exact population gradients $\mathop{\mathbb{E}}_{\mathcal{D}_S}[y x_i \sin(w^\top x/\sqrt{n})]$ induced by sinusoidal activations, for random initializations $w$. In these small cases, they are comparable the the Fourier gaps of majority (red dashed line). \emph{Right:} Fourier gaps of a sinusoidal neuron's population gradients along the SGD optimization path (10 trials shown). The Fourier gap is consistently positive at initialization, but somewhat smaller than that of majority (red dashed line). Interestingly, it amplifies through the course of training.}
    \label{fig:other-fourier-gaps}
\end{figure}

\paragraph{Boolean functions along the SGD path.} Finally, to further close the gap between Theorem~\ref{thm:mlp-main} and the empirical results, it is necessary to address the fact that SGD accumulates gradients with respect to \emph{time-varying} iterates, while our analysis approximates this using a large-batch gradient at a static iterate $w_0$. In fact, SGD seems to help in some cases: Figure~\ref{fig:other-fourier-gaps} \emph{(right)} shows that when training a sinusoidal neuron, SGD amplifies the initial Fourier gap.

\subsection{Counterintuitive roles of the building blocks of deep learning}
\label{subsec:building-blocks}

Even in this simple problem setting, the simultaneous computational and statistical considerations lead to counterintuitive consequences for the optimal configurations of architectures and algorithms for this setting. We encountered the following, in the search for architecture configurations for the empirical study:

\begin{itemize}[leftmargin=*]
\item \emph{Activation functions.} This mechanism of features emerging via Fourier gaps (see Definition~\ref{def:fourier-gap}) is strongest with non-smooth activations such as the ReLU, whose derivatives are discontinuous threshold functions. This is an orthogonal consideration to representational capacity and mitigation of local minima (under which one might conclude that degree-$k$ polynomial activations are optimal). In summary, in feature learning settings where the Fourier gaps and low-complexity solutions are simultaneous relevant, there is a \emph{sharpness-smoothness tradeoff} for the activation function.
\item \emph{Biases.} The symmetry of the majority function (as well as all unbiased LTFs) causes its even-degree Fourier coefficients to be zero; thus, certain variants of the setups in Section~\ref{sec:empirics} fail for odd $k$. Bias terms (trainable or fixed) are necessary to break this symmetry, in theory and practice. Simultaneously, biases serve the more conventional role of shifting the loss surface; see Section~\ref{subsubsec:width-1-details} for how this affects the details of how the biases were chosen in the experiments.
\item \emph{Initializations.} The role of the initialization distribution is similarly twofold in this setting: $w_0$ should be close to the desired solution $w^*$, but it must also be selected such that SGD will successfully amplify the Fourier gap. A third consideration, which we do not attempt to study in this work, is that multiple randomly-initialized neurons will tend to learn the correct features at different times (see the weight trajectory visualizations in Figure~\ref{fig:hidden-progress} and Figure~\ref{fig:progress-extra}, as well the staircase-like training curves seen for MLPs in Figure~\ref{fig:intro} \emph{(left)}). We expect this \emph{symmetry breaking} phenomenon to be present in more complex feature learning settings. Finally, as shown in the training curves from setting (vi) in Figure~\ref{fig:extra-training-curves}, and in more detail in Appendix~\ref{subsec:appendix-poly-no-plateau}, the choice of activation function influences the qualitative behavior of the training curves: namely, whether the plateaus disappear at large widths and batch sizes.
\item \emph{Batch sizes and learning rates.} The empirical and theoretical results both suggest that SGD uses independent samples to gradually amplify a signal containing the correct features in the initial population gradient of the correlation $\left( \nabla_w \mathop{\mathbb{E}}[ -yx_i \sigma'(w^\top x) ] \right) \vert_{w=w_0}$. However, it would be truer to this mechanism to stay at the initialization $w_0$ until the algorithm has accumulated enough data to discern the correct indices (equivalently, scale up the batch size); in contrast, standard training takes gradients with respect to $w_t$ along the SGD path.
We hypothesize that the bias incurred by this drifting $w_t$ (and thus drifting population gradient) accounts for the degradations seen in Figure~\ref{fig:scaling-laws} \emph{(right)} and Figure~\ref{fig:extra-scaling-curves}. However, Figure~\ref{fig:other-fourier-gaps} \emph{(right)} shows that the movement of SGD can be helpful, amplifying the Fourier gap.
\end{itemize}

\subsection{Hidden progress measures}
\label{subsec:appendix-hidden-progress}

We provide some supplemental plots for the experiments outlined in Figure~\ref{fig:hidden-progress}.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/hidden_progress_sin.pdf}
    \vspace{3mm}
    \includegraphics[width=0.98\linewidth]{figures/hidden_progress_mlp.pdf}
    \caption{Supplementary plots for visualizations of the optimization trajectory $w_t$, and the hidden progress measure $\rho(w_{0:t})$. In the MLP plots, the single scalar $\rho$ is the $\infty$-norm of the entire first layer $W$, and weights are color-coded by row (i.e. neuron). With small batch sizes $B \in \{1,4\}$, per-iteration losses are averaged over a short window (lengths $16$ and $4$, respectively).}
    \label{fig:progress-extra}
\end{figure}

\paragraph{Weight visualizations.} Figure~\ref{fig:hidden-progress} \emph{(left, center)} directly visualizes the evolution of the weights $w_t$ for a single sinusoidal neuron (with a bias, but no second layer). Figure~\ref{fig:progress-extra} supplements these plots from the main paper with additional plots of weight trajectories, at different batch sizes $B$, as well as a width-10 MLP architecture.

\paragraph{$\ell_\infty$ path length.} We present an example of a \emph{single} measurement of the optimization path which captures hidden progress in this setting, which can be plotted alongside loss and accuracy curves. For any iterates of a neuron's weights $w_t \in \RR^n$, we choose the $\infty$-norm of the movement from initialization: $\rho(w_{0:t}) := \norm{w_t - w_0}_\infty$. We present a brief intuitive sketch of the motivation for this choice of $\rho(\cdot)$, and some additional visualizations.

From the theoretical analysis, under the approximation $\nabla \ell(w_t) \approx \nabla \ell(w_0)$ (so that feature learning is performed by estimating the initial population gradient to high precision), we can think of the $i$-th coordinate of $w_t$ as a biased random walk with constant variance $\sigma^2$; the Fourier gap condition entails that biases $\beta_i$ of these random walks are large when $i \in S$. Then, this choice of $\rho$ is an estimate for the drift term $t \cdot \max_i |\beta_i|$, which is larger than the $\sigma \sqrt{t}$ contribution of the variance for sufficiently large $t$.

This progress measure is shown alongside the loss curves in Figure~\ref{fig:progress-extra}, in red. We do not attempt to characterize the dynamics of $\rho$; we only note that they are clearly distinguishable from the maximum of $n$ unbiased random walks, even when SGD appears to make no progress in terms of loss and accuracy.
Studying hidden progress measures in deep learning more quantitatively, as well as in more general settings, presents a fruitful direction for future work.

\subsection{Convergence time vs. width}
\label{subsec:appendix-width}

We provide supplementary plots for the experiment outlined in Figure~\ref{fig:extra-experiments} \emph{(left)}, which probes whether extremely large widths ($r \gg n^k$) afford factor-$r$ parallel speedups of the parity learning mechanism (as one would expect from random search). On 3 parity instances $n \in \{30,40,50\}, k=3$, we varied the width $r \in \{1, 2, 3, \ldots, 9, 10, 30, 100, 300, \ldots, 10^6, 3\times 10^6\}$, keeping all other parameters the same $(B = 128, \eta = 0.1)$.

\paragraph{Results.} We did not find evidence of such parallel speedups over $1000$ runs in each setting; see Figure~\ref{fig:width-extra}. This serves as further evidence that the mechanism by which standard training solves parity is best understood as deterministic and sequential, rather than behaving like random search over size-$k$ subsets. A benefit of width appears to be \emph{variance reduction}: the upper tail of long convergence times is mitigated by a large number of randomly-initialized neurons.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/width_full.pdf}
    \caption{Number of iterations for MLPs (with standard initialization and training) to converge on sparse parity problems, in terms of width $r$. Boxes denote interquartile ranges over 1000 runs; whiskers denote $\pm 1.5 \cdot \mathrm{IQR}$; $\times$ markers denote minimum and maximum outliers. Underparameterized models ($r < k$) are shown in red, and considered ``converged'' at $55\%$ accuracy. Scaling $r$ this way does not lead to $r \times$ parallel speedups, like the expected success time for $r$ copies of random search (shown in green for comparison). }
    \label{fig:width-extra}
\end{figure}

\subsection{Learning and grokking in the finite-sample case}
\label{subsec:appendix-grok}

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/grokking_all.pdf}
    \caption{Supplementary plots for the finite-sample setting. The same configuration (width-$100$ ReLU MLP, $B = 32, \eta = 0.1$) varying the sample size $m$ (decreasing from left to right) and weight decay $\lambda$ (increasing from top to bottom). When $m$ is sufficiently large (much larger than the statistical threshold $\Theta(k \log n)$), generalization error is negligible. When $m$ is too small, the model fails to train. In between, we observe eventual convergence to the correct solution, with training curves exhibiting the \emph{grokking} phenomenon. Weight decay governs a statistical-computational tradeoff in this setting: larger $\lambda$ improves generalization, but can cause optimization to fail (bottom row).}
    \label{fig:grok-extra}
\end{figure}

We provide some supplementary plots for the experiments outlined in Figure~\ref{fig:extra-experiments} \emph{(right)}. In these settings, a fixed architecture (width-$100$ MLP with ReLU activations) is trained with minibatch SGD in an otherwise fixed configuration (hinge loss, learning rate $\eta = 0.1$, batches of size $B=32$) on a finite training sample of size $m$. We also vary a \emph{weight decay} parameter $\lambda$.

As shown in Figure~\ref{fig:grok-extra}, the weight decay parameter $\lambda$ modulates a delicate computational-statistical tradeoff: it improves generalization (expanding the range of $m$ for which training eventually finds the correct solution), but the model fails to train at large values of $\lambda$. For small $m$ and appropriately tuned $\lambda$, we observe \emph{grokking}: the model initially overfits the training data, but finds a classifier that generalizes after a large number of iterations.

\subsection{Learning noisy parities}
\label{subsec:noise-experiments}

The other empirical results in this work focus on noiseless parity distributions $\mathcal{D}_S$, to reduce the number of sources of variance and degrees of freedom. However, the setting of \emph{random classification noise} is important for several reasons. In this section, we briefly demonstrate that our results extend to this case. Let $\mathcal{D}_S^{(\epsilon)}$ denote the \emph{$(n,k,\epsilon)$-noisy parity distribution}, defined by flipping the labels in the $(n,k)$-parity distribution $\mathcal{D}_S$ independently with probability $\frac{1}{2} - \epsilon$. Note that when $\epsilon = 0$, the labels are completely random (thus, $S$ cannot be learned). By a standard PAC-learning argument, when $0 < \epsilon \leq \frac{1}{2}$, the statistical limit for identifying $S$ from i.i.d. samples from $\mathcal{D}_S$ scales as
$\Theta\left( \frac{k \log n}{\epsilon^2} \right)$.

\paragraph{Motivations.} First, learning parities from \emph{noisy} samples is the true ``emblematic computationally-hard distribution''. Without noise, there is a non-SQ algorithm which avoids the exponential-in-$k$ computational barrier: Gaussian elimination can identify $S$ in $O(n^3)$ time and $\Theta(n)$ samples. Second, viewing parities as an idealized setting in which to understand training dynamics, resource scaling, and emergence in deep learning, it is important to see that this phenomenon is robust to label noise.

\paragraph{Theory.} It is easy to incorporate label noise into the theoretical analysis, which works with correlations of the form $\mathop{\mathbb{E}}_{\mathcal{D}_S}[y \, f(x)]$; each coordinate of the population gradient of the correlation loss is a quantity of this form. In the noisy case, these quantities are replaced with
\[\mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S^{(\epsilon)}}[y \, f(x)] = \epsilon \cdot \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S}[y \, f(x)].\]
In particular, when architecture's population gradient has a Fourier gap with parameter $\gamma$ in the noiseless case implies a Fourier gap with parameter $\epsilon \cdot \gamma$.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/noisy.pdf}
    \caption{Training curves (over $5$ random seeds) for a width-$100$ ReLU MLP on noisy $(50,3,\epsilon)$-parity learning problems at various batch sizes $B$ and noise levels (flipping labels with probability $p$, so that $\epsilon = 1 - 2p$). The models learn the features and converge successfully (measured by accuracy on the noiseless distribution), even with $49\%$ of labels flipped randomly (i.e. $\epsilon = 0.01$). This is a preliminary illustration that the phenomena investigated in this paper are robust with respect to i.i.d. label noise. Note that the scale of $t$ is much larger for small batches and high noise. }
    \label{fig:noisy}
\end{figure}

\paragraph{Experiments.} We find that the experimental findings are robust to label noise, in the sense that models are able to obtain nontrivial (and sometimes $100\%$) accuracy;
see Figure~\ref{fig:noisy} for some training curves under various settings of $\epsilon$. This provides concrete evidence against the (already extremely dubious) hypothesis that neural networks, with standard initialization and training, learn noiseless parities by implicitly simulating an efficient algorithm such as Gaussian elimination. Note that with a constant learning rate (here, $\eta = 0.1$) and label noise, the iterates of SGD do not always converge to $100\%$ accurate solutions.

\subsection{Counterexample for layer-by-layer learning}
\label{subsec:deep-quad}

\paragraph{Notation.} Consider an $L$-layer MLP with activation $\sigma$, parameterized by weights and biases
\[\theta = (W_1, b_1, \ldots, W_{L-1}, b_{L-1}, u),\]
and defined by
\[ f_\mathsf{mlp}(x; \theta) := (f_1 \circ \ldots \circ f_{L-1} \circ f_L)(x), \]
where $f_i$ denotes the function $z \mapsto \sigma(W_i z + b_i)$ for $1 \leq i \leq L-1$, and $f_L$ denotes $z \mapsto u^\top z$.
The shapes of the parameters $W_i, b_i, u$ are selected such that each function composition is well-defined.
Let the \emph{intermediate activations} at layer $i$ be denoted by
\[ z_i(x; \theta) := (f_1 \circ \ldots \circ f_{i})(x). \]
Finally, $r_i$ (the \emph{width at layer $i$}) refers to the dimensionality of $z_i$ as defined above.

\paragraph{Construction where layer-by-layer learning is impossible.} Notice that when $\sigma$ is a degree-$2$ polynomial (say, $\sigma(z) = z^2$), an $L$-layer MLP can represent parities up to degree $2^{L-1}$-- for example, a $3$-layer MLP (which composes quadratic activations twice) can represent a $4$-sparse parity as a $2$-sparse parity of $2$-sparse parities. However, Equation~(\ref{eq:parity-orth}) implies the following:
\begin{itemize}[leftmargin=*]
    \item An individual layer cannot represent a parity of $k > 2$ inputs.
    \item The population gradient (as in Equation~(\ref{eq:pop-grad}) is zero (since every coordinate of the gradient is the correlation between a $k$-wise parity and a polynomial of degree $2$).
\end{itemize}
Thus, this setting serves as an idealized counterexample for layer-by-layer learning: if SGD succeeds on parities with higher degree than the architecture's polynomial activations, it must do so by an end-to-end mechanism. Intuitively, earlier layers can only make progress by \emph{knowing how their outputs will be used downstream}. Concretely, consider the population gradient of the correlation loss, with respect to a first-layer neuron's weights $w := (W_1)_{j,:}$. With layer-by-layer training, this gradient contains no information:
\[ \nabla_w \mathop{\mathbb{E}}[ -y x_i \, \underbrace{\sigma'(w^\top x)}_{\text{degree $1$}} \, u_j ] = 0.\]
However, in end-to-end training, the presence of downstream layers removes this barrier:
\[ \nabla_w \mathop{\mathbb{E}}\Big[ -y x_i \, \underbrace{\sigma'(w^\top x)}_{\text{degree $1$}} \, \underbrace{\frac{\partial f_2 \circ \ldots \circ f_L}{\partial (z_1)_j}}_{\text{degree $2^{L-2} - 1$}} \Big], \]
giving the gradient greater representation capacity (in terms of polynomial degree). The question remains of whether end-to-end training works in this setting, which we resolve positively in small experiments.

\paragraph{Results: end-to-end training works empirically.} We empirically observed successful training (to $100\%$ accuracy) in a few settings (with SGD, learning rate $\eta = 0.01$, batch size $B=32$, and default uniform initialization as described in Appendix~\ref{subsec:arch-configs}):
\begin{itemize}[leftmargin=*]
\item $L=3, n \in \{10, 20, 30\}, k \in \{1, 2, 3, 4\}$. Small widths suffice: $(r_1, r_2) = (2, 1)$. Over 10 random seeds, all models converged within $20000$ iterations.
\item $L=4, n \in \{10, 20, 30\}, k \in \{1, 2, 3, 4, 5, 6\}$. Widths were chosen to be slightly larger for stability: $(r_1, r_2, r_3) = (10, 10, 1)$. Over 10 random seeds, all models converged within $50000$ iterations. Additionally, models trained on $(n,k) \in \{ (10,7), (20, 7), (30, 7), (10, 8) \}$ converged within $500000$ iterations.
\end{itemize}
As a sanity check, the models failed to converge in experimental setups where $k > 2^{L-1}$: $(L = 2, k \geq 3)$ and $(L = 3, k \geq 5)$.

\paragraph{Discussion.} This construction serves as a simple counterexample to the \emph{``deep only works if shallow is good''} principle of \cite{malach2019deeper}, demonstrating a case where a deep network can get near-perfect accuracy even when greedy layerwise training (e.g. \citep{belilovsky2019greedy}) cannot beat trivial performance. It remains to characterize these positive empirical results theoretically, as well as to investigate whether there are pertinent analogues in real data distributions.

\subsection{Lack of plateaus for wide polynomial-activation MLPs}
\label{subsec:appendix-poly-no-plateau}

An interesting qualitative observation from the training curves in Figure~\ref{fig:extra-training-curves} is that the validation accuracy curves in setting (vi) (width-1000 polynomial-activation MLPs) do not follow the same ``plateau'' or ''staircase'' pattern as the others. Figure~\ref{fig:poly-no-plateau} shows a few additional examples of training curves for polynomial-activation MLPs, varying the width $r$ and batch size $B$. We find that the rate of descent of the validation error increases with both of these parameters; note that this does not occur with ReLU activations (where there are sharp phase transitions between plateaus at all batch sizes).

This constitutes an exception to this paper's theme of ``hidden progress'' behind flat loss (or error) curves: with enough overparameterization \emph{and} ``over-sampling'', the continuous progress of SGD in this setting is no longer hidden, and manifests in the training curves. This phenomenon seems to be specific to certain activation functions; we leave it for future work to understand why and when it occurs, as well as potential practical implications.

\begin{figure}
    \centering
    \includegraphics[width=0.98\linewidth]{figures/poly_no_plateau.pdf}
    \caption{Additional training curves for polynomial-activation MLP architectures (iv), (v), (vi), on the (50, 3)-sparse parity problem. Unlike the other architectures, these settings exhibit continuous progress when the width $r$ and batch size $B$ are large.}
    \label{fig:poly-no-plateau}
\end{figure}
\section{Details for all experiments}
\label{sec:experiment-details}

\subsection{Deep learning configurations}
\label{subsec:arch-configs}

\paragraph{Losses.}
Our ``robust space'' of empirical results use the following loss functions:
\begin{itemize}[leftmargin=*]
    \item Hinge: $\ell(y, \hat y) := (1 - y\hat y)_+$.
    \item Square: $\ell(y, \hat y) := (y - \hat y)^2$.
    \item Cross entropy: $\ell(y, \hat y) := - \log \frac{e^{y \hat y}}{1 + e^{y \hat y}}$.
\end{itemize}
Additionally, the theoretical analysis considers the \emph{correlation loss} $\ell(y, \hat y) := - y\hat y$.

In the configurations corresponding to all of the figures and convergence time experiments, we used the hinge loss. This was a relatively arbitrary choice (i.e. they appeared to be interchangeable upon running small experiments); an advantage of the hinge and square losses over cross entropy is that for architectures that can realize the parity function, there is a zero-loss solution with finite weights.

\paragraph{Initializations.}
Our empirical results use the following i.i.d. weight initializations:
\begin{itemize}[leftmargin=*]
    \item Uniform on the interval $[-c, c]$, where the scale $c$ is chosen for all affine transformation parameters using the ``Xavier initialization'' convention \citep{glorot2010understanding}. The experiments are quite tolerant to the particular choice of $c$ (as these are not deep networks); this choice, which is the default in deep learning packages, emphasizes that our positive empirical results hold under a \emph{standard} initialization scheme.
    \item Gaussian with mean $0$ and variance $\sigma^2$, selected using the ``Kaiming initialization'' convention \citep{he2015delving}.
    \item Bernoulli (i.e. random sign) initialization: the discrete distribution $\mathrm{Unif}(\{-c, c\})$, for the same choice of $c$ as for the uniform distribution.
\end{itemize}

\subsubsection{2-layer MLPs}
We consider $2$-layer MLPs $f(x; W, b, u) = u^\top \sigma(Wx + b)$ for two choices of activations:
\begin{itemize}[leftmargin=*]
    \item ReLU: $\sigma(z) := (z)_+ = \max(0,z)$.
    \item Degree-$k$ polynomial: $\sigma(z) := z^k$.
\end{itemize}
In both cases, whenever $r \geq k$ (and, in the case of polynomial activations, choosing the degree to be $k$), there exists a width-$r$ MLP which can represent $k$-sparse parities: for all $(n,k)$ and $|S|=k$, there is a setting of $W, b, u$ such that $f(x; W, b, u) = \chi_S(x)$.

Note that if the output $f(x; \theta)$ is a degree-$k' < k$ polynomial in $x$ (e.g. an MLP with $\sigma(z) = z^{k'}$ activations), the architecture is incapable of representing a parity of $k$ inputs. In fact, it is incapable of representing any function that has a nonzero correlation with parity; this follows from orthogonality (Equation (\ref{eq:parity-orth})).

\subsubsection{Single neurons}
\label{subsubsec:width-1-details}

To explore the limits of concise parameterization for architectures capable of learning parities, we propose a variety of non-standard activation functions which allow a single neuron to learn sparse parities. These constructions leverage the fact that the parity is a nonlinear function of the sum of its inputs $w_S^\top x$, where $w_S := \sum_{i \in S} e_i$.

\begin{figure}
    \centering
    \includegraphics[width=0.16\linewidth]{figures/activations-left.pdf}
    \includegraphics[width=0.64\linewidth]{figures/activations-center.pdf}
    \includegraphics[width=0.16\linewidth]{figures/activations-right.pdf}
    \caption{Visualizations of all activation functions considered in this work. As discussed in Section~\ref{subsubsec:width-1-details}, care must be taken to ensure that the architectures can realize sparse parities with the same bias $b$ (in particular, $b = 0$) across varying $k$. Multiple $k$-dependent displacements are shown for the $\infty$-zigzag and sinusoidal activations.}
    \label{fig:funny-activations}
\end{figure}

\paragraph{$k$-zigzag activation.} $\sigma(\cdot)$ is the piecewise linear function which interpolates the $k+1$ points $\{(k,+1), (k-2,-1), (k-4,+1), \ldots, (-k+2, \pm 1), (-k, \mp 1) \}$ with $k$ linear regions $\{(-\infty,-k], [-k,-k+2], \ldots, [k-2,k], [k, +\infty) \}$. Then, $\sigma(w_S^\top x) = \chi_S(x)$.
\paragraph{Oscillating polynomial activation.} $\sigma(\cdot)$ is the degree-$k$ polynomial which interpolates the same points as above.
\paragraph{$\infty$-zigzag activation.} The \emph{infinite} extension of the zigzag activation is the \emph{triangle wave} function $\sigma(\cdot)$ which linearly interpolates the infinite set of points $\bigcup_{i \in \ZZ} \{ (2i, +1), (2i+1, -1) \}$. This can express parities of \emph{arbitrary} degree. The $+1$ and $-1$ can be swapped (resulting in an activation which is equivalent when shifted by a bias term). However, in our experiments, we choose the sign convention depending on $k$ such that $\sigma(w_S^\top x) = +1$. This allows different convergence time curves to be more directly comparable across different $k$, since it removes the effects of the bias of the global minimizer alternating with $k$.
\paragraph{Sinusoidal activation.} $\sigma(z) := \sin(z)$. The sinusoidal neuron $\sin(w^\top x + b)$ can also express parities of arbitrary degree, since it can interpolate the same set of points as the $\infty$-zigzag activation. In the experiments, we pick a shift $\beta$ and use the activation $\sigma(z) := \sin(\frac{\pi}{2} z + \beta)$, such that $\sigma(z)$ interpolates the same points as the sign convention selected for the $\infty$-zigzag activation. In the experiments in Section~\ref{fig:hidden-progress}, the sinusoidal activation is additionally scaled by a factor of 2 ($z \mapsto \sigma(2z)$); this is interchangeable with scaling the learning rate and initialization, and is done to obtain more robust convergence in the particular setting of $(n,k)=(50,3)$.

Figure~\ref{fig:funny-activations} visualizes all families of activations considered in this paper.

\subsubsection{Parity Transformer}
\label{subsubsec:parity-transformer}

The Transformer experiments use a slightly simplified version of the architecture introduced in \citep{vaswani2017attention}. In particular, it omits dropout, layer normalization, tied input/output embedding weights, and a positional embedding on the special $\texttt{[CLS]}$ token. Including these does not change the results significantly; they are all present in the preliminary findings in \citep{edelman2021inductive} (in which an ``off-the-shelf'' Transformer implementation successfully learns sparse parities). We specify the architecture below.

Our \emph{Parity Transformer} has the following hyperparameters: sequence length $n$, token embedding dimension $d_\mathrm{emb}$, attention embedding dimension $d_\mathrm{attn}$, feedforward embedding dimension $d_\mathrm{mlp}$, and number of heads $H$. Its trainable parameters (together denoted by $\theta$) are:

\begin{itemize}[leftmargin=*]
    \item Token embeddings $E_{-1}, E_{+1}, E_{\texttt{[CLS]}} \in \RR^{d_\mathrm{emb}}$ and position embeddings $P_1, \ldots, P_n \in \RR^{d_\mathrm{emb}}$. Let $\theta_\mathsf{emb}$ denote this subset of parameters.
    \item Attention head matrices: $W_Q^{[h]}, W_K^{[h]}, W_V^{[h]} \in \RR^{d_\mathrm{emb} \times d_\mathrm{attn}}$ and $W_\mathrm{out}^{[h]} \in \RR^{d_\mathrm{attn} \times d_\mathrm{emb}}$, for $h = 1, \ldots, H$. Let $\theta_\mathsf{attn}$ denote this subset of parameters.
    \item MLP weights and biases: $W_1 \in \RR^{d_\mathrm{mlp} \times d_\mathrm{emb}}, b_1 \in \RR^{d_\mathrm{mlp}}, W_2 \in \RR^{d_\mathrm{mlp} \times d_\mathrm{emb}}, b_2 \in d_\mathrm{emb}$. Let $\theta_\mathsf{mlp}$ denote this subset of parameters.
    \item Classification head: $u \in \RR^{d_\mathrm{emb}}$.
\end{itemize}

Then,
\[ f_\mathsf{tf}(x ; \theta) := u^\top \left( f_\mathsf{mlp}( f_\mathsf{attn}( f_\mathsf{emb}(x; \theta_\mathsf{emb}) ; \theta_\mathsf{attn} ) ; \theta_\mathsf{mlp} ) \right), \]
where these submodules are defined by

\begin{itemize}[leftmargin=*]
    \item Embedding $f_\mathsf{emb} : \{\pm 1\}^n \rightarrow \RR^{(n+1) \times d_\mathrm{emb}}$: $f_\mathsf{emb}(x ; \theta_\mathsf{emb})_{i,:} := E_{x_i} + P_i$ for $i \in [n]$. We will include an the extra index $\texttt{[CLS]}$, for which $f_\mathsf{emb}(x ; \theta_\mathsf{emb})_{\texttt{[CLS]},:} := E_{\texttt{[CLS]}}$ (with no positional embedding). \texttt{[CLS]} stands for ``classification'', as in ``use the output at this position to classify the sequence''. This is a standard construction which makes the classifier permutation-invariant.
    \item Attention block $f_\mathsf{attn} : \RR^{(n+1) \times d_\mathrm{emb}} \rightarrow \RR^{d_\mathrm{emb}}$:
    \[f_\mathsf{attn}(X ; \theta_\mathsf{attn}) = X_{\texttt{[CLS]}, :} + \sum_{h=1}^H \mathsf{softmax}\left( \frac{1}{\sqrt{d_\mathrm{attn}}} X_{\texttt{[CLS]},:}W_Q^{[h]} (X W_K^{[h]})^\top \right) W_V^{[h]} W_\mathrm{out}^{[h]},\]
    where $\mathsf{softmax}(z) := \exp(z) / \mathbf{1}^\top\exp(z)$. Note that we have specialized this architecture to a single output, at the \texttt{[CLS]} position.
    \item MLP $f_\mathsf{mlp} : \RR^{d_\mathrm{emb}} \rightarrow \RR^{d_\mathrm{emb}}$:
    \[f_\mathsf{mlp}(z ; \theta_\mathsf{mlp}) := z + W_2 \sigma(W_1x + b_1) + b_2,\]
    where $\sigma(\cdot) = \mathsf{GeLU}(\cdot)$ (the Gaussian error linear unit) is the standard choice in Transformers.
\end{itemize}

\paragraph{Training.} Each matrix-shaped parameter was initialized using PyTorch's default ``Xavier uniform'' convention. Unlike the other settings considered in this paper, we were unable to observe successful convergence beyond a few small $(n,k)$ using standard SGD. As is common practice when training Transformers, we used Adam \citep{kingma2014adam} with default adaptive parameters $\beta_1 = 0.9, \beta_2 = 0.999, \epsilon = 10^{-8}$ in our experiments. While there are more fine-grained accounts of why Adam outperforms vanilla SGD \citep{zhang2020adaptive,agarwal2020disentangling}, finding the optimal optimizer configuration and investigating ablations of this optimizer are outside the scope of this work. In this work, we only tune Adam's learning rate $\eta$.

\subsubsection{PolyNet}
For positive integers $n, k$, the PolyNet architecture is parameterized by weights and biases $\theta := \{(w_i \in \RR^n, b_i \in \RR)\}_{i=1}^k$, and is defined by
\[f_\mathsf{PolyNet}(x; \theta) := \prod_{i=1}^k (w_i^\top x + b_i). \]
Even with all biases $b_i$ set to $0$, this architecture can realize a $k$-wise parity, by setting $\{w_i\} = \{e_j : j \in S\}$ in any permutation.

\subsubsection{Details for Figure~\ref{fig:intro} (left)}
\label{subsec:fig1-details}

Figure~\ref{fig:intro} (left) shows training curves from 8 representative configurations, with online i.i.d. samples from the same distribution, corresponding to the $(n=50, k=3)$-sparse parity problem. The first row encompasses various MLP settings with standard activations:

\begin{itemize}[leftmargin=*]
    \item Setting (i): width-$10$ MLP with ReLU activation ($B = 32, \eta = 0.5$).
    \item Setting (i): width-$10$ MLP with ReLU activation, with large batches ($B = 1024, \eta = 0.05$).
    \item Setting (ii): width-$100$ MLP with ReLU activation, with tiny batches ($B = 1, \eta = 0.05$).
    \item Setting (iv): width-$10$ MLP with polynomial $\sigma(z) = z^3$ activation ($B = 32, \eta = 0.05$).
\end{itemize}

The second row shows other settings:
\begin{itemize}[leftmargin=*]
    \item Setting (vii): width-$1$ MLP with a piecewise linear $k$-zigzag activation ($B = 32, \eta = 0.2$).
    \item Setting (x): width-$1$ MLP with a sinusoidal activation (scaled and shifted for $k = 3$; see the discussion in Section~\ref{subsubsec:width-1-details}) ($B = 32, \eta = 0.05$).
    \item Setting (*iii): Parity Transformer, with $d_\mathrm{emb} = 1024, d_\mathrm{attn} = 8, H = 128$ ($B = 32, \eta = 5 \times 10^{-4}$).
    \item Setting (xv): degree-$3$ PolyNet ($B = 32, \eta = 0.07$).
\end{itemize}

\subsubsection{Details for Figure~\ref{fig:extra-training-curves}}
\label{subsec:extra-training-curves-details}

The first row uses the width-$10$ ReLU MLP configuration $(ii)$, holding $B = 32$ and $\eta = 0.1$ while varying the task difficulty across 6 settings: $(n,k) \in \{(30,3),(60,3),(90,3),(30,4),(30,5),(30,6)\}$. The remaining plots are all for the $(50,3)$ setting.

The second row uses the $k=3$ PolyNet configuration $(xv)$, varying $(B,\eta) \in \{(1,0.005), (4,0.01), (16,0.1), \allowbreak (64,0.1), (256,0.1), (1024,0.1)\}$.

The third row uses the minimally-wide configurations (*i), (*ii), (vii), (viii), (xi), (xii) (thus presenting an example for each non-standard activation), holding batch size $B=1$. $\eta = 0.1$ in each of the cases except (*ii), where $\eta = 0.01$.

The fourth row uses three large architectures: settings (iii), (vi), and (*iii), with $(B,\eta) \in \{(1,0.1), (1024,0.1), \allowbreak (1,0.001), (1024,0.01), (32,0.0003), (1024,0.0003).\}$ (*iii) uses the Adam optimizer instead of SGD.

\subsubsection{Details for Figure~\ref{fig:extra-scaling-curves}}
\label{subsec:scaling-details}

Figure~\ref{fig:extra-scaling-curves} contains scaling plots in various settings for the median convergence time $t_c$. Below, we give comprehensive details about these settings. For each of these runs, we chose $B = 32$ (settings with smaller batch sizes exhibited additional variance; with larger batch sizes, the models were slower to converge), as well as the hinge loss. We used SGD with constant learning rate $\eta$ (enumerated below), except in setting (*iii).

The top row shows MLP settings (i) through (vi). From left to right:
\begin{itemize}[leftmargin=*]
    \item Setting (i): width-$10$ MLP with ReLU activation ($\eta = 1$).
    \item Setting (ii): width-$100$ MLP with ReLU activation ($\eta = 1$).
    \item Setting (iii): width-$1000$ MLP with ReLU activation ($\eta = 1$).
    \item Setting (iv): width-$10$ MLP with $\sigma(z) = z^k$ activation ($\eta = 0.01$).
    \item Setting (v): width-$100$ MLP with $\sigma(z) = z^k$ activation ($\eta = 0.01$).
    \item Setting (vi): width-$1000$ MLP with $\sigma(z) = z^k$ activation ($\eta = 0.01$).
\end{itemize}

The bottom row shows miscellaneous settings. From left to right:
\begin{itemize}[leftmargin=*]
    \item Setting (vii): width-$1$ MLP with degree-$k$ oscillating polynomial activation interpolating the parity function ($\eta = 0.01$).
    \item Setting (xiv): single sinusoidal neuron with no second layer ($\eta = 0.01$).
    \item Setting (*iii): Parity Transformer, with $d_\mathrm{emb} = 1024, d_\mathrm{attn} = 8, H = 128$ ($B = 32, \eta = 3 \times 10^{-4}$).
    \item Setting (iv): degree-$k$ PolyNet ($\eta = 0.05$).
    \item Setting (ii): width-$100$ MLP with ReLU activation ($\eta = 1$), showing an expanded range of $n$ for smaller $k$.
    \item Setting (xv): width-$100$ MLP with $\sigma(z) = z^k$ activation ($\eta = 1$), showing an expanded range of $n$ for smaller $k$.
\end{itemize}

\subsection{Training curves and convergence time plots}
For all example training curves in all figures (in Sections~\ref{sec:empirics} and \ref{sec:discussion}, as well as the appendix), population losses and accuracies are approximated using a batch of size $8192$, sampled once at the beginning of training from the same distribution $\mathcal{D}_S$. All plots of single representative training runs use a fixed random seed (\texttt{torch.manual\_seed(0)}); when $R$ training runs are shown, seeds $0, \ldots, R-1$ are used.

In Figures~\ref{fig:converge-time-table} and \ref{fig:converge-prob-table}, validation accuracies were recorded every $10$ iterations, and a run was recorded as converged if it reached $100\%$ accuracy within $10^5$ iterations; we report the $10^\text{th}$ percentile over 25 random seeds, to reduce variance arising from the more initialization-sensitive settings. In Figure~\ref{fig:scaling-table}, coarse-grained scaling estimates for the ($10^\text{th}$ percentile) convergence time are computed as follows: for $n \in \mathcal{N} := \{10, 20, 30\}$, the smallest $\alpha$ is chosen such that $t_c \leq c \cdot (n - n_0)^{\alpha}$, choosing $n_0 = \min \mathcal{N} - 1 = 9$, so that $c = t_c$ at $n = 10$. These estimates are calculated to give quantitative order-of-magnitude upper bounds for the convergence time. Indeed, the power-law convergence times do not extrapolate at a constant learning rate; see Figure~\ref{fig:scaling-laws} \emph{(right)}, the ``larger $n$'' plots in Figure~\ref{fig:extra-scaling-curves}, and the discussion on batch sizes and learning rates in Appendix~\ref{subsec:building-blocks}.

To reduce computational load, for the larger-scale probes of convergence times $t_c$, validation accuracies were instead checked on a sample of size $128$. For the underparameterized networks (i.e. unable to represent parity, but can still get a meaningful gradient signal), this threshold was changed to $10$ consecutive batches with accuracy at least $55\%$.
Note that for parity learning in particular, a weak learner can be converted into a strong learner: there is an efficient algorithm \citep{goldreich1989hard,kushilevitz1993learning} which, given a classifier which achieves $1/2 + \epsilon$ accuracy on $\mathcal{D}_S$ for a constant $\epsilon > 0$, outputs $S$ with high probability.

In the median convergence time plots in Figure~\ref{fig:intro} \emph{(right)}, Figure~\ref{fig:scaling-laws} \emph{(right)}, and Figure~\ref{fig:extra-scaling-curves}, error bars for median convergence times in all plots are $95\%$ confidence intervals, computed from $100$ bootstrap samples. Each point on the each curve corresponds to $1000$ random trials. Halted curves signify more than $50\%$ of runs failing to converge within $T=10^5$ iterations (hence, infinite medians).

\subsection{Implementation, hardware, and compute time}
All training experiments were implemented using PyTorch \citep{paszke2019pytorch}.

Although most of the networks considered in the main empirical results are relatively small, a large ($\sim 10^8$) total number of models were trained to certify the ``robust space'' of results and obtain precise scaling curves. These individual experiments were not large enough to benefit from GPU acceleration; on an internal cluster, the CPU compute expenditure totaled approximately $1500$ CPU hours.

A subset of these experiments stood to benefit from GPU acceleration: width $r \geq 100$ MLPs; scaling behaviors for $n \geq 100$; all experiments involving Transformers. These were performed with NVIDIA Tesla P100, Tesla P40, and RTX A6000 GPUs on an internal cluster, consuming a total of approximately $200$ GPU hours.
\section{Additional background, preliminaries, and related work} \label{app:preliminaries}

\subsection{Parities, orthogonality, and hardness}
For each integer $n \geq 1$ and nonempty subset of indices $S \subseteq [n]$, define the $(n,S)$-\emph{parity function} $\chi_S(x) = \prod_{i \in S} x_i$, i.e. the parity of the bits at the indices in $S$. We define the $(n,S)$-\emph{parity distribution} $\mathcal{D}_S$ over examples $(x,y)$ as follows: the features $x \sim \mathrm{Unif}(\{\pm 1\}^n)$ are uniform random bits, with labels $y \in \{\pm 1\}$ given by the parity function $y = \chi_S(x)$. For $0 \leq k \leq n$, the corresponding $(n,k)$-\emph{parity learning problem} is the task of identifying an unknown size-$k$ set $S$ (chosen at random), using samples from $\mathcal{D}_S$. With knowledge of $k$ but not $S$, a learner must use the labels to distinguish between $\binom{n}{k}$ possible ``relevant feature sets''; thus, the statistical limit for this problem is $\log \binom{n}{k} = \Theta(k \log n)$ samples.

This work leverages the parity problem as a ``computationally hard case'' for identifying the features $S$ which are relevant to the label. Observe that for any $S' \in [n]$, it holds that
\begin{equation}
\mathop{\mathbb{E}}_{x \sim \mathrm{Unif}(\{\pm 1\}^n)} \left[\chi_S(x) \chi_{S'}(x)\right] = \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S} \left[ \chi_{S'}(x) \, y \right] =
\begin{cases}
1 & S' = S \\
0 & \text{otherwise}
\end{cases}.
\label{eq:parity-orth}
\end{equation}
That is, a learner who guesses indices $S'$ cannot use correlations as feedback to reveal which (or how many) indices in $S'$ are correct, unless $S'$ is \emph{exactly} the correct subset. In this sense, for the $(n,k)$-parity problem, the $\binom{n}{k} - 1$ wrong answers are indistinguishable from each other. Thus, the structure of this problem forces this form of learner (but not necessarily all learning algorithms) to perform exhaustive search over subsets.

Property~\ref{eq:parity-orth} (a.k.a. the orthogonality of parities under the correlation inner product) implies that any function $f : \{\pm 1\}^n \rightarrow \RR$ has a unique Fourier expansion (see, e.g. \citep{o2014analysis}):
\begin{equation}
f(x) = \sum_{S \subseteq [n]} \hat f(S) \chi_S(x), \qquad \hat f_S = \mathop{\mathbb{E}}_{x \sim \mathrm{Unif}(\{\pm 1\}^n)} \left[f(x) \chi_S(x) \right].
\label{eq:fourier}
\end{equation}
In the statistical query (SQ) model \citep{kearns1998efficient}, Property~(\ref{eq:parity-orth}) implies computational lower bounds. In this model, a learner, rather than having access to examples drawn from the distribution, can query an oracle, which responds with noisy estimates of the query over the distribution. Namely, each iteration the learner outputs a query $q_i : \{\pm 1\}^n \to [-1,1]$, and the oracle returns some value $v_i$ satisfying $|v_i - \mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}}\left[q_i(x,y)\right]| \le \tau$, for some tolerance parameter $\tau > 0$. It can be shown that Equation~\ref{eq:parity-orth} implies that each query will have a non-trivial correlation only with a small fraction (namely, $1/\tau^2$) of the possible parities, which then implies that an SQ algorithm which solves the $(n,k)$-parity problem using $T$ queries of tolerance $\tau$ must satisfy $T / \tau^2 \ge \Omega(n^k)$. This constitutes a lower bound on the number of queries (and/or on the tolerance), which indicates that essentially, SQ algorithms cannot do much better than exhaustive search (going over all the possible choices of size-$k$ subsets).

It should be mentioned that the $(n,k)$-parity problem can be solved efficiently by a learning algorithm that has access to examples (i.e., an algorithm that does not operate in the SQ framework). Specifically, this problem can be solved by the Gaussian elimination algorithm. Moreover, it has been shown that the (Stochastic) Gradient Descent algorithm, discussed in the next section, can also be utilized for solving parities, given accurate enough estimates of the gradient and a very particular choice of neural network architecture \cite{abbe2020poly}. That said, when the accuracy of the gradients is not sufficient, GD suffers from the same SQ lower bound mentioned above (i.e., GD is essentially an SQ algorithm \cite{abbe2021power}).

Learning sparse \emph{noisy} parities, even at a very small noise level (i.e., $o(1)$ or $n^{-\delta}$) is believed to be computationally hard. This was first explicitly conjectured by \citet{alekhnovich2003more}, and has been the basis for several cryptographic schemes (e.g., \citep{ishai2008cryptography, applebaum2009fast, applebaum2010public}). For noiseless sparse parities, \citet{kol2017time} show time-space hardness in the setting where $k = \omega(1)$. We present some experiments with noisy parities in Appendix~\ref{subsec:noise-experiments}, finding that our empirical results (and theoretical analysis) are robust to $\Theta(1)$ noise.

\subsection{Neural networks and standard training}
\label{subsec:prelims-nn}

Next, we establish notation for the standard neural network training pipeline. Our main results are presented in the online learning setting, with a stream of i.i.d. batches of examples. At each iteration $t = 1, \ldots, T$, a learning algorithm receives a batch of $B$ examples $\{(x_{t,i},y_{t,i})\}_{i = 1}^{B}$ drawn i.i.d. from $\mathcal{D}_S$, then outputs a classifier $\hat y_t : \{\pm 1\}^n \rightarrow \{\pm 1\}$. If $\mathop{\mathbb{E}}_{(x,y) \sim \mathcal{D}_S}[ \mathbbm{1}[\hat y_{t}(x) \neq y] ] \leq \epsilon$ (i.e. $\hat y_{t}$ agrees with the correct parity on at least a $(1-\epsilon)$ fraction of inputs), the learner is said to have solved the parity problem with error $\epsilon$ in $t$ iterations ($tB$ samples); the smallest $t$ for which this is true is the \emph{convergence time} $t_c$. A learner may also output an initial classifier $\hat y_0$ before observing any data.

This formulation permits \emph{improper} function classes (i.e. other than parities over subsets $S'$) for the parity learning problem. In particular, we will focus on hypothesis classes of continuous functions $f : \{\pm 1\}^n \times \Theta \rightarrow \RR$, which map to classifiers $\hat y(x) = \mathop{\mathrm{sign}}(f(x;\theta))$\footnote{When $f(x;\theta) = 0$ in practice (e.g. with sign initialization), we break the tie arbitrarily. We ensure in the theoretical analysis that this does not happen.}. When $\Theta$ is a vector space over $\RR$, a ubiquitous online learning algorithm is gradient descent (GD). For a choice of loss function $\ell : \{\pm 1\} \times \RR \rightarrow \RR$, initialization $\theta_0$,
learning rate schedule $\{\eta_t\}_{t=1}^T \subseteq \RR$ and weight-decay schedule $\{\lambda_t\}_{t=1}^T \subseteq \RR$, GD defines iterative updates
\begin{equation}
\theta_{t+1} \leftarrow (1-\lambda_t)\theta_t - \eta_t \nabla_\theta \left( \frac{1}{B} \sum_{i=1}^B \ell( y_{t,i}, f(x_{t,i}; \theta_t) ) \right),
\label{eq:gd}
\end{equation}
where $f$ (the \emph{architecture}) and $\ell$ are assumed to be such that this gradient (more generally, subgradient) is well-defined.
In this context, \emph{online} and \emph{stochastic} gradient descent (OGD/SGD) are equivalent names for the update rule (\ref{eq:gd}).

A fundamental object of study in deep learning is the multi-layer perceptron (MLP). In this setting, a 2-layer MLP with \emph{width} $h$ and activation function $\sigma : \RR \rightarrow \RR$, parameterized by $W \in \RR^{r \times n}, b \in \RR^r, u \in \RR^{r}$, specifies the function
\[ f(x ; W, b, u) = u^\top \sigma(Wx+b), \]
where $\sigma(\cdot)$ is applied entrywise. It is standard to use GD to jointly update the network's parameters.
Our results include positive results about ``single neurons'': MLPs with width $r=1$. We note that for our theoretical analysis, when training MLPs with GD, we allow for different learning rate and weight decay schedule for the different layers.

Finally, we will analyze randomized learning algorithms, such as GD with random initialization $\theta_0$, whose iterates $\theta_t$ (and thus classifiers $\hat y$) are random variables even when the samples are not. A learning algorithm has \emph{permutation symmetry} if, for all sequences of data $\{(x_{t,i},y_{t,i})\}$, the classifiers $\hat y_t \circ \pi$ resulting from feeding $\{(\pi(x_{t,i}),y_{t,i})\}$ to the learner have identical distributions as $\pi$ ranges over all permutations of indices. The neural architectures and initializations (and thus, SGD) considered in this work are permutation-symmetric; for this reason, it is convenient for notation to choose $S = [k]$ as the canonical $(n,k)$-parity learning problem, without loss of generality.

\subsection{Additional related work}\label{app:relwork}

\paragraph{Feature learning using GD on neural networks.} A line of recent work has focused on understanding the feature learning ability of gradient descent dynamics on neural networks. These analyses go beyond the Neural Tangent Kernel (NTK) regime, where they show a separation between learning with fixed features versus GD on neural networks, for these problems. Several of these works assume structure (often "sparse") in the input data which is useful for the prediction task, and helps avoid computational hardness. In contrast, our work focuses on studying hard problems at their computational limit. Here we discuss the most relevant works in detail:

A line of work (\cite{diakonikolas2020approximation,yehudai2020learning, frei2020agnostic}) focuses on learning a single non-linearity $y = \sigma(w^\top x)$ (typically $\sigma(\cdot)$ is the ReLU or sigmoid) using gradient-based methods. These works obtain polynomial-time convergence guarantees when the distribution satisfies a \emph{spread} condition. These results do not extend to the Boolean hypercube.

\cite{daniely2020learning} also study the problem of learning sparse parities using neural networks. One key difference from our work is that they assume a modified version of the problem, where the input distribution is not uniform over the hypercube, but instead leaks information about the label. In particular, the distribution ensures that the relevant parity bits always have the same value. \cite{shi2021theoretical} generalize this setting by considering a setting where labels are generated based on certain class specific patterns and the data itself is generated using these patterns with some extra background patterns. This also embeds information in the data itself regarding the label, unlike our setting, where the labels are uncorrelated with the input features. Under these structural assumption, the papers study how GD on a two-layer network can learn useful features in polynomial time. Both these analysis also exploits the first gradient step to find useful features. \cite{shi2021theoretical} additionally require a second step to refine the features.

\cite{ba2022high} show how the first gradient step is important for feature leaning. In particular, they show that first update is essentially rank-1 and aligns with the linear component of the underlying function. The functions we consider (parity) do not have a linear component.

\cite{abbe2022initial} define a notion of initial alignment between the network at initialization and the target function and show that it is essential to get polynomial time learnability with noisy gradients on a fully connected network. Our MLP results also exploit the correlation between the gradient and the label to ensure that the gradient update gives us signal.

\cite{frei2022random} also study learnability of a parity-like function with $k=2$ under noisy labels. The paper analyzes early stopping GD for learning the underlying labeling function. Our setup is quite different from theirs and can handle $k>2$.

In concurrent work, \cite{damian2022neural} consider the problem of learning polynomials which depend on few relevant directions using gradient descent on a two-layer network. They assume that the distributional Hessian of the ground truth function spans exactly the subspace of the relevant direction. Using this, they show that gradient descent can learn the relevant subspace with sample complexity scaling as $d^2$ and not $d^p$ where $p$ is the degree of the underlying polynomial as long as the number of relevant directions is much less than $d$. Their proof technique is similar to our two-layer MLP result which also exploits correlation in the first gradient step. However, for our setting, the distributional Hessian has rank 0 and does not satisfy their assumptions.