# To grok or not to grok: Disentangling generalization and memorization on
corrupted algorithmic datasets

Darshil Doshi
&Aritra Das33footnotemark: 3
&Tianyu He22footnotemark: 2 33footnotemark: 3
&Andrey Gromov 22footnotemark: 2 33footnotemark: 3
Corresponding author<EMAIL_ADDRESS>Matter Theory Center, University
of Maryland, College ParkDepartment of Physics, University of Maryland,
College ParkMeta AI

###### Abstract

Robust generalization is a major challenge in deep learning, particularly when
the number of trainable parameters is very large. In general, it is very
difficult to know if the network has memorized a particular set of examples or
understood the underlying rule (or both). Motivated by this challenge, we
study an interpretable model where generalizing representations are understood
analytically, and are easily distinguishable from the memorizing ones. Namely,
we consider two-layer neural networks trained on modular arithmetic tasks
where ($\xi\cdot 100\%$) of labels are corrupted (_i.e._ some results of the
modular operations in the training set are incorrect). We show that (i) it is
possible for the network to memorize the corrupted labels _and_ achieve
$100\%$ generalization at the same time; (ii) the memorizing neurons can be
identified and pruned, lowering the accuracy on corrupted data and improving
the accuracy on uncorrupted data; (iii) regularization methods such as weight
decay, dropout and BatchNorm force the network to ignore the corrupted data
during optimization, and achieve $100\%$ accuracy on the uncorrupted dataset;
and (iv) the effect of these regularization methods is (“mechanistically”)
interpretable: weight decay and dropout force all the neurons to learn
generalizing representations, while BatchNorm de-amplifies the output of
memorizing neurons and amplifies the output of the generalizing ones. Finally,
we show that in the presence of regularization, the training dynamics involves
two consecutive stages: first, the network undergoes the _grokking_ dynamics
reaching high train _and_ test accuracy; second, it unlearns the memorizing
representations, where train accuracy suddenly jumps from $100\%$ to
$100(1-\xi)\%$. 111Code to reproduce our results is available at
https://github.com/d-doshi/Grokking.git

##  1 Introduction

The astounding progress of deep learning in the last decade has been
facilitated by massive, high-quality datasets. Annotated real-world datasets
inevitably contain noisy labels, due to biases of annotation schemes (Paolacci
et al., 2010; Cothey, 2004) or inherent ambiguity (Beyer et al., 2020). A key
challenge in training large models is to prevent overfitting the noisy data
and attain robust generalization performance. On the other hand, in large
models, it is possible for memorization and generalization to coexist (Zhang
et al., 2017; 2021). By and large, the tussle between memorization and
generalization, especially in the presence of label corruption, remains poorly
understood.

In generative language models the problem of memorization is even more
nuanced. On the one hand, some factual knowledge is critical for the language
models to produce accurate information. On the other hand, verbatim
memorization of the training data is generally unwanted due to privacy
concerns (Carlini et al., 2021). The full scope of the conditions that affect
memorization in large language models (LLMs) is still under investigation
(Carlini et al., 2022; Tirumala et al., 2022). In particular, larger models as
well as repeated data (Hernandez et al., 2022) favor memorization.

In this work, we take an approach to disentangling generalization and
memorization motivated by theoretical physics. Namely, we focus on simple,
analytically tractable, (“mechanistically”) interpretable models trained on
algorithmic datasets. More concretely, we study two-layer Multilayer
Perceptrons (MLP) trained on modular arithmetic tasks. In this case, the
network is trained to learn a rule, such as $z=(m+n)\%p\equiv
m+n\,\,\textrm{mod}\,\,p$ from a number of examples. These tasks are phrased
as classification problems. We then force the network to memorize $\xi\cdot
100\%$ of the examples by corrupting the labels, _i.e._ injecting data points
where $(m+n)\%p$ is mapped to $z^{\prime}\neq z$. In this setting, the
generalizing representations are understood qualitatively Nanda et al. (2023);
Liu et al. (2022) as well as quantitatively, thanks to the availability of an
analytic solution for the representations described in Gromov (2023). We
empirically show that sufficiently large networks can memorize the corrupted
examples while achieving nearly perfect generalization on the test set. The
neurons responsible for memorization can be identified using inverse
participation ratio (IPR) and pruned away, leading to perfect accuracy on the
un-corrupted dataset. We further show that regularization dramatically affects
the network’s ability to memorize the corrupted examples. Weight decay,
Dropout and BatchNorm all prevent memorization, but in different ways: Weight
decay and Dropout eliminate the memorizing neurons by converting them into
generalizing ones, while BatchNorm de-amplifies the signal coming from the
memorizing neurons without eliminating them.

###  1.1 Preliminaries

Figure 1: Training curves in various phases. All plots are made for networks
trained with data-fraction $\alpha=0.5$ and corruption-fraction $\xi=0.35$.
(a) No regularization: _Coexistence_ of generalization and memorization – both
train and test accuracies are high. (b)(c) Adding weight decay: The network
generalizes on the test data but does not memorize the corrupted training
data, resulting in a negative generalization gap! Remarkably, the network
predicts the “true” labels for the corrupted examples. We term these phases
_Partial Inversion_ and _Full Inversion_ , based on the degree of memorization
of corrupted data. (“inversion” refers to test accuracy being higher than
train accuracy.)

##### Grokking modular arithmetic.

We consider the modular addition task $(m+n)\,\%\,p$. It can be learned by a
two-layer MLP. Explicitly, the network function takes form

$\displaystyle{\bm{f}}(m,n)={\bm{W}}\phi\left({\bm{W}}_{in}({\bm{e}}_{m}\oplus{\bm{e}}_{n})\right)={\bm{W}}\phi\left({\bm{U}}\,{\bm{e}}_{m}+{\bm{V}}\,{\bm{e}}_{n}\right)\,.$
(1)

Here, ${\bm{e}}_{m},{\bm{e}}_{n}\in\mathbb{R}^{p}$ are one_hot encoded numbers
$m,n$. “$\oplus$” denotes concatenation of vectors
$({\bm{e}}_{m}\oplus{\bm{e}}_{n}\in\mathbb{R}^{2p})$.
${\bm{W}}_{in}\in\mathbb{R}^{N\times 2p}$ and ${\bm{W}}\in\mathbb{R}^{p\times
N}$ are the first and second layer weight matrices, respectively. $\phi$ is
the activation function. ${\bm{W}}_{in}$ is decomposed into two $N\times P$
blocks: ${\bm{U}},{\bm{V}}\in\mathbb{R}^{N\times p}$. ${\bm{U}},{\bm{V}}$
serve as embedding vectors for $m,n$ respectively.
${\bm{f}}(m,n)\in\mathbb{R}^{p}$ is the network-output on one example pair
$(m,n)$. The targets are one_hot encoded answers ${\bm{e}}_{(m+n)\,\%\,p}$.

The dataset consists of $p^{2}$ examples pairs; from which we randomly pick
$\alpha p^{2}$ examples for training222Throughout the text, networks are
trained with Full-batch AdamW and MSE loss, unless otherwise stated., and use
the other $(1-\alpha)p^{2}$ examples as a test set. We will refer to $\alpha$
as the _data fraction_.

This setup exhibits grokking: delayed and sudden occurrence of generalization,
long after memorization (Power et al., 2022; Nanda et al., 2023; Liu et al.,
2022; Gromov, 2023).

##### Periodic weights.

Choosing quadratic activation function $\phi(x)=x^{2}$ makes the problem
analytically solvable. Upon grokking, the trained network weights are
qualitatively similar to the analytical solution presented in Gromov
(2023).333Periodicity of weights is approximate in trained networks. The
following analytical expression for the network weights leads to $100\%$
generalization (for sufficiently large width)

$\displaystyle\begin{aligned}
&{U}_{ki}=\left(\frac{2}{N}\right)^{-\frac{1}{3}}\cos{\left(\frac{2\pi}{p}i\sigma(k)+\phi^{(u)}_{k}\right)}\,,&&{V}_{kj}=\left(\frac{2}{N}\right)^{-\frac{1}{3}}\cos{\left(\frac{2\pi}{p}j\sigma(k)+\phi^{(v)}_{k}\right)}\,,&\\\
&{W}_{qk}=\left(\frac{2}{N}\right)^{-\frac{1}{3}}\cos{\left(-\frac{2\pi}{p}q\sigma(k)-\phi^{(u)}_{k}-\phi^{(v)}_{k}\right)}\,,\end{aligned}$
(2)

where $\sigma(k)$ denotes a random permutation of $k$ in $S_{N}$ – reflecting
the permutation symmetry of the hidden neurons. The phase $\phi_{k}^{(u)}$ and
$\phi_{k}^{(v)}$ are uniformly i.i.d. sampled between $(-\pi,\pi]$. We refer
the reader to Appendix C for further details about the analytical solution.

##### Inverse participation ratio.

To characterize the generalizing representations quantitatively, we introduce
a quantity familiar from the physics of localization: the _inverse
participation ratio (IPR)_. It’s role is to detect periodicity in the weight
matrix.

Let ${\bm{U}}_{k\cdot},{\bm{V}}_{k\cdot}({\bm{W}}_{k\cdot})$ denote the
$k^{th}$ row(column) vectors of the weights ${\bm{U}},{\bm{V}}({\bm{W}})$.
Consider the discrete Fourier transforms of these vectors, denoted by
$\widetilde{\bm{U}}_{k\cdot},\widetilde{\bm{V}}_{k\cdot}(\widetilde{\bm{W}}_{k\cdot})$.
The Fourier decompositions of the periodic rows(columns) of these weights are
highly localized. We leverage this to quantify the similarity of trained
weights to the analytic solution (equation 2). To that end, we define the
_Inverse Participation Ratio (IPR)_ (Girvin & Yang, 2019; Pastor-Satorras &
Castellano, 2016) for these vectors: 444In general IPR is defined as $(\lVert
U_{k\cdot}\rVert_{2r}/\lVert U_{k\cdot}\rVert_{2})^{2r}$. We set $r=2$, a
common choice.

$\displaystyle\text{IPR}^{(u)}_{k}:=\left(\frac{\lVert\widetilde{\bm{U}}_{k\cdot}\rVert_{4}}{\lVert\widetilde{\bm{U}}_{k\cdot}\rVert_{2}}\right)^{4}\,;$
$\displaystyle\text{IPR}^{(v)}_{k}:=\left(\frac{\lVert\widetilde{\bm{V}}_{k\cdot}\rVert_{4}}{\lVert\widetilde{\bm{V}}_{k\cdot}\rVert_{2}}\right)^{4}\,;$
$\displaystyle\text{IPR}^{(w)}_{k}:=\left(\frac{\lVert\widetilde{\bm{W}}_{\cdot
k}\rVert_{4}}{\lVert\widetilde{\bm{W}}_{\cdot k}\rVert_{2}}\right)^{4}\,;$ (3)

where $\lVert\cdot\rVert_{P}$ denotes the $L^{P}$-norm of the vector. One can
readily see from equation 3 that $IPR\in[1/p,1]$, with higher values for
periodic vectors and lower values for non-periodic ones. It is useful to
quantify the periodicity of each neuron in the hidden layer of the network
(indexed by $k$). We define per-neuron IPR by averaging over the weight-
vectors connected to the neuron:

$\displaystyle\text{IPR}_{k}:=\frac{1}{3}\left(\text{IPR}^{(u)}_{k}+\text{IPR}^{(v)}_{k}+\text{IPR}^{(w)}_{k}\right)\,.$
(4)

Figure 2: Grokking the modular arithmetic task over $\mathbb{Z}_{97}$ with
2-layer MLP, trained with AdamW. (a) Sharp transition in test accuracy, long
after overfitting.
$\overline{\text{IPR}}\coloneqq\mathbb{E}_{k}\left[\text{IPR}_{k}\right]$
monotonically increases over time, indicating periodic representations. (b)
Example row vector $(U_{k\cdot})$ before and after grokking. Generalization is
achieved through periodic weights (equation 2). (c) Histogram of per-neuron
IPRs before and after grokking – the distribution shifts to high IPR.

In trained networks, a larger population of high-IPR neurons leads to better
generalization. To quantify the overall similarity of the trained network to
the analytical solution (equation 2) we average $\text{IPR}_{k}$ over all
hidden neurons:
$\overline{\text{IPR}}\coloneqq\mathbb{E}_{k}\left[\text{IPR}_{k}\right]$.

##### Label corruption.

To introduce label corruption, we choose a fraction $\xi$ of the training
examples, and replace the true labels with random labels. The corrupted labels
are generated from a uniform random distribution over all the labels. This is
often called _symmetric label noise_ in literature (Song et al., 2022). This
type of label corruption does not introduce label asymmetry.

###  1.2 Related Works

##### Grokking

Grokking was first reported by Power et al. (2022) for Transformers trained on
modular arithmetic datasets. Liu et al. (2022; 2023) presented explanations
for grokking in terms of quality of representation learning as well as weight-
norms. They also found grokking-like behaviours on other datasets. Nanda et
al. (2023); Zhong et al. (2023) reverse-engineered grokked Transformer models
and revealed the underlying learned algorithm. Barak et al. (2022) observed
grokking on the sparse-parity problem, which they attributed to the
amplification of the so-called Fourier gap. Merrill et al. (2023) examined
grokking in the sparse-parity setup and attributed it to the emergence of
sparse subnetworks. Žunkovič & Ilievski (2022) discuss solvable models for
grokking in a teacher-student setup. Recently, Notsawo et al. (2023)
investigated the loss landscape at early training time and used it to predict
the occurrence of grokking. Varma et al. (2023) demonstrated new grokking-
behaviours by investigating the critical amount of data needed for grokking.

##### Label Noise

Zhang et al. (2017) showed that neural networks can easily fit random labels.
A large body of work has been dedicated to developing techniques for robust
training with label noise: explicit regularization(Hinton et al., 2012; Ioffe
& Szegedy, 2015), noise-robust loss functions (Ghosh et al., 2017), sample
selection (Katharopoulos & Fleuret, 2018; Paul et al., 2021), noise modeling
and dedicated architectures (Xiao et al., 2015; Han et al., 2018; Yao et al.,
2019) etc. We refer the reader to recent reviews by Song et al. (2022); Liang
et al. (2022) for a comprehensive overview.

##### Regularization and Pruning

Morcos et al. (2018) investigated the interplay between label corruption and
regularization in vision tasks. They showed that BatchNorm and Dropout have
qualitatively different effects. Dropout effectively reduces the size of the
network but does not discourage over-reliance on a few special directions in
the weights space; while BatchNorm spreads the generalizing features over many
dimensions. Bartoldson et al. (2020) found that there is a trade-off between
stability and generalization when pruning trained networks.

##  2 Grokking with label noise

Grokking MLP on modular addition dataset is remarkably robust to label
corruption. Even without explicit regularization, the model can generalize to
near $100\%$ accuracy with sizable label corruption. In many cases, the
network surprisingly manages to “correct” some of the corrupted examples,
resulting in _Inversion_ (i.e. test accuracy $>$ training accuracy). We
emphasise that this is in stark contrast to the common belief that grokking
requires explicit regularization. Adding regularization makes grokking more
robust to label corruption, with stronger _Inversion_. To quantify these
behaviours, we observe the training and test performance of the network with
varying sizes of training dataset, amount of label corruption and
regularization. We summarize our findings in the form of empirical phase
diagrams (Figure 3), featuring the following phases. More detailed phase
diagrams along with training and test performances are presented in Appendix
G.

##### Coexistence.

In the absence of explicit regularization the network generalizes on test data
_and_ memorizes the corrupted training data. In other words, both train and
test accuracy are close to $100\%$ despite label corruption.

##### Partial Inversion.

In this phase, the network generalizes on the test data but only memorizes a
fraction of the corrupted training data. Remarkably, the network predicts the
“true” labels on the remaining corrupted examples. In other words, the network
_corrects_ a fraction of the corrupted data. Consequently, we see $<100\%$
train accuracy but near-$100\%$ test accuracy, resulting in a _negative_
generalization gap (Figure 1(b))! We term this phenomenon _Partial Inversion_
; where “inversion” refers to the test accuracy being higher than train
accuracy. Remarkably, partial inversion occurs even in the absence of any
explicit regularization, but only when there is ample training data (left-most
panels in Figure 3(a,b)).

##### Full Inversion.

Upon adding regularization, the network often generalizes on the test set but
does not memorize any of the corrupted examples. This results in training
accuracy plateauing close to $100(1-\xi)\%$ ($\xi$ is the _corruption
fraction_); while the test accuracy is at $100\%$ (Figure 1(c)). Increasing
regularization enhances this phase (middle columns in Figure 3(a,b)).

##### Memorization.

With very high label corruption and/or very low training data-fraction the
network memorizes all of the training data (corrupted as well as non-
corrupted), but does not generalize. We get close to $100\%$ accuracy on
training data, but random guessing accuracy on test data in this phase.

##### Forgetting.

With very high weight decay, the network performs poorly on both training and
test data (right-most panel in Figure 3(a)). However, upon examination of the
training curves (Figure 7(a)) we find that the accuracies plummet _after_ the
network has generalized! This performance drop is a result of the norms of
weights collapsing to arbitrarily small values (Figure 7(b)). Forgetting
exclusively occurs in grokked networks with activation functions of degree
higher than 1 (e.g. we use $\phi(x)=x^{2}$) and with adaptive optimizers. (See
Appendix D for a detailed discussion.)

##### Confusion.

With very high dropout the network becomes effectively narrow, resulting in
reduced capacity to memorize. Consequently, even in the absence of grokking,
the network does not manage to memorize all the corrupted labels (right-most
panel in Figure 3(b)).

(a) Regularization with weight decay. (wd denotes the weight decay value)

(b) Regularization with Dropout.

Figure 3: Modular Addition phase diagrams with various regularization methods. A larger data fraction leads to more “correct” examples, leading to higher corruption-robustness. Increasing regularization, in the form of weight decay or dropout, enhances robustness to label corruption and facilitates better generalization. | Coexistence |  Partial Inversion |  Full Inversion | Memorization |  Forgetting / Confusion
---|---|---|---|---|---
Memorizes corrupted train data | Yes | Partially | No | Yes | No
Generalizes on test data | Yes | Yes | Yes | No | No
Table 1: Qualitative characterization of phases.

###  2.1 Generalizing and Memorizing Sub-networks

(a) Regularization with weight decay. (wd denotes the weight decay value)

(b) Regularization with Dropout. (dp denotes the Dropout probability)

Figure 4: Distribution of per-neuron IPR for trained networks in various
phases. Coexistence phase has a bimodal IPR distribution, where the high and
low IPR neurons facilitate generalization and memorization, respectively.
Regularization with weight decay or Dropout shifts the IPR distribution
towards higher values – Generalizing neurons get more populous compared to
memorizing ones; resulting in more robust generalization and _Inversion_. All
plots are made for networks trained with data-fraction $\alpha=0.5$ and
corruption-fraction $\xi=0.35$; with various regularization strengths.

The generalizing features in our setup are “mechanistically” interpretable
(equation 2). As a result, we can leverage the measures of periodicity such as
IPR to quantitatively characterize various phases presented above.
Furthermore, it also allows us to analyze the effect of various regularization
schemes. To that end, we state the following hypothesis and extensively
examine it with experiments.

######  Hypothesis 2.1.

For a two-layer MLP with quadratic activation trained on (one-hot encoded)
modular addition datasets with label corruption, the high IPR neurons
facilitate generalization via feature-learning whereas the low IPR neurons
cause memorization of corrupted training data.

To test Hypothesis 2.1 we study the distribution of per-neuron IPR in various
phases, with different levels of regularization. In the _Coexistence_ phase,
the network simultaneously generalizes and memorizes the corrupted train data.
From Hypothesis 2.1, we expect the presence of both generalizing and
memorizing neurons. Indeed, we find a bi-modal distribution of IPRs in trained
networks (left column of Figure 4). Remarkably, the network not only
distinguishes between the tasks of feature-learning and memorization, it also
assigns distinct and independent sub-networks to each.

As mentioned in the previous section, adding regularization in the form of
weight decay or Dropout makes the generalization more robust and facilitates
_Inversion_. Since these phases see decreased memorization of corrupted
labels, we expect to see a decrease in the number of memorizing neurons and a
corresponding increase in the generalizing ones. Indeed, we observe a
distribution-shift towards higher IPR in the _Partial Inversion_ phase (middle
column of Figure 4). In the _Full Inversion_ phase (induced by weight decay or
dropout), we see that memorizing neurons are almost entirely eliminated in
favour of generalizing ones (right column of Figure 4). This provides a
mechanism by which the network recovers the “true” labels even for the
corrupted data – in the absence of memorizing sub-network, its predictions are
guided by the underlying rule and not corrupted samples.

BatchNorm also has a regularizing effect on training (Figure 8), albeit,
through a qualitatively distinct mechanism. In this case, the low IPR neurons
persist even in the _Full Inversion_ phase (Figure 5). However, the BatchNorm
parameters de-amplify the output of the memorizing neurons (Figure 5), biasing
the network towards generalizing representations. BatchNorm effectively prunes
the network and it learns to do this directly from the data! Note that this
effect of BatchNorm acts in conjunction with the implicit regularization from
the mini-batch noise in this case. We find that decreasing the batch-size
enhances _Inversion_ even without explicit regularization (Figure 9).

Figure 5: IPR distributions of networks with/without BatchNorm; the
correlation between BatchNorm weights and IPR. The IPR distributions in the
two models do not exhibit significant differences. However, BatchNorm helps
generalization by assigning higher weights ($\gamma_{k}$) to high IPR neurons
compared to low IPR neurons. Both models are trained at data fraction
$\alpha=0.65$ and corruption fraction $\xi=0.2$, with batch-size = 64.

To further isolate the effect of individual neurons on performance, we perform
the two complementary pruning experiments (Figure 6) : Starting from a trained
network in the Coexistence phase, we gradually prune out neurons, one-at-a-
time, (i) starting from the lowest IPR neuron (ii) starting from the highest
IPR neuron, while keeping track of training and test accuracies. To
distinguish network’s ability to memorize, we also track the accuracies on
corrupted and uncorrupted parts of the training dataset. Before pruning, the
network has near-$100\%$ accuracy on both train and test data. In case (i)
(Figure 6(left)), we see a monotonic decrease in accuracy on corrupted
examples as more and more low IPR neurons are pruned. The accuracy on test
data as well as ucorrupted training data remains high. The resulting (overall)
training accuracy plateaus around $100(1-\xi)\%$, corresponding to the
fraction of uncorrupted examples. In case (ii) (Figure 6(right)), as more and
more high IPR neurons are pruned, we see a decline in test accuracy (as well
as the accuracy on the uncorrupted training data). The network retains its
ability to memorize, so the accuracy on corrupted training data remains high.

We emphasize that, although each generalizing neuron learns the useful
periodic features, the overall computation in the network is emergent: it is
performed collectively by _all_ generalizing neurons simultaneously.

Figure 6: Progressively pruning neurons from a trained model (_Coexistence_
phase), based on IPR. (data fraction $\alpha=0.5$, corruption fraction
$\xi=0.35$). (Left) Pruning out neurons starting from the lowest IPR. The
accuracy on corrupted training data decreases, providing evidence that low IPR
neurons are responsible for memorization. Test accuracy as well as accuracy on
uncorrupted training data remains high. (Right) Pruning out neurons starting
from the highest IPR. Test accuracy decreases providing evidence that high IPR
neurons are responsible for generalization. Accuracy on corrupted training
data remains high.

##  3 Effect of Regularization

In the previous Section, we demonstrated that the network can identify the
correct features and construct a generalizing sub-network even in the presence
of label corruption. We also showed that explicit regularization facilitates
robust generalization by strengthening the generalizing representations. In
this section, we offer a further interpretation of how regularization
techniques such as weight decay, Dropout and BatchNorm aid in building
generalizing sub-networks.

###  3.1 Weight decay

Weight decay has been recognized as an effective technique to enhance
generalization. It is widely accepted that weight decay helps limit a
network’s _capacity_ , thereby preventing overfitting (Neyshabur et al., 2015;
Bartlett, 1998). Despite this intuition, the exact effect of weight decay on
general architectures and tasks is far from well-understood. Moreover, the
definition of capacity also remains an open question (Neyshabur et al., 2017).
Liu et al. (2023) discussed the effect of weight-norms and weight decay on
grokking tasks, showing a distinction in the weight norms of generalizing and
memorizing solutions. They showed that weight decay naturally favours
generalization by decreasing the weight norms.

The difference between generalizing and memorizing solutions is made sharper
in our setting by means of label corruption. As shown in Figure 3(a),
increasing the weight decay has two notable effects: (i) Generalization with
less data: The phase boundary between the _Memorization_ and _Coexistence_
phases shifts toward lower data-fraction. Thus, weight decay encourages
generalization even when training data is scarce; in line with classical ideas
(Bartlett, 1998). (ii) Emergence of the _Full Inversion_ phase: In this phase,
the network “corrects” the labels for corrupted data points. This observation
corroborates the idea that weight decay prevents memorization.

The effect of weight decay on representations is elucidated by the experiment
presented in Figure 4(a). For a fixed amount of data and label corruption,
increased weight decay results in a diminished proportion of low IPR neurons
which constitute the memorizing sub-network; which brings the network closer
to the generalizing solution (Equation 2).

###  3.2 Dropout

Dropout, as introduced by Hinton et al. (2012), is a technique designed to
prevent overfitting by randomly deactivating a fraction of neurons during each
optimization step. The phase diagrams (Figure 3(b)) indeed support the
regularizing effect of Dropout. We see enhanced generalization and the onset
of _Full Inversion_ as the dropout probabilities are increased.

It is believed that dropping neurons during training prevents complex co-
adaptations between neurons in different layers, encouraging each neuron to
learn useful features (Hinton et al., 2012; Srivastava et al., 2014). This
intuition becomes more crisp in our setup. Note that the generalizing solution
constitutes of periodic neurons, where each neuron corresponds to weight-
vectors of a given “frequency” $\sigma(k)$. For such a solution to work, we
need $p$ different frequencies (see Equation C). In a wide network ($N=500$ in
our setup), there is room for redundancy.555Although, the analytical solution
becomes exact only at large $N$. Random phase approximation gives sub-leading
corrections of order $1/\sqrt{N}$. Due to this redundancy, generalizing
solutions can be more robust to dropped neurons. As a result, training with
Dropout gets attracted to such periodic solutions. Memorizing solutions, on
the other hand, do not seem to have such a feature. Another common heuristic
is that dropout makes the network effectively narrower, resulting in a
decreased capacity to memorize. This intuition is in line with results about
the similarity between Dropout and weight decay Wager et al. (2013).666Wager
et al. (2013) showed the equivalence for _linear_ models. In our setup, one
can find empirical corroboration to this idea by observing the similarity
between the phase diagrams for Dropout and weight decay (Figure 3 (a) and
(b)).

The picture is supplemented by studying the IPR distributions (Figure 4(b)).
It shifts towards high IPR values upon increasing Dropout, explaining the
regularizing effect of Dropout. Notably, the similarity between Dropout and
weight decay also carries over to the IPR distributions.

Extremely high dropout probabilities behave differently from Extremely high
weight decay. Unlike the _Forgetting_ phase, where elevated weight decay
values disrupts the optimization process and destroys grokking; we observe
that high dropout probabilities do not eliminate grokking. They do, however,
limit the network’s capacity to memorize. As a result, we find _Confusion_
phase at high dropout probabilities and high label noise. (right-most column
in Figure 3)

###  3.3 BatchNorm

BatchNorm (Ioffe & Szegedy, 2015) was introduced to reduce “internal covariate
shift”, and is a form of regularization. It acts component-wise on $N$-dim
input $\mathbf{z}=(z_{1},\cdots,z_{N})$

$\displaystyle\hat{z}_{k}=\gamma_{k}\frac{z_{k}-\mathbb{E}_{B}[z_{k}]}{\sqrt{\mathrm{Var}_{B}[z_{k}]}}+b_{k}\,,$
(5)

where the mean and variance are taken over the mini-batches. $\gamma_{k}$ and
$b_{k}$ are trainable weights and bias of the BatchNorm layer. A subsequent
study by Morcos et al. (2018) showed that BatchNorm has greater resilience to
label noise compared to Dropout on vision tasks. They argued that the effects
of Dropout and BatchNorm are qualitatively different. We clarify this in our
setting.

To that end, we repeat our experiments for networks with BatchNorm applied
after activation functions (post-BN) on the corrupted modular addition data
(no other regularization).

$\displaystyle{\bm{f}}(m,n)={\bm{W}}\texttt{BN}\left(\phi\left({\bm{U}}\,{\bm{e}}_{m}+{\bm{V}}\,{\bm{e}}_{n}\right)\right)\,.$
(6)

In Figure 5, we find that the BatchNorm layer does not affect the IPR
distribution of fully connected layers significantly. Instead, BatchNorm
adjusts its own weights to de-amplify the low IPR (memorizing) neurons, which
leads to a strong correlation between $\text{IPR}_{k}$ and the corresponding
$\gamma_{k}$ of the BatchNorm layer. This correlation manifests as two
distinct clusters in the scatter-plot (Figure 5 right), with the cluster with
high IPR values being shifted towards larger weights and vice versa. Bias
$b_{k}$ in this setting remains almost unchanged, and the clustering behavior
is not observed between bias $b_{k}$ and $\text{IPR}_{k}$.

##  4 Conclusion

We have investigated the interplay between memorization and generalization in
two-layer neural networks trained on algorithmic datasets. The memorization
was induced by creating a dataset with corrupted labels, forcing the network
to memorize incorrect examples and learn the rule at the same time. We have
found that the memorizing neurons can be explicitly identified and pruned away
leading to perfect generalization.

Next we showed that various regularization methods such as weight decay,
Dropout and BatchNorm encourage the network to ignore the corrupted labels and
lead to perfect generalization. We have leveraged the interpretability of the
representations learnt on the modular addition task to quantitatively
characterize the effects of these regularization methods. Namely, we have
found that the weight decay and dropout act in a similar way by reducing the
number of memorizing neurons; while BatchNorm de-amplifies the output of
regularizing neurons without eliminating them, effectively pruning the
network.

##  5 Limitations and Future work

Although the trained models in our setup are completely interpretable by means
of analytical solutions, the training dynamics that lead to these solutions
remains an open question. Solving the dynamics should also shed more light on
the effects of regularization and label corruption.

While our work is limited to the exactly solvable models and algorithmic
datasets, it paves the way to quantitative investigation of generalization and
memorization in more realistic settings.

#### Acknowledgments

A.G.’s work at the University of Maryland was supported in part by NSF CAREER
Award DMR-2045181, Sloan Foundation and the Laboratory for Physical Sciences
through the Condensed Matter Theory Center.

## References

  * Barak et al. (2022) Boaz Barak, Benjamin L. Edelman, Surbhi Goel, Sham M. Kakade, eran malach, and Cyril Zhang.  Hidden progress in deep learning: SGD learns parities near the computational limit.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=8XWP2ewX-im.
  * Bartlett (1998) P.L. Bartlett.  The sample complexity of pattern classification with neural networks: the size of the weights is more important than the size of the network.  _IEEE Transactions on Information Theory_ , 44(2):525–536, 1998.  doi: 10.1109/18.661502.
  * Bartoldson et al. (2020) Brian R. Bartoldson, Ari S. Morcos, Adrian Barbu, and Gordon Erlebacher.  The generalization-stability tradeoff in neural network pruning, 2020\.
  * Beyer et al. (2020) Lucas Beyer, Olivier J. Hénaff, Alexander Kolesnikov, Xiaohua Zhai, and Aäron van den Oord.  Are we done with imagenet?, 2020.
  * Carlini et al. (2021) Nicholas Carlini, Florian Tramer, Eric Wallace, Matthew Jagielski, Ariel Herbert-Voss, Katherine Lee, Adam Roberts, Tom Brown, Dawn Song, Ulfar Erlingsson, et al.  Extracting training data from large language models.  In _30th USENIX Security Symposium (USENIX Security 21)_ , pp. 2633–2650, 2021.
  * Carlini et al. (2022) Nicholas Carlini, Daphne Ippolito, Matthew Jagielski, Katherine Lee, Florian Tramer, and Chiyuan Zhang.  Quantifying memorization across neural language models.  _arXiv preprint arXiv:2202.07646_ , 2022.
  * Cothey (2004) Viv Cothey.  Web-crawling reliability.  _Journal of the American Society for Information Science and Technology_ , 55(14):1228–1238, 2004.  doi: https://doi.org/10.1002/asi.20078.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/asi.20078.
  * Ghosh et al. (2017) Aritra Ghosh, Himanshu Kumar, and P Shanti Sastry.  Robust loss functions under label noise for deep neural networks.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 31, 2017.
  * Girvin & Yang (2019) Steven M. Girvin and Kun Yang.  _Modern Condensed Matter Physics_.  Cambridge University Press, 2019.  doi: 10.1017/9781316480649.
  * Gromov (2023) Andrey Gromov.  Grokking modular arithmetic, 2023.  URL https://arxiv.org/abs/2301.02679.
  * Han et al. (2018) Bo Han, Jiangchao Yao, Gang Niu, Mingyuan Zhou, Ivor Tsang, Ya Zhang, and Masashi Sugiyama.  Masking: A new perspective of noisy supervision.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 31, 2018.
  * Hernandez et al. (2022) Danny Hernandez, Tom Brown, Tom Conerly, Nova DasSarma, Dawn Drain, Sheer El-Showk, Nelson Elhage, Zac Hatfield-Dodds, Tom Henighan, Tristan Hume, et al.  Scaling laws and interpretability of learning from repeated data.  _arXiv preprint arXiv:2205.10487_ , 2022.
  * Hinton et al. (2012) Geoffrey E. Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan R. Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors, 2012.
  * Ioffe & Szegedy (2015) Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift, 2015.
  * Katharopoulos & Fleuret (2018) Angelos Katharopoulos and Francois Fleuret.  Not all samples are created equal: Deep learning with importance sampling.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 2525–2534. PMLR, 10–15 Jul 2018.
  * Liang et al. (2022) Xuefeng Liang, Xingyu Liu, and Longshan Yao.  Review–a survey of learning from noisy labels.  _ECS Sensors Plus_ , 1(2):021401, jun 2022.  doi: 10.1149/2754-2726/ac75f5.  URL https://dx.doi.org/10.1149/2754-2726/ac75f5.
  * Liu et al. (2022) Ziming Liu, Ouail Kitouni, Niklas S Nolte, Eric Michaud, Max Tegmark, and Mike Williams.  Towards understanding grokking: An effective theory of representation learning.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), _Advances in Neural Information Processing Systems_ , volume 35, pp. 34651–34663. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/dfc310e81992d2e4cedc09ac47eff13e-Paper-Conference.pdf.
  * Liu et al. (2023) Ziming Liu, Eric J Michaud, and Max Tegmark.  Omnigrok: Grokking beyond algorithmic data.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=zDiHoIWa0q1.
  * Merrill et al. (2023) William Merrill, Nikolaos Tsilivis, and Aman Shukla.  A tale of two circuits: Grokking as competition of sparse and dense subnetworks, 2023.
  * Morcos et al. (2018) Ari S. Morcos, David G. T. Barrett, Neil C. Rabinowitz, and Matthew Botvinick.  On the importance of single directions for generalization, 2018.
  * Nanda et al. (2023) Neel Nanda, Lawrence Chan, Tom Lieberum, Jess Smith, and Jacob Steinhardt.  Progress measures for grokking via mechanistic interpretability, 2023\.
  * Neyshabur et al. (2015) Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  Norm-based capacity control in neural networks.  In Peter Grünwald, Elad Hazan, and Satyen Kale (eds.), _Proceedings of The 28th Conference on Learning Theory_ , volume 40 of _Proceedings of Machine Learning Research_ , pp. 1376–1401, Paris, France, 03–06 Jul 2015. PMLR.  URL https://proceedings.mlr.press/v40/Neyshabur15.html.
  * Neyshabur et al. (2017) Behnam Neyshabur, Srinadh Bhojanapalli, David Mcallester, and Nati Srebro.  Exploring generalization in deep learning.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper_files/paper/2017/file/10ce03a1ed01077e3e289f3e53c72813-Paper.pdf.
  * Notsawo et al. (2023) Pascal Jr. Tikeng Notsawo, Hattie Zhou, Mohammad Pezeshki, Irina Rish, and Guillaume Dumas.  Predicting grokking long before it happens: A look into the loss landscape of models which grok, 2023.
  * Paolacci et al. (2010) Gabriele Paolacci, Jesse Chandler, and Panagiotis G. Ipeirotis.  Running experiments on amazon mechanical turk.  _Judgment and Decision Making_ , 5(5):411–419, 2010.  doi: 10.1017/S1930297500002205.
  * Pastor-Satorras & Castellano (2016) Romualdo Pastor-Satorras and Claudio Castellano.  Distinct types of eigenvector localization in networks.  _Scientific Reports_ , 6(1):18847, 2016.  doi: 10.1038/srep18847.  URL https://doi.org/10.1038/srep18847.
  * Paul et al. (2021) Mansheej Paul, Surya Ganguli, and Gintare Karolina Dziugaite.  Deep learning on a data diet: Finding important examples early in training.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 20596–20607. Curran Associates, Inc., 2021.
  * Power et al. (2022) Alethea Power, Yuri Burda, Harri Edwards, Igor Babuschkin, and Vedant Misra.  Grokking: Generalization beyond overfitting on small algorithmic datasets, 2022.
  * Song et al. (2022) Hwanjun Song, Minseok Kim, Dongmin Park, Yooju Shin, and Jae-Gil Lee.  Learning from noisy labels with deep neural networks: A survey.  _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–19, 2022.  doi: 10.1109/TNNLS.2022.3152527.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(56):1929–1958, 2014.  URL http://jmlr.org/papers/v15/srivastava14a.html.
  * Tirumala et al. (2022) Kushal Tirumala, Aram Markosyan, Luke Zettlemoyer, and Armen Aghajanyan.  Memorization without overfitting: Analyzing the training dynamics of large language models.  _Advances in Neural Information Processing Systems_ , 35:38274–38290, 2022.
  * Varma et al. (2023) Vikrant Varma, Rohin Shah, Zachary Kenton, János Kramár, and Ramana Kumar.  Explaining grokking through circuit efficiency, 2023.
  * Wager et al. (2013) Stefan Wager, Sida Wang, and Percy S Liang.  Dropout training as adaptive regularization.  _Advances in neural information processing systems_ , 26, 2013.
  * Xiao et al. (2015) Tong Xiao, Tian Xia, Yi Yang, Chang Huang, and Xiaogang Wang.  Learning from massive noisy labeled data for image classification.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2015.
  * Yao et al. (2019) Jiangchao Yao, Jiajie Wang, Ivor W. Tsang, Ya Zhang, Jun Sun, Chengqi Zhang, and Rui Zhang.  Deep learning from noisy image labels with quality embedding.  _IEEE Transactions on Image Processing_ , 28(4):1909–1922, 2019.  doi: 10.1109/TIP.2018.2877939.
  * Zhang et al. (2017) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  In _International Conference on Learning Representations_ , 2017.  URL https://openreview.net/forum?id=Sy8gdB9xx.
  * Zhang et al. (2021) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning (still) requires rethinking generalization.  _Communications of the ACM_ , 64(3):107–115, 2021\.
  * Zhong et al. (2023) Ziqian Zhong, Ziming Liu, Max Tegmark, and Jacob Andreas.  The clock and the pizza: Two stories in mechanistic explanation of neural networks, 2023.
  * Žunkovič & Ilievski (2022) Bojan Žunkovič and Enej Ilievski.  Grokking phase transitions in learning local rules with gradient descent, 2022.

##  Appendix A Further Experimental Details

In every plot presented, unless explicitly stated otherwise, we use the
following default setting:

##### Default setting

A 2-layer fully connected network with quadratic activation function and width
$N=500$, trained on $p=97$ Modular Addition dataset. The networks are trained
for 2000 optimization steps, using full-batch, MSE loss, AdamW optimizer,
learning rate $\eta=0.01$ and $\beta=(0.9,0.98)$.

###  A.1 Figures in Main Text

Figure 1: Default setting with $\alpha=0.5,\xi=0.35$ and (a) $wd=0$ (b) $wd=5$
and (c) $wd=15$.

Figure 2: Default setting with $\alpha=0.5,\xi=0.0,wd=5.0$.

Figure 3: Each plot is scanned over $17$ different data fractions $\alpha$
ranging from $0.1$ to $0.9$ and $19$ noise levels $\xi$ from $0.0$ to $0.9$
using the default setting.

Figures 4 and 6: Default setting with $\alpha=0.5,\xi=0.35$.

Figure 5: 2-layer fully connected network with quadratic activation function
and width $N=500$, with and without BatchNorm. We train for $\sim 8000$ steps,
with MSE loss, Adam optimizer with minibatch of size 64, and learning rate
$\eta=0.005$. Data fraction $\alpha=0.65$, corruption fraction $\xi=0.2$.

General setting for minibatch training: We train the network with batch size
$256$ for $2000$ steps. For batch sizes $8$ and $64$, we scale the number of
total steps so as to keep the number of epochs roughly equal. (In other words,
all the networks see the same amount of data.) The learning rate is scaled as
$0.01\cdot\sqrt{\text{batch size}/256}$.

###  A.2 Figures in Appendix

Figure 7: Default setting with $\alpha=0.5,\xi=0.35$ and $wd=20$

Figures 8, 9 and 10: Same as Figure 5. Except for the full batch training
receipt, which is the same as the default setting.

Figure 11: Default setting with different training steps shown in each phase.

Figures 12, 13 and 14: All networks are trained with batch-size = 64,
corruption fraction $\xi=0.3$ and various data fractions.

Figures 15 and 16: Same as Figure 3.

Figure 17: Used CrossEntropy loss, other settings remain the same as the
default.

Figures 18 and 19: Replace quadratic activation by ReLU, other settings remain
default.

Figures 20 and 21: Default setting, except width $N=5000$ and the networks are
trained for 10000 epochs.

Figure 22: Same setting as Figure 3. Trained on modular multiplication data
$z=x\cdot y\text{ mod }p$, with $p=97$.

##  Appendix B Details of resources used

Figure 3 requires 72 hours on a single NVIDIA RTX 3090 GPU.

All the remaining figures in the main text are either further experiments
performed on the data from the Figure 3 experiment; or are very quick to
perform on any machine.

We used $~{}300$ GPU hours for early-stage exploration and making phase
diagrams in the appendix.

###  B.1 Quantitative Characterization of phases

All phase diagrams are plotted using the empirical accuracies of trained
networks.

Trained networks with $\geq 90\%$ test accuracies get classified into
_Coexistence_ , _Partial Inversion_ and _Full Inversion_ based on their
training accuracy (specifically, their ability to memorize corrupted
examples). _Coexistence_ phase has $\geq 90\%$ training accuracy, since it
memorizes most of the corrupted training examples. _Partial Inversion_ phase
has $<90\%$ but $\geq 100(1.05-\xi)\%$ training accuracy, where $\xi$ is the
corruption fraction. Note that we subtract $\xi$ from $1.05$ instead of $1.00$
to account for the $1-2\%$ random guessing accuracy as well as statistical
fluctuations around this value. _Full Inversion_ phase has $<100(1.05-\xi)\%$
training accuracy, since it does not memorize any of the corrupted examples.
_Memorization_ phase has $<90\%$ test accuracy but $\geq 90\%$ training
accuracy, indicating that the network memorizes all of the training data,
including the corrupted and uncorrupted examples. _Forgetting_ or _confusion_
phase has $<90\%$ test as well as training accuracies. the network neither
generalizes nor manages to memorize all the training examples. We distinguish
_Forgetting_ from _confusion_ using the training curve. In the _Forgetting
phase_ the grokked network loses its performance due to steadily decreasing
weight-norms. Whereas in the _confusion_ phase, the network never reaches high
performance.

| Coexistence |  Partial Inversion |  Full Inversion | Memorization |  Forgetting / Confusion
---|---|---|---|---|---
Test accuracy | $\geq 90\%$ | $\geq 90\%$ | $\geq 90\%$ | $<90\%$ | $<90\%$
Training accuracy | $\geq 90\%$ |  $<90\%$ and $\geq 100(1.05-\xi)\%$ | $<100(1.05-\xi)\%$ | $\geq 90\%$ | $<90\%$
Table 2: Quantitative characterization of phases.

##  Appendix C Analytical Solution

We show that the network with periodic weights stated in Equation 2 is indeed
the analytical solution for the modular arithmetic task. The network output is
given by

$\displaystyle f_{q}(m,n)$
$\displaystyle=\sum_{k=1}^{N}W_{qk}(U_{km}+V_{kn})^{2}$
$\displaystyle=\frac{2}{N}\cos{\left(-\frac{2\pi}{p}\sigma(k)q-(\phi^{(u)}_{k}+\phi^{(v)}_{k})\right)}\cdot$
$\displaystyle\qquad\qquad\cdot\left(\cos{\left(\frac{2\pi}{p}\sigma(k)m+\phi^{(u)}_{k}\right)}\cos{\left(\frac{2\pi}{p}\sigma(k)n+\phi^{(v)}_{k}\right)}\right)^{2}$
$\displaystyle=\frac{2}{N}\sum_{k=1}^{N}\left\\{\frac{1}{4}\cos{\left(\frac{2\pi}{p}\sigma(k)(2m-q)+\phi^{(u)}_{k}-\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{4}\cos{\left(\frac{2\pi}{p}\sigma(k)(2m+q)+3\phi^{(u)}_{k}+\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{4}\cos{\left(\frac{2\pi}{p}\sigma(k)(2n-q)-\phi^{(u)}_{k}+\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{4}\cos{\left(\frac{2\pi}{p}\sigma(k)(2n+q)+\phi^{(u)}_{k}+3\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\boxed{\frac{1}{2}\cos{\left(\frac{2\pi}{p}\sigma(k)(m+n-q)\right)}}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{2}\cos{\left(\frac{2\pi}{p}\sigma(k)(m+n+q)+2\phi^{(u)}_{k}+2\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{2}\cos{\left(\frac{2\pi}{p}\sigma(k)(m-n-q)-2\phi^{(v)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\frac{1}{2}\cos{\left(\frac{2\pi}{p}\sigma(k)(m-n+q)+2\phi^{(u)}_{k}\right)}\right.$
$\displaystyle\left.\qquad\qquad\quad+\cos{\left(-\frac{2\pi}{p}\sigma(k)q-\phi^{(u)}_{k}-\phi^{(v)}_{k}\right)}\right\\}\,.$
(7)

We have highlighted the term that will give us the desired output with a box.
Note that the desired term is the only one that does not have additive phases
in the argument. Recall that the phases $\phi_{k}^{(u)}$ and $\phi_{k}^{(v)}$
randomly chosen – uniformly i.i.d. sampled between $(-\pi,\pi]$. Consequently,
as $N$ becomes large, all other terms will vanish due to random phase
approximation. The only surviving term will be the boxed term. We can write
this boxed term in a more suggestive form to make the analytical solution
apparent.

$\displaystyle
f_{q}(m,n)={\frac{1}{N}\sum_{k=1}^{N}\cos{\left(\frac{2\pi}{p}\sigma(k)(m+n-q)\right)}}\approx\delta^{p}(m+n-q)\,,$
(8)

where we have defined the _modular Kronecker Delta function_
$\delta^{p}(\cdot)$ as Kronecker Delta function up to integer multiples of the
modular base $p$.

$\displaystyle\delta^{p}(x)=\begin{cases}1&\;x=rp\qquad(r\in\mathbb{Z})\\\
0&\;\text{otherwise}\end{cases}\,,$ (9)

where $\mathbb{Z}$ denotes the set of all integers. Note that
$\delta^{p}(m+n-q)$ are the desired one_hot encoded labels for the Modular
Arithmetic task, by definition. Thus our network output with periodic weights
is indeed a solution.

##  Appendix D Forgetting phase

Figure 7: _Forgetting_ phase (data fraction $\alpha=0.5$, corruption fraction
$\xi=0.35$, weight decay = 20). (a) Training curves. (b) Frobenius norms of
network weights with training.

We saw in Figure 3(a) (right-most panel, wd=20.0) that extremely high weight
decay leads to _Forgetting_ phase. This phase has random guessing accuracies
on both training and test datasets at the end of training. Naively, this would
seem similar to the _Confusion_ phase. However, Upon close inspection of their
training dynamics, we find rich behaviour. Namely, the network initially
undergoes dynamics similar to _Full Inversion_ , i.e., the test accuracy
reaches $100\%$ while the train accuracy plateaus around $100(1-\xi)\%$.
However, both accuracies subsequently decay down to random guessing. This
curious phenomenon is explained by the effect of the high weight decay on the
network weights. We notice that the Frobenius norms of the weights of the
network steadily decay down to arbitrarily small values.

This curious behaviour only occurs when we use activation functions of degree
higher than 1 (e.g. $\phi(x)=x^{2}$ in this case). (More precisely, we mean
that the first term in the series expansion should be of degree at least 2.)
This can be intuitively explained by examining the gradients carefully. With
quadratic activation functions and MSE loss, gradients take the following form

$\displaystyle\begin{aligned} \frac{\partial\mathcal{L}}{\partial
U_{ab}}&=\frac{4}{|D|p}\sum_{(b,n)\in
D}\sum_{q}^{p}\left(\sum_{k=1}^{N}W_{qk}\left(U_{kb}+V_{kn}\right)^{2}-\delta^{p}(b+n-q)\right)W_{qa}\left(U_{ab}+V_{an}\right)\\\
\frac{\partial\mathcal{L}}{\partial V_{ab}}&=\frac{4}{|D|p}\sum_{(m,b)\in
D}\sum_{q}^{p}\left(\sum_{k=1}^{N}W_{qk}\left(U_{km}+V_{kb}\right)^{2}-\delta^{p}(b+n-q)\right)W_{qa}\left(U_{am}+V_{ab}\right)\\\
\frac{\partial\mathcal{L}}{\partial W_{cd}}&=\frac{4}{|D|p}\sum_{(m,n)\in
D}\left(\sum_{k=1}^{N}W_{ck}\left(U_{km}+V_{kn}\right)^{2}-\delta^{p}(m+n-q)\right)\left(U_{dm}+V_{dn}\right)^{2}\,,\end{aligned}$
(10)

where $D$ is the training dataset and $|D|$ is its size.

For non-generalizing solutions, the term inside the bracket is of degree 0 in
weights, due to the delta function term. But close to generalizing solutions,
the term inside the bracket becomes effectively degree 3 in weights, since the
delta function term cancels with the network output. Consequently, around the
generalizing solution, the weight gradients have a higher degree dependence on
the weights themselves. This makes such solutions more prone to the collapse
due to high weight decay. This serves as a heuristic justification as to why
only grokked networks seemed to be prone to _Forgetting_ , while memorizing
networks seem to be immune.

If we use vanilla Gradient Descent optimizer, the updates would plateau once
the balance between gradient updates and weight decay contribution is reached.
But this is not the case in Adaptive optimizers. Indeed, we find that this
phenomenon exclusively occurs with Adaptive optimizers (e.g. Adam in this
case).

We believe that it is possible to get further insight into this phase by
closely examining the optimization landscape and network dynamics. Since this
behaviour is exclusive to grokked networks, understanding this phase would
lead to further understanding of grokking dynamics. Alternatively, a better
understanding of the dynamics of grokking should give us insights into this
curious phenomenon. We defer this analysis for future work and welcome an
examination by our readers.

##  Appendix E Comparison of BatchNorm, LayerNorm and minibatch w/o
Normalization

Here we examine the effects of BatchNorm, LayerNorm and mini-batch noise. We
find that generalization and _Inversion_ are indeed enhanced by decreasing the
batch-size. This demonstrates that implicit regularization from mini-batch
noise favours simpler solutions in our setup.

Additionally, upon comparing Figure 8, Figure 9 and Figure 10 we find enhanced
_Inversion_ in cases with normalization layers (especially BatchNorm). For
BatchNorm, we provide an explanation in the main text, with the the BatchNorm
weights de-amplifying memorizing neurons in favour of generalizing ones. For
LayerNorm, however, we do not find such an effect.

Figure 8: Phase diagram with BatchNorm (post-BN) Figure 9: Phase diagram with
mini-batch training, without any explicit regularization or normalization
layers. Figure 10: Phase diagram with LayerNorm (post-LN)

##  Appendix F Additional Training Curves

Figure 11: Training curves with dropout (a) dp=0.0 (b) dp=0.3 (c) dp=0.6 (d)
dp=0.9 Figure 12: Training curves in various phases with BatchNorm. All
networks are trained with corruption fraction $\xi=0.3$ and batch-size = 64;
with variable data fractions. (a) Coexistence ($\alpha=0.4$) (b) Partial
Inversion ($\alpha=0.5$) (c) Full Inversion ($\alpha=0.8$) Figure 13: Training
curves in various phases with LayerNorm. All networks are trained with
corruption fraction $\xi=0.3$ and batch-size = 64; with variable data
fractions. (a) Coexistence ($\alpha=0.4$) (b) Partial Inversion ($\alpha=0.5$)
(c) Full Inversion ($\alpha=0.8$) Figure 14: Training curves in various phases
without any normalization, but with mini-batches. All networks are trained
with corruption fraction $\xi=0.3$ and batch-size = 64; with variable data
fractions. (a) Coexistence ($\alpha=0.4$) (b) Partial Inversion ($\alpha=0.5$)
(c) Full Inversion ($\alpha=0.8$)

##  Appendix G Detailed Phase Diagrams

Here we replot the phase diagrams in Figure 3 together with train and test
accuracy.

Figure 15: Modular Addition with weight decay. The first two rows of color-
plots are train and test accuracies, respectively. Figure 16: Modular Addition
with Dropout. The first two rows of color-plots are train and test accuracies,
respectively.

##  Appendix H Choice of Loss Function

CrossEntropy loss is known to be sensitive to label noise (Ghosh et al.,
2017). By adding regularization, we recover similar phase diagrams. The
requirement for regularization is much greater for CrossEntropy loss compared
to MSE. We find that CrossEntropy training is very sensitive to weight decay,
with _Inversion_ occuring for a narrow band of weight decay values, beyond
which we find _Confusion_. We do not observe a _Forgetting_ phase in this
case. See Figure 17. With Dropout, we find that CrossEntropy training requires
very high dropout probabilities to grok and/or attain _Inversion_.

(a) Weight decay

(b) Dropout

Figure 17: Modular Addition trained with CrossEntropy Loss. We see that the
robustness is getting weaker compared to MSE.

##  Appendix I Choice of Activation Function

To check if the story depends on the activation function, we did the modular
addition experiments for 2-layer ReLU network, see Figures 18 and 19. The
requirement for regularization for grokking is higher compared to quadratic
activation function. We noticed that the phase diagram for weight decay
behaves similarly to quadratic activation functions, except for the
disappearance of the _Forgetting_ phase. Dropout does not seem as effective
for ReLU networks, with the network entering _Confusion_ phase for even
moderate dropout probabilites.

Figure 18: Phase diagram with weight decay, ReLU. Figure 19: Phase diagram
with Dropout, ReLU.

##  Appendix J Effect of of width

Figure 20: Phase diagram for width $N=5000$ fully connected network with
quadratic activation. Regularization with weight decay. Compared to $N=500$
network in Figure 3(a), the network generalizes marginally poorly. Note: the
absence of _Partial Inversion_ phase and smaller _Coexistence_ phase with
$wd=0.0$. This can be attributed to the increased capacity to memorize due to
abundance of Neurons. Figure 21: Phase diagram for width $N=5000$ fully
connected network with ReLU activation. Regularization with weight decay. In
this case, the network performance is drastically improved compared to $N=500$
(Figure 18). The confusion phase is largely eliminated, due to the increased
capacity of the network, as a result of abundance of neurons.

##  Appendix K Other Modular Arithmetic Datasets

In this section, we check if the phase diagram story is unique to modular
addition. We repeat our experiments with modular multiplication
($z=(mn)\%p\equiv mn\,\,\textrm{mod}\,\,p$) and obtain almost identical phase
diagrams to modular addition (see Figure 22).

Figure 22: Phase diagram for Modular Multiplication, weight decay. Note that
the phase diagram is identical to Modular Arithmetic case – representative of
the similarity in the algorithms learnt by the network.


# Third-order intrinsic anomalous Hall effect with generalized semiclassical
theory

Longjun Xiang  College of Physics and Optoelectronic Engineering, Shenzhen
University, Shenzhen 518060, China  Chao Zhang  Department of Physics, The
Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong, China  Luyang Wang  College of Physics and Optoelectronic Engineering,
Shenzhen University, Shenzhen 518060, China  Jian Wang<EMAIL_ADDRESS>College of Physics and Optoelectronic Engineering, Shenzhen University,
Shenzhen 518060, China  Department of Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong, China

###### Abstract

The linear intrinsic anomalous Hall effect (IAHE) and second-order IAHE have
been intensively investigated in time-reversal broken systems. However, as one
of the important members of the nonlinear Hall family, the investigation of
third-order IAHE remains absent due to the lack of an appropriate theoretical
approach, although the third-order extrinsic AHE has been studied within the
framework of first- and second-order semiclassical theory. Herein, we
generalize the semiclassical theory for Bloch electrons under the uniform
electric field up to the third-order using wavepacket method and based on
which we predict that the third-order IAHE can also occur in time-reversal
broken systems. Same as the second-order IAHE, we find the band geometric
quantity, the second-order field-dependent Berry curvature arising from the
second-order field-induced positional shift, plays a pivotal role to observe
this effect. Moreover, with symmetry analysis, we find that the third-order
IAHE, as the leading contribution, is supported by $15$ time-reversal broken
3D magnetic point groups (MPGs), corresponding to a wide class of
antiferromagnetic (AFM) materials. Guided by the symmetry arguments, a two-
band model is chosen to demonstrate the generalized theory. Furthermore, the
generalized third-order semiclassical theory depends only on the properties of
Bloch bands, implying that it can also be employed to explore the IAHE in
realistic AFM materials, by combining with first-principles calculations.

††preprint: APS/123-QED

Introduction.— As one of the important phenomena in condensed matter physics,
the charge Hall effect and its variants, such as the anomalous Hall effect
Karplus ; RevModPhys.82.1539  and spin Hall effect SCZhang ; Sinova2004 , have
been studied continuously since its discovery in 1879 Hall . For a long time,
it has been believed that the observation of charge Hall current requires
broken time-reversal ($\mathcal{T}$) symmetry landau1999statistical . However,
within the framework of the first-order semiclassical theory
RevModPhys.82.1959 , Sodemann and Fu predicted PhysRevLett.115.216806  that
the second-order nonlinear Hall effect (NLHE) can exist even in
$\mathcal{T}$-invariant materials, driven by the band geometric quantity Berry
curvature dipole (BCD). Recently, this extrinsic (proportional to the
relaxation time $\tau$) NLHE has been observed experimentallyNature2019 ;
Raffaele ; NRP2021 .

In fact, slightly before investigating NLHEs within $\mathcal{T}$-invariant
materials, the second-order nonlinear intrinsic (free of $\tau$ and depends
only on the band topology) anomalous Hall effect (IAHE) in magnetic systems
has been predicted by Y. Gao et al.PhysRevLett.112.166601 , using the second-
order semiclassical theory developed by them. Very recently, with this second-
order semiclassical theory, two research groups further predicted that the
second-order nonlinear IAHE can be observed in $\mathcal{PT}$-symmetric
compensated antiferromagnetic (AFM) materials CuMnAs wang2021intrinsic  and
Mn2Au liu2021intrinsic . Furthermore, using the second-order semiclassical
theory, the extrinsic third-order NLHE, have been investigated in materials
lai2021third ; PhysRevB.105.045118 ; TaIrTe4  with $\mathcal{T}$ and
$\mathcal{P}$ symmetries, suggesting that the third-order Hall signal can be
dominant in some materials. Besides, the extrinsic third-order nonlinear AHE,
induced by Berry curvature quadrupole, can also be the leading contribution in
some $\mathcal{T}$-broken materials, predicted by Zhang et al.
zhang2020higherorder  with the first-order semiclassical theory. However, as
one of the important members of the nonlinear intrinsic anomalous Hall family,
the third-order IAHE remains unexplored, as far as we know.

In this work, we generalize the semiclassical theory for Bloch electrons under
the uniform electric field to third-order using the wavepacket method and
based on which we predict the existence of third-order IAHE in
$\mathcal{T}$-broken systems. Like the second-order semiclassical theory
PhysRevLett.112.166601 , we find the band geometric quantity, the second-order
field dependent Berry curvature (BC) arising from the second-order field-
induced positional shift PhysRevLett.112.166601 , plays a key role in
observing this effect. Moreover, from symmetry analysis, we find that the
third-order IAHE, as the leading contribution, can be hosted by $15$
$\mathcal{T}$-broken 3D magnetic point groups (MPGs), corresponding to a wide
class of AFM materials. Following the symmetry arguments, a two-band model is
chosen to illustrate our generalized theory.

The third-order semiclassical theory.— Using the semiclassical wavepacket
approach RevModPhys.82.1959 ; PhysRevLett.112.166601 , the equation of motion
(EOM) describing the dynamics of the center of wavepacket under the uniform
electric field $\bm{E}$ can be obtained:

$\displaystyle\dot{\bm{r}}_{c}$
$\displaystyle=\partial_{\bm{k}_{c}}\bar{\epsilon}-\dot{\bm{k}}_{c}\times\bm{\Omega}(\bm{E})$
(1) $\displaystyle\dot{\bm{k}}_{c}$ $\displaystyle=-\bm{E}$ (2)

where $\bm{k}_{c}/\bm{r}_{c}$ represents the momentum/position center of the
wavepacket, for a paticular band $n$ $\bar{\epsilon}$ is the band energy
including the correction due to the external electric field,
$\partial_{\bm{k}_{c}}\equiv\partial/\partial\bm{k}_{c}$, and
$\bm{\Omega}(\bm{E})\equiv\nabla_{\bm{k}_{c}}\times\bm{\mathcal{A}}_{nn}({\bm{E}})$
is the field-dependent BC with $\bm{\mathcal{A}}_{nn}(\bm{E})$ the field-
dependent Berry connection, or field-induced positional shift. In particular,
by calculating the first-order field-induced positional shift for the position
center $\bm{r}_{c}$ with wavepacket approach, Y. Gao et al.
PhysRevLett.112.166601  has successfully developed the second-order
semiclassical theory. Herein, we further calculate the second-order positional
shift to formulate the third-order semiclassical theory.

To that purpose, we focus on the $0$th band and construct the wavepacket
PhysRevLett.112.166601  as follows:

$\displaystyle|\Psi\rangle$
$\displaystyle=\int_{\bm{p}}e^{i\bm{p}\cdot\bm{r}}\left(C_{0}(\bm{p})|u_{0}(\bm{p})\rangle+\sum_{n\neq
0}C_{n}(\bm{p})|u_{n}(\bm{p})\rangle\right)$ (3)

where $|u_{n}(\bm{p})\rangle$ is the Bloch state, $C_{0}$ is the zeroth-order
amplitude with $|C_{0}|^{2}=\delta(\bm{p}-\bm{p}_{c})$ to normalize the
wavepacket up to the first-orderRevModPhys.82.1959 ,
$C_{n}=C_{n}^{(1)}+C_{n}^{(2)}$ with $C_{n}^{(1)}$/$C^{(2)}_{n}$ the first-
order/second-order amplitude in terms of the external electric field $\bm{E}$.
Note that both $C_{n}^{(1)}$ and $C^{(2)}_{n}$ are related to $C_{0}$, for
example PhysRevLett.112.166601 ,

$\displaystyle C_{n}^{(1)}$ $\displaystyle=\mathcal{M}_{n0}C_{0}$ (4)

where
$\mathcal{M}_{n0}=\bm{E}\cdot\bm{\mathcal{A}}_{n0}/(\epsilon_{0}-\epsilon_{n})$
with $\bm{\mathcal{A}}_{n0}=\langle u_{n}|i\partial_{\bm{k}_{c}}|u_{0}\rangle$
the interband Berry connection and $\epsilon_{0/n}$ the unperturbed band
energy. Eq.(4) is obtained by solving the time-dependent Schrödinger equation,
therefore the relation between $C_{n}^{(2)}$ and $C_{0}$ can be identified
similarly note1 :

$\displaystyle C^{(2)}_{n}$
$\displaystyle=\dfrac{\bm{E}\cdot[(i\partial_{\bm{p}}-\bm{r}_{c})C_{0}]\mathcal{M}_{n0}}{\epsilon_{0}-\epsilon_{n}}+\dfrac{1}{2}\sum_{m\neq
0}\mathcal{M}_{0m}\mathcal{M}_{m0}C_{0}\delta_{0n}$
$\displaystyle+\left(\sum_{m\neq
0}\frac{\bm{E}\cdot\bm{\mathcal{A}}_{nm}\mathcal{M}_{m0}}{\epsilon_{0}-\epsilon_{n}}+\frac{\bm{E}\cdot
i\partial_{\bm{p}}\mathcal{M}_{n0}}{\epsilon_{0}-\epsilon_{n}}\right)C_{0}$
(5)

where the second term is due to the normalization of wavepacket up to the
second-order.

Once the wavepacket is obtained, the position center $\bm{r}_{c}$ can be
expressed as GaoYangLT :

$\displaystyle\bm{r}_{c}\equiv\langle\Psi|\bm{r}|\Psi\rangle=\partial_{\bm{k}_{c}}\gamma+\bm{\mathcal{A}}_{0}+\bm{\mathcal{A}}_{0}^{(1)}+\bm{\mathcal{A}}_{0}^{(2)}$
(6)

where $\gamma=-\text{arg}(C_{0})$, $\bm{\mathcal{A}}_{0}=\langle
u_{0}|i\partial_{\bm{k}_{c}}|u_{0}\rangle$ is the intraband Berry connection,
$\bm{\mathcal{A}}^{(1)}_{0}=2\text{Re}\sum_{n\neq
0}\bm{\mathcal{A}}_{0n}\mathcal{M}_{n0}$ is the first-order positional shift,
firstly derived by Y. Gao et al.PhysRevLett.112.166601 . At this stage, by
taking the curl of $\bm{\mathcal{A}}_{0}$ and $\bm{\mathcal{A}}_{0}^{(1)}$,
the conventional BC and the first-order field-dependent BC
$\bm{\Omega}^{(1)}(\bm{E})$ can be obtained and hence the first-order and
second-order semiclassical theory are formulated. Furthemore, with the second-
order wavepacket constructed in Eq.(3), we obtain the second-order positional
shift note2 :

$\displaystyle\bm{\mathcal{A}}_{0}^{(2)}$ $\displaystyle=\text{Re}\sum_{n\neq
0}^{m\neq
0}\left[\dfrac{2\bm{\mathcal{A}}_{0n}\bm{E}\cdot\bm{\mathcal{A}}_{nm}\mathcal{M}_{m0}}{\epsilon_{0}-\epsilon_{n}}-\mathcal{M}_{0n}\bm{\mathcal{A}}_{nm}\mathcal{M}_{m0}\right]$
$\displaystyle-\text{Re}\sum_{n\neq
0}\left[\dfrac{2\bm{E}\cdot\bm{\mathcal{A}}_{00}\mathcal{M}_{n0}\bm{\mathcal{A}}_{0n}}{\epsilon_{0}-\epsilon_{n}}-\bm{\mathcal{A}}_{00}\mathcal{M}_{0n}\mathcal{M}_{n0}\right]$
$\displaystyle-\text{Re}\sum_{n\neq
0}\left[\mathcal{M}_{0n}i\partial_{\bm{k}_{c}}\mathcal{M}_{n0}-\dfrac{\bm{\mathcal{A}}_{0n}\bm{E}\cdot
i\partial_{\bm{k}_{c}}\mathcal{M}_{n0}}{\epsilon_{0}-\epsilon_{n}}\right]$
$\displaystyle-\text{Re}\sum_{n\neq 0}\mathcal{M}_{n0}\bm{E}\cdot
i\partial_{\bm{k}_{c}}\left(\dfrac{\bm{\mathcal{A}}_{0n}}{\epsilon_{0}-\epsilon_{n}}\right)$
(7)

Note that we have always replaced the momentum center $\bm{p}_{c}$ with
$\bm{k}_{c}$ at the final results. Importantly, by writting
$\mathcal{A}^{\alpha,(2)}_{0}\equiv
T^{0}_{\alpha\beta\gamma}E_{\beta}E_{\gamma}$, we find

$\displaystyle T^{0}_{\alpha\beta\gamma}$ $\displaystyle=\text{Re}\sum_{n\neq
0}\left[2\mathcal{A}_{00}^{\beta}M_{0n}^{\alpha}M_{n0}^{\gamma}+\mathcal{A}_{00}^{\alpha}M_{0n}^{\beta}M_{n0}^{\gamma}\right]$
$\displaystyle-\text{Re}\sum_{n\neq 0}^{m\neq
0}\left[2M_{0n}^{\alpha}\mathcal{A}_{nm}^{\beta}M_{m0}^{\gamma}+M_{0n}^{\beta}\mathcal{A}_{nm}^{\alpha}M_{m0}^{\gamma}\right]$
$\displaystyle-\text{Re}\sum_{n\neq
0}\left[M_{0n}^{\alpha}i\partial_{\beta}M_{n0}^{\gamma}-M_{n0}^{\gamma}i\partial_{\beta}M_{0n}^{\alpha}\right]$
$\displaystyle-\text{Re}\sum_{n\neq
0}M_{0n}^{\beta}i\partial_{\alpha}M_{n0}^{\gamma}$ (8)

which is the second-order Berry-connection polarizability tensor (BPT) for
band $0$. Here $\partial_{\alpha}=\partial/\partial k_{c}^{\alpha}$ and
$M_{mn}^{\alpha}=\mathcal{A}^{\alpha}_{mn}/(\epsilon_{n}-\epsilon_{m})$. Same
as the physical implication of
$\bm{\mathcal{A}}_{0}^{(1)}(\bm{E})$PhysRevLett.112.166601 ,
$\bm{\mathcal{A}}_{0}^{(2)}(\bm{E})$ stands for a second-order correction to
the conventional Berry connection $\bm{\mathcal{A}}_{0}$ of the unperturbed
band, which means that the wavepacket also acquires a shift
$\bm{\mathcal{A}}_{0}^{(2)}(\bm{E})$ in its center $\bm{r}_{c}$ of mass
position. In addition, we note that $\bm{\mathcal{A}}^{(2)}_{0}(\bm{E})$ also
respects the periodicity of the lattice due to the uniform external field,
hence it does not cause any macroscopic charge density gradient and also will
not affect the electron chemical potential profile PhysRevLett.112.166601 ;
GaoYang2019 . The second-order BPT is the central concept of our third-order
semiclassical theory.

By taking the curl of the second-order field-dependent Berry-connection, we
obtain the field-dependent BC $\bm{\Omega}^{(2)}(\bm{E})$ at the same order.
Furthermore, substituting $\bm{\Omega}^{(2)}(\bm{E})$ into Eq.(1), the third-
order semiclassical theory is established when the band energy
$\bar{\epsilon}$ is corrected to third-ordernote3 . Like the second-order
semiclassical theory, the third-order EOM and the second-order field-dependent
BC, plays an essential role in investigating the IAHE in AFM materials, as
will be illustrated below, especially when the linear and second-order IAHE
signals are forbidden by symmetry. The third-order semiclassical theory, with
the second-order field-dependent BC originating from the second-order field-
induced positional shift as the band geometric quantity, is our first main
result.

The third-order IAHE. — Within the framework of third-order semiclassical
transport theory, if we ignore the scattering effects arising from impurities,
the third-order intrinsic Hall current density can be expressed as
$J_{\alpha}=\chi^{\text{int}}_{\alpha\beta\gamma\eta}E_{\beta}E_{\gamma}E_{\eta}$,
where

$\displaystyle\chi_{\alpha\beta\gamma\eta}^{\text{int}}=\int_{k}\Lambda_{\alpha\beta\gamma\eta}f_{0}$
(9)

is the third-order intrinsic Hall conductivity, which is a rank-$4$ tensor. In
Eq.(9), $f_{0}$ is the equilibrium Fermi distribution function, and

$\displaystyle\Lambda_{\alpha\beta\gamma\eta}=\sum_{n}\left[\partial_{\beta}T_{\alpha\gamma\eta}^{n}-\partial_{\alpha}T_{\beta\gamma\eta}^{n}\right]$
(10)

is the integrand for the third-order intrinsic Hall conductivity, which is
antisymmetric in its first two indices and symmetric in its last two indices,
namely
$\chi^{\text{int}}_{\alpha\beta\gamma\eta}=-\chi^{\text{int}}_{\beta\alpha\gamma\eta}$
and
$\chi^{\text{int}}_{\alpha\beta\gamma\eta}=\chi^{\text{int}}_{\beta\alpha\eta\gamma}$.
By performing an integration by parts for Eq.(9), we find that the third-order
IAHE is also a Fermi liquid property PhysRevLett.115.216806 ;
PhysRevLett.112.166601 ; wang2021intrinsic ; liu2021intrinsic ; Haldane , as
expected.

Table 1:  The classification for $90$ $\mathcal{T}$-broken MPGs in 3D with the third-, and fourth-order IAHEs, as the leading contribution. The order of IAHE |  The MPGs for IAHE
---|---
Third-order |  $mmm$, $4^{\prime}/m$, $4/mmm$, $4^{\prime}/mm^{\prime}m$, $\bar{3}m$, $6^{\prime}$, $6^{\prime}/m^{\prime}$, $6^{\prime}22^{\prime}$, $6^{\prime}mm^{\prime}$, $\bar{6}m2$ $6/mmm$, $6^{\prime}/m^{\prime}mm^{\prime}$, $m\bar{3}$, $4^{\prime}32^{\prime}$, $m\bar{3}m^{\prime}$
Fourth-order |  $6^{\prime}/m$, $6^{\prime}/mmm^{\prime}$, $\bar{4}3m$, $m^{\prime}\bar{3}^{\prime}m$

Symmetry analysis.— Next, we investigate what kind of symmetry will host a
non-vanishing third-order intrinsic Hall current and when this response
becomes the leading contribution. It has been well known that the number of
independent components of a physical quantity such as conductivity tensor is
dictated by magnetic point group (MPG) symmetry of the system, as required by
Neumann’s principle Neumann ; anisotropy . Under $\mathcal{T}$-symmetry, the
field dependent BC is $\mathcal{T}$-odd and hence there is no intrinsic Hall
signal in $\mathcal{T}$-invariant systems. However, for materials with
$\mathcal{T}$-broken MPGs, the BC is nonzero and hence we can observe the
intrinsic Hall signal in these systems.

In particular, for the rank-4 IAHE conductivity tensor
$\chi_{\alpha\beta\gamma\eta}^{\text{int}}$, the constraint imposed by MPG
symmetry operations $R$ and $R\mathcal{T}$ can be expressed as anisotropy :

$\displaystyle\chi^{\text{int}}_{\alpha\beta\gamma\eta}=\eta_{T}R_{\alpha\alpha^{\prime}}R_{\beta\beta^{\prime}}R_{\gamma\gamma^{\prime}}R_{\eta\eta^{\prime}}\chi^{\text{int}}_{\alpha^{\prime}\beta^{\prime}\gamma^{\prime}\eta^{\prime}}$
(11)

where $\eta_{T}=1~{}(\eta_{T}=-1)$ is for $R(R\mathcal{T})$ and
$R_{\alpha\alpha^{\prime}}$ is the matrix element of the spatial point group
operation $R$. Starting from Eq.(11), in principle we can classify the $90$ 3D
MPGs without $\mathcal{T}$ (ruled out $32$ grey MPGs in total $122$ MPGs in
3D) using the linear, second-order, third-order and other high-order IAHEs.
For example, considering the MPG $6^{\prime}$ with the generator
$C_{6}\mathcal{T}$, we find
$\chi^{\text{int}}_{xzxx}=-\chi^{\text{int}}_{zxxx}=\chi^{\text{int}}_{zxyy}=-\chi^{\text{int}}_{xzyy}\neq
0$ and
$\chi^{\text{int}}_{xzxy}=-\chi^{\text{int}}_{zxxy}=\chi^{\text{int}}_{xzyx}=-\chi^{\text{int}}_{zxyx}\neq
0$, whereas both the rank-$2$ and rank-3 conductivity tensors
$\sigma^{\text{int}}_{\alpha\beta}$ and
$\sigma^{\text{int}}_{\alpha\beta\gamma}$ vanish for linear IAHE and second-
order IAHE, respectively.

On the other hand, we can also define the Jahn’s notations Jahn
$a\\{V^{2}\\}$, $a\\{V^{2}\\}V$, $a\\{V^{2}\\}[V^{2}]$ and
$a\\{V^{2}\\}[V^{3}]$ for linear, second-order, third-order and fourth-order
IAHE conductivity tensors, respectively, and then use Bilbao Crystallographic
Server Bilbao  to find the 3D MPGs with a non-vanishing Hall signal note4  and
identify the leading contribution, as shown in TABLE I. We find that the the
linear, second-order, third-order, and fourth-order IAHEs are supported by
$31$, $39$, $15$, and $4$ 3D MPGs, respectively, which means that one can
observe Hall signal in almost all of $\mathcal{T}$-broken crystals (The
leading contribution for linear and second-order IAHEs can be found in TABLE I
of supporting materials, in which a full clsssifiction rather than leading-
order is also given.) It should be noted that the MPGs hosting
linear/nonlinear IAHE, are/aren’t compatible with ferromagnetism, therefore,
the high-order Hall current response, as the leading contribution, can only be
observed in AFM materials. We also note that the linear IAHE are normally
assumed to be proportional to magnetization such as in FM metals, but
theoretical predictions and experimental observations have recognized that
large Hall effects can also occur in non-collinear or collinear AFM crystals
Chenhua ; Christoph ; Satoru ; Suzuki , offered by neither of the global
magnetic-dipole symmetry-breaking mechanisms SciAdv ; naturereview . The
symmetry arguments to search for the third-order IAHE in $\mathcal{T}$-broken
systems, as the leading contribution, is our second main result.

Two-band model.— In this section, we employ a two-band model to illustrate our
general theory. To that purpose, let us first consider the typical transport
platform with a planar geometry, in which the applied electric field $\bm{E}$
and the generated current both within the plane (denoted as $xy$ plane without
loss of generality), and assuming the $\bm{E}$ field forming an angle $\theta$
with the principal axis $x$ of the crystal, namely
$\bm{E}=E(\cos\theta,\sin\theta)$, we find that the in-plane intrinsic third-
order anomalous Hall current density can be calculated as note5 :

$\displaystyle J^{(3)}_{\text{AH}}$
$\displaystyle=E^{3}\left[\chi^{\text{int}}_{yxxx}\cos^{2}\theta-\chi^{\text{int}}_{xyyy}\sin^{2}\theta\right]$
$\displaystyle+\dfrac{1}{2}E^{3}\sin{2\theta}[\chi^{\text{int}}_{yxxy}+\chi^{\text{int}}_{yxyx}]$
(12)

Meanwhile, we note that the in-plane intrinsic third-order parallel current
$J^{(3)}_{||}$ is zero note5 . For simplicity, below we will take $\theta=0$
and hence $J_{\text{AH}}^{(3)}=\chi_{yxxx}^{\text{int}}E^{3}$.

Following the full classification in the supporting materials note4 , the two-
band model respecting MPG $4^{\prime}m^{\prime}m$

$\displaystyle
H(\bm{k})=tk^{2}+v(k_{y}\sigma_{x}-k_{x}\sigma_{y})+m(k_{x}^{2}-k_{y}^{2})\sigma_{z}$
(13)

is chosen to demonstrate the third-order IAHE, which is a low energy effective
Hamiltonian to describe the band behavior around $\Gamma$ point of monolayer
SrMnBi2 zhang2020higherorder . In Eq.(13), $\sigma_{i}(i=x,y,z)$ represent the
Pauli matrices acting on the spin, and $k^{2}=k_{x}^{2}+k_{y}^{2}$. The band
dispersion for this model is $\epsilon_{\pm}=tk^{2}\pm h$ with
$h^{2}=v^{2}k^{2}+m^{2}(k_{x}^{2}-k_{y}^{2})$ even in $k_{x}$ and $k_{y}$, as
shown in FIG.1(a), which resembles the Rashba dispersion but with a
$\mathcal{T}$-broken second-order warping term zhang2020higherorder .
Interestingly, this effective model also has been utilized to domonstrate the
existence of the third-order extrinsic anomalous Hall effect, as the leading
contribution induced by Berry quadrupole zhang2020higherorder .

Figure 1:  (a) The band structure along $k_{x}$ direction for the two-band
model. (b-d) The $\bm{k}$-resolved distribution of $T_{xxx}$, $T_{yxx}$, and
$\Lambda_{yxxx}$ for the lower band of the two-band model, with $t=3,v=1,m=2$.
(e) The intrinsic third-order anomalous Hall conductivity
$\chi^{\text{int}}_{yxxx}$ versus the Fermi energy $\mu$ and different
parameters $m$.

For this model, the conventional Berry curvature for the lower band is found
to be:

$\displaystyle\Omega_{z}=-\dfrac{mv^{2}(k_{x}^{2}-k_{y}^{2})}{2h^{3}}$ (14)

which is even in $k_{x}$ and $k_{y}$. However, due to the mirror symmetry
$\mathcal{M}_{x+y}$, the integral for $\Omega_{z}$ in the Brillouin zone
vanishes and hence there is no linear IAHE. Similarly, the first-order field-
dependent BC for the lower band is given as:

$\displaystyle\Omega^{(1)}_{z}$
$\displaystyle=-\dfrac{v^{2}k_{y}\left[v^{2}+m^{2}(3k_{x}^{2}+k_{y}^{2})\right]}{2h^{5}}E_{x}$
(15)

which is odd in $k_{y}$ and hence makes no contribution to IAHE at the second-
order. Note that this result is not inconsistent with our full classification,
in which the system with MPG $4^{\prime}m^{\prime}m$ can hold a leading-order
second-order IAHE. In fact, by adding a term
$k_{+}^{3}\sigma_{+}+k_{-}^{3}\sigma_{-}$ Fu2009  with MPG $41^{\prime}$ which
maintains the symmetry of the system, we indeed obtain a nonzero second-order
intrinsic anomalous Hall signal. Here we focus on the third-order nonlinear
transport behavior of this model. In terms of our third-order semiclassical
theory, the second-order field-dependent BC can be calculated as (assume
$E_{x}\equiv 1$):

$\displaystyle\Omega^{(2)}_{z}=\partial_{x}T_{yxx}-\partial_{y}T_{xxx}=\dfrac{g_{1}m^{3}v^{2}+g_{2}mv^{4}}{8h^{7}}$
(16)

with $T_{yxx/xxx}$ the second-order BPTs and $g_{1/2}$ even functions of
$k_{x}$ and $k_{y}$ note6 . Therefore, we conclude that the second-order
field-dependent BC will lead to a nonzero third-order Hall current response.
More intuitively, we present the $\bm{k}$-resolved distribution for $T_{xxx}$,
$T_{yxx}$, and $\Lambda_{yxxx}$ (same as the $\Omega^{(2)}_{z}$) for the lower
band, as shown in FIG.1(b-d), which determine the final result of
$\chi_{yxxx}^{\text{int}}$. Particularly, we find that both $T_{xxx}$ and
$T_{yxx}$ exhibit a dipole pattern in momentum space, respectively, and the
resultant integrand $\Lambda_{yxxx}$ approximately features a negative ellipse
landscape around $\Gamma$ point. In FIG. 1(e), we calculate the third-order
IAHE conductivity $\chi^{\text{int}}_{yxxx}$ at zero temperature, as a
function of the scaled chemical potential for different parameters $m$, from
which we conclude that the third-order response is nonzero and becomes
significant when the chemical potential is close to zero energy, manifesting
the fact that the third-order intrinsic conductivity is concentrated around
the small-gap region liu2021intrinsic . As a check, we find the third-order
IAHE conductivity $\chi^{\text{int}}_{yxxx}$ will disappear when $m=0$, namely
when the $\mathcal{T}$-symmetry is recovered, as can also be easily seen from
Eq.(16).

Summary.— In this work, we develop the third-order semiclassical theory for
Bloch electrons under the uniform external electric field with the
semiclassical wavepacket approach. As one of the important applications, we
predict that the third-order IAHE, driven by the band geometric quantity—the
second-order BC arising from the second-order field-induced positional shift,
can occur in $\mathcal{T}$-broken systems. It should be emphasized that the
third-order IAHE, as an important member of the nonlinear Hall family, has not
been explored so far due to the lack of an appropriate theoretical approach.
Furthermore, with symmetry arguments, we find that almost all the MPGs without
time-reversal symmetry can be classified by linear, second-order, third-order
and also fourth-order IAHEs. Importantly, we find the intrinsic third-order
nonlinear anomalous Hall signal, as the leading contribution, can be
accommodated by $15$ 3D MPGs and hence fill the gap in previous studies,
especially in AFM spintronics. Finally, a two-band toy model is employed to
demonstrate the generalized theory.

Last but not least, we note that our third-order semiclassical theory relies
only on the properties of Bloch bands, which indicates that our theory can be
combined with first-principles calculations to explore the IAHE of the
realistic AFM materials. For example, following our symmetry analysis, the AFM
materials MnTe MnTe  and CoNb3S6 CoNb3S6  should exhibit a leading-order
three-order IAHE signal, which will be explored in futural works.

## Acknowledgements

This work was financially supported by the Natural Science Foundation of China
(Grant No.12034014 and No. 12004442) and Guangdong Basic and Applied Basic
Research Foundation (Grants No. 2021B1515130007).

## References

  * (1) R. Karplus and J. M. Luttinger, Phys. Rev. 95, 1154 (1954).
  * (2) N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).
  * (3) S. Murakami, N. Nagaosa, S.-C. Zhang, Science 301, 1348 (2003).
  * (4) J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, Phys. Rev. Lett. 92, 126603 (2004).
  * (5) E. Hall, American Journal of Mathematics 2, 287 (1879).
  * (6) L.D. Landau, E.M. Lifshitz and L.P. Pitaevskii, in Statistical Physics, Course of Theoretical Physics Vol. 5 3rd. ed. (Pergamon Press, Oxford, 1999)
  * (7) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * (8) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * (9) Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A.M.M. Valdivia, S. Wu, Z. Du, C.-H. Hsu, S. Fang, Q.D. Gibson, K. Watanabe, T. Taniguchi, R.J. Cava, E. Kaxiras, H.-Z. Lu, H. Lin, L. Fu, N. Gedik, and P. Jarillo-Herrero, Nature 565, 337 (2019).
  * (10) R. Battilomo, N. Scopigno, and C. Ortix, Phys. Rev. Lett. 123, 196403 (2019).
  * (11) Z. Z. Du, H.-Z. Lu, and X. C. Xie, Nat. Rev. Phys. 3, 744 (2021).
  * (12) Y. Gao, S.A. Yang, and Q. Niu, Phys. Rev. Lett. 112, 166601 (2014).
  * (13) C. Wang, Y. Gao, and D. Xiao, Phys. Rev. Lett. 127, 277201 (2021).
  * (14) H. Liu, J. Zhao, Y.-X. Huang, W. Wu, X.-L. Sheng, C. Xiao, and S. A. Yang, Phys. Rev. Lett. 127, 277202 (2021).
  * (15) S. Lai, H. Liu, Z. Zhang, J. Zhao, X. Feng, N. Wang, C. Tang, Y. Liu, K. S. Novoselov, S. A. Yang, and W. B. Gao, Nat. Nanotechnol. 16, 869 (2021).
  * (16) C. Wang, R.-C. Xiao, H. Liu, Z. Zhang, S. Lai, C. Zhu, H. Cai, N. Wang, S. Chen, Y. Deng, Z. Liu, S.A. Yang, and W.-B. Gao, National Science Review nwac020, 2095 (2022).
  * (17) H. Liu, J. Zhao, Y.-X. Huang, X. Feng, C. Xiao, W. Wu, S. Lai, W.-B. Gao, and S.A. Yang, Phys. Rev. B 105, 045118 (2022).
  * (18) C.-P. Zhang, X.-J. Gao, Y.-M. Xie, H.C. Po, and K.T. Law, arXiv:2012.15628.
  * (19) See section I (A) in the supporting materials for details.
  * (20) Y. Gao, Low. Temp. Phys. Lett. 41, 0241 (2019).
  * (21) See section I (B) in the supporting materials for details.
  * (22) Y. Gao, Frontiers of Physics 14, 1 (2019).
  * (23) See section I (D) in the supporting materials for details.
  * (24) F. D. M. Haldane, Phys. Rev. Lett. 93, 206602 (2004).
  * (25) Neumann’s principle can be stated as: if a crystal is invariant with respect to certain symmetry, any of its physical properties must also be invariant with respect to the same symmetry. As a consequence of this principle, symmetry imposes constraint on the physical quantity, making the relevant effect disappear.
  * (26) R. E. Newnham, Properties of Materials: Anisotropy, Symmetry, Structure (Oxford University Press, 2005).
  * (27) H. A. Jahn, Acta Crystallogr. 2, 30 (1949).
  * (28) S. V. Gallego, J. Etxebarria, L. Elcoro, E. S. Tasci, and J. M. Perez-Mato, Acta Crystallogr. Sect. A 75, 438 (2019).
  * (29) See section II in the supporting materials for the full classification.
  * (30) H. Chen, Q. Niu, and A. H. MacDonald, Phys. Rev. Lett. 112, 017205 (2014).
  * (31) C. Sürgers, G. Fischer, P. Winkel, and H. v. Löhneysen, Nat. Commun. 5, 3400 (2014).
  * (32) S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212 (2015).
  * (33) T. Suzuki, R. Chisnell, A. Devarakonda, Y.-T. Liu, W. Feng, D. Xiao, J. W. Lynn, and J. G. Checkelsky, Nat. Phys. 12, 1119 (2016).
  * (34) L. Smejkal, R. Gonzalez-Hernandez, T. Jungwirth, and J. Sinova, Sci. Adv. 6, eaaz8809 (2020).
  * (35) L. Smejkal, A.H. MacDonald, J. Sinova, S. Nakatsuji, T. Jungwirth, Nat. Rev. Mater. 7, 482 (2022).
  * (36) L. Fu, Phys. Rev. Lett. 103, 266801 (2009). Under $\mathcal{T}$-symmetry, we have $k_{+}\rightarrow k_{-}$ and $\sigma_{+}\rightarrow\sigma_{-}$ with $k_{\pm}=k_{x}\pm ik_{y}$ and $\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$; under $N$-fold rotation symmetry, we have $k_{\pm}\rightarrow e^{\pm i2\pi/N}k_{\pm}$ and $\sigma_{\pm}\rightarrow e^{\pm i2\pi/N}\sigma_{\pm}$.
  * (37) See section III in the supporting materials for the details.
  * (38) We have $T_{yxx}=k_{x}mv^{2}/(4h^{5})$ and $T_{xxx}=-3k_{y}mv^{2}/(8h^{5})$ from which we obtain $g_{1}=-15k_{x}^{4}+40k_{x}^{2}k_{y}^{2}-25k_{y}^{4}$ and $g_{2}=-5k_{x}^{2}-10k_{y}^{2}$. Since $g_{2}$ is negative definite, its contribution to Eq.(16) is nonzero and hence the third-order IAHE is nonzero for this model.
  * (39) D. Kriegner, H. Reichlova, J. Grenzer, W. Schmidt, E. Ressouche, J. Godinho, T. Wagner, S. Y. Martin, A. B. Shick, V. V. Volobuev, G. Springholz, V. Holý, J. Wunderlich, T. Jungwirth, and K. Výborný, Phys. Rev. B 96, 214418 (2017).
  * (40) N.J. Ghimire, A.S. Botana, J.S. Jiang, J. Zhang, Y.-S. Chen, and J. F. Mitchell, Nat. Commun. 9, 3280 (2018).