# Probabilistic Stability of Stochastic Gradient Descent

Liu Ziyin1 Botao Li2 Tomer Galanti3 Masahito Ueda1,4
1Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan
2Laboratoire de Probabilités, Statistique et Modélisation/Université Paris
Cité, Paris, France
3 Center of Brains, Minds and Machines (CBMM), Massachusetts Institute of
Technology, Cambridge, MA, USA
4RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama 351-0198,
Japan

###### Abstract

Characterizing and understanding the stability of Stochastic Gradient Descent
(SGD) remains an open problem in deep learning. A common method is to utilize
the convergence of statistical moments, esp. the variance, of the parameters
to quantify the stability. We revisit the definition of stability for SGD and
propose using the convergence in probability condition to define the
probabilistic stability of SGD. The probabilistic stability sheds light on a
fundamental question in deep learning theory: how SGD selects a meaningful
solution for a neural network from an enormous number of possible solutions
that may severely overfit. We show that only through the lens of probabilistic
stability does SGD exhibit rich and practically relevant phases of learning,
such as the phases of the complete loss of stability, incorrect learning where
the model captures incorrect data correlation, convergence to low-rank
saddles, and correct learning where the model captures the correct
correlation. These phase boundaries are precisely quantified by the Lyapunov
exponents of the dynamics. The obtained phase diagrams imply that SGD prefers
low-rank saddles in a neural network when the underlying gradient is noisy,
thereby influencing the learning performance.

##  1 Introduction

Stochastic gradient descent (SGD) is the primary workhorse for optimizing
neural networks. As such, an important problem in deep learning theory is to
characterize and understand how SGD selects the solution of a deep learning
model, which often exhibits remarkable generalization capability. At the heart
of this problem lies the stability of SGD because the models trained with SGD
stay close to the attractive solutions where the dynamics is stable and moves
away from unstable ones. Solving this problem thus hinges on having a good
definition of the stability of SGD. The stability of SGD is often defined as a
function of the variance of the model’s parameters or gradients during
training. The hidden assumption behind this mainstream idea is that when the
variance diverges, the training becomes unstable (Wu et al., 2018; Zhu et al.,
2018; Liu et al., 2020, 2021; Ziyin et al., 2022b). In some sense, the idea
that the variance of the parameters matters the most is also an underlying
assumption in the deep learning optimization literature, where the utmost
important quantity is how fast the variance and the expected distance of the
parameters decay to zero (Vaswani et al., 2019; Gower et al., 2019). We
revisit this perspective and show that a variance-based notion of stability is
insufficient to understand the empirically observed stability of training of
SGD. In fact, we demonstrate a lot of natural learning settings where the
variance of SGD diverges, yet the model still converges with high probability.

In this work, we study the convergence in probability condition to understand
the stability of SGD. We then show that this stability condition can be
quantified with the Lyapunov exponent (Lyapunov, 1992) of the optimization
dynamics of SGD, a quantity deeply rooted in the study of dynamical systems
and has been well understood in physics and control theory (Eckmann & Ruelle,
1985). The main contribution of this work is to propose a new notion of
stability that sheds light on how SGD selects solutions and to discover
multiple deep-learning phenomena that can only be understood in terms of this
notion. Perhaps the most important implication of our theory is the
characterization of the highly nontrivial and practically important phase
diagram of SGD in a neural network loss landscape.

##  2 Probabilistic Stability

In this section, we introduce probabilistic stability, a rather general
concept that appears in many scenarios in control theory (Khas’ minskii, 1967;
Eckmann & Ruelle, 1985; Teel et al., 2014). The definition of probabilistic
stability relies on the notion of convergence in probability. We use
$\|\cdot\|_{n}$ to denote the $n-$norm of a matrix or vector and $\|\cdot\|$
to denote the case where $n=2$. The notation $\to_{p}$ indicates convergence
in probability.

######  Definition 1.

A sequence of random variables $\\{\theta_{t}\\}_{t}^{\infty}$ is
probabilistically stable at $\theta^{*}$ if $\theta_{t}\to_{p}\theta^{*}$.

Here, the notation $\to_{p}$ denotes convergence in probability. A sequence
$z_{t}$ converges in probability to $\theta^{*}$ if
$\lim_{t\to\infty}\mathbb{P}(\|\theta_{t}-\theta^{*}\|<\epsilon)=0$ for any
$\epsilon>0$.

We will see that this notion of stability is especially suitable for studying
the stability of the stationary points in SGD, be it a local minimum or a
saddle point. In contrast, the popular type of stability is based on the
convergence of statistic moments, which we also define below.

######  Definition 2.

A sequence of random variables $\\{\theta_{t}\\}_{t}^{\infty}$ is $L_{p}$-norm
stable at $\theta^{*}$ if
$\lim_{t\to\infty}\mathbb{E}\|\theta_{t}-\theta^{*}\|_{p}^{p}\to 0$.

For deep learning, the sequence of $\theta_{t}$ is the model parameters
obtained by the iterations of the SGD algorithm for a neural network. For
dynamical systems in general and deep learning specifically, it is impossible
to analyze the convergence behavior of the dynamics starting from an arbitrary
initial condition. Therefore, we have to restrict ourselves to the
neighborhood of a given stationary point and consider the linearized dynamics
around it.

The type of dynamics we consider in this work are all of the following form:

$\theta_{t+1}=\theta_{t}-\lambda\hat{H}(x)(\theta_{t}-\theta^{*})+O(||\theta_{t}-\theta^{*}||^{2}),$
(1)

where $\lambda$ is the learning rate, $\theta^{*}$ is the critical point under
consideration, $\hat{H}(x)$ is a random symmetric matrix that is a function of
the random variables $x$, which can stand for both a single data point or a
minibatch of data points. In this work, $\theta^{*}$ is said to be a local
minimum if $H:=\mathbb{E}_{x}[\hat{H}(x)]$ is positive semi-definite (PSD) and
is a saddle point otherwise.

Note that $\hat{H}(x)$ does not have to be tied to a single data point but can
also be the Hessian of a minibatch of data points. Essentially, it is the data
distribution of $\hat{H}$ that matters. When we have the same data
distribution but a different batch size $S$, the distribution of $\hat{H}$ is
different. This linearized dynamics is especially suitable for studying two
types of stationary points that appear in modern deep learning: (1)
interpolation minima and (2) symmetry-induced saddle points. Let us first
consider the stability of the interpolation minimum.

A minimum is said to be an interpolation minimum if the loss function reaches
zero for all data points. This means that close to such a minimum, the per-
sample loss functions all have vanishing first-order derivatives and positive-
semidefinite Hessians (Wu et al., 2018):

$\ell(\theta,x)=(\theta-\theta^{*})^{T}\hat{H}(x)(\theta-\theta^{*})+O(||\theta||^{3}).$
(2)

The dynamics of $w$ thus obeys Eq. (1). This type of minimum is of great
importance in deep learning because modern neural networks are often
“overparametrized,” and overparametrized networks under gradient flow are
observed to reach these interpolation minima easily.

The dynamics we consider is more general than that around an interpolation
minimum because the Hessians $\hat{H}$ in Eq. (1) are allowed to have
eigenvalues of both positive and negative signs, whereas Eq. (2) only allows
for positive semidefinite $\hat{H}$. Therefore, the general solution of Eq.
(1) also helps us understand Eq. (2) once we restrict the study to PSD
Hessians. The types of saddle points that obey Eq.(1) have been found to
widely exist when there is any type of loss function symmetry, such as
permutation symmetry (Fukumizu & Amari, 2000; Simsek et al., 2021; Entezari et
al., 2021; Hou et al., 2019), rescaling symmetry (Dinh et al., 2017; Neyshabur
et al., 2014), and rotation symmetry (Ziyin et al., 2023b). Recently, it has
been shown that every symmetry in the loss function leads to a critical point
of this type, and, more importantly, the subset of parameters relevant to the
symmetry, up to the leading order, has a dynamics completely independent of
the rest of the parameters (Ziyin, 2023). This type of critical point is
mostly likely a saddle point, but can easily become local minima or global
minima when weight decay is used. This means that Eq. (1) can be seen as an
effective description for only a small subset of all parameters in the model
under consideration and is not as restrictive as it naively seems to be.

##  3 The Probabilistic Stability of SGD

This section presents the main theoretical results. We first show that the
dynamics is exactly solvable for a rank-$1$ dynamics. We then prove a result
showing that no moment-based stability can understand the stability of SGD
around the saddle. Lastly, we prove that the probabilistic stability of SGD in
high-dimension is equivalent to a condition on the sign of the Lyapunov
exponent of the SGD dynamics.

###  3.1 Rank-1 Dynamics

Let us first consider the case in which $\hat{H}(x)=h(x)nn^{T}$ is rank-$1$
for a random scalar function $h(x)$, and a fixed unit vector $n$ for all data
points $x$. Thus, the dynamics simplifies to a one-dimensional dynamics, where
$h(x)\in\mathbb{R}$ is the corresponding eigenvalue of $\hat{H}(x)$:

$\theta_{t+1}=\theta_{t}-\lambda h(x)(\theta_{t}-\theta^{*}).$ (3)

######  Theorem 1.

Let $\theta_{t}$ follow Eq. (3). Then, for any distribution of $h(x)$,

$n^{T}(\theta_{t}-\theta^{*})\to_{p}0$ (4)

if and only if111This condition generalizes to the case when the batch size
$S$ is larger than $1$, where $h(x)$ becomes the per-batch Hessian, and the
expectation is taken over all possible batches.

$\mathbb{E}_{x}[\log|1-\lambda h(x)|]<0.$ (5)

Figure 1: SGD exhibits a complex phase diagram through the lens of
probabilistic stability. Left: For a matrix factorization saddle point, the
dynamics of SGD can be categorized into at least five different phases. Phase
I, II, and IV correspond to a successful escape from the saddle. Phase III is
where the model converges to a low-rank saddle point. Phase I corresponds to
the case $w_{t}\to_{p}u_{t}$, which signals correct learning. In phase Ib, the
model also converges in variance. Phase II corresponds to stable but incorrect
learning, where $w_{t}\to_{p}-u_{t}$. Phase IV corresponds to complete
instability. Right: the phases of SGD can quantified by the sign of the
Lyapunov exponent $\Lambda$. Where $\Lambda<0$, SGD collapses to a saddle
point; when $\Lambda>0$, SGD escapes the saddle and enters a escaping phase.
The two escaping phases are qualitatively different. For a small learning
rate, the model is in a learning phase due to the repulsiveness of the saddle
point at a small learning rate, and the model is likely to converge to local
minima close to the saddle. For a very large learning rate, SGD escapes the
saddle due to the dynamical instability of SGD, and the model will move far
away from the saddle. Besides, the magnitude of the Lyapunov exponent can also
quantity the speed of the learning dynamics. See Appendix A.1 for numerical
details of this example.

The condition (5) is a sharp characterization of when a stationary condition
or an interpolation minimum becomes attractive. It also works with weight
decay. When weight decay is present, the diagonal terms of $\hat{H}$ are
shifted by $\gamma$, and so $h=h^{\prime}+\gamma$.

At what learning rate is the condition violated? To leading orders in
$\lambda$, this can be identified by expanding the logarithmic term up to the
second order in $\lambda$:

$\mathbb{E}_{x}[\log|1-\lambda
h(x)|]=-\lambda\mathbb{E}[h(x)]-\frac{1}{2}\lambda^{2}\mathbb{E}_{x}[h(x)^{2}]+(\lambda^{3}).$
(6)

Ignoring the second-order term, we see that the dynamics always follow the
sign of $\mathbb{E}[h(x)]$, in agreement with the GD algorithm. If
$\mathbb{E}[h]>0$, the condition is always violated. When the second-order
term is taken into consideration, the fluctuation of $h(x)$ now decides the
stability of the stationary condition. The stationary condition is attractive
if

$\lambda>2\frac{-\mathbb{E}[h(x)]}{\mathbb{E}[h(x)^{2}]}.$ (7)

This result implies that the stationary condition can be attractive even if
$h<0$. Comparing this critical learning rate with the continuous-time analysis
in Ziyin et al. (2023a), we see that the two results agree. The r.h.s. of the
condition also has a natural interpretation as a signal-to-noise ratio in the
gradient. The denominator is the Hessian of the original loss function, which
determines the signal in the gradient. The denominator is the strength of the
gradient noise in the minibatch, and has been termed as “gradient non-
uniformity” in previous literature Wu et al. (2018). An illustration of this
solution is given in Figure 1. We show the probabilistic stability conditions
for a rank-$1$ saddle point with a rescaling symmetry (see Section 4). The
loss function is $\ell(u,w)=-xyuw+o(u^{2}+w^{2})$. Here, the data points
$xy=1$, and $xy=a$ are sampled with equal probability. These saddles appear
naturally in matrix factorization problems and also in recent sparse learning
algorithms (Poon & Peyré, 2021, 2022; Ziyin & Wang, 2023; Kolb et al., 2023).

When the dynamics is high-dimensional, the problem becomes harder to solve
because for each realization of $x$, $\hat{H}(x)$ do not commute with each
other. While an analytical solution to the stability condition in high-
dimension is unlikely to exist, we can say something quite general about them.

###  3.2 Insufficiency of Norm-Stability

Probabilistic stability offers a new tool for understanding the stability of
the SGD algorithm around a given solution. Theorem 1 provides a perfect
example to compare the probabilistic stability with the norm-stability, which
is a conventional way to study the stability of SGD (Wu et al., 2018; Zhu et
al., 2018; Liu et al., 2020, 2021; Ziyin et al., 2022b; Vaswani et al., 2019;
Gower et al., 2019). The following rather trivial proposition shows that if
SGD converges to a point in $L_{p}$-norm, it must converge in probability.

######  Proposition 1.

If $\theta_{t}$ is stable at $\theta^{*}$ in $L_{p}$ norm, then it is stable
at $\theta^{*}$ in probability.

The proof follows from the well-known fact that convergence in $L_{p}$ norm
implies the convergence in probability. This implies that norm stability is a
more restricted notion than probabilistic stability. In many cases, the two
types of stabilities agree. However, we will see that for SGD, the two types
of stability conditions can offer dramatically different predictions, which is
constructively established by the following proposition.

######  Proposition 2.

Let $\theta_{t}$ follow Eq. (1) around a critical point $\theta^{*}$. Then,

  1. 1.

for any fixed $\lambda$, there exists a data distribution such that
$\theta_{t}$ is probabilistically stable but $L_{p}$-stable;

  2. 2.

if $\theta^{*}$ is a saddle point and $p\geq 1$, the set of $\theta_{0}$ that
is $L_{p}$-stable has Lebesgue measure zero.

Therefore, this means that the $L_{p}$-stability is not useful in
understanding the stability of SGD close to saddle points. One reason is that
the outliers strongly influence the $L_{p}$ norm in the data, whereas the
probabilistic stability is robust against such outliers.

###  3.3 Lyapunov Exponent and Probabilistic Stability

Extending the probabilistic stability to high-rank dynamics is nontrivial
because the stability of SGD is generally initialization dependent, unlike the
rank-$1$ case, where the condition is found to take an exact form and is
initialization independent. It is thus useful to consider the worst-case
initialization. Here, the crucial quantity is the Lyapunov exponent of a point
$\theta^{*}$:222More appropriately, this should be called the maximum Lyapunov
exponent, which is initialization-independent. One can also consider the
initialization-dependent Lyapunov exponent. If $\bar{H}$ is a $d$-by-$d$
matrix, a well-known fact is that the initialization-dependent Lyapunov
exponent takes at most $d$ distinctive values. Conceptually, this means that
there can be, at most, $d$ collapses at different critical learning rates.

$\Lambda=\max_{\theta_{0}}\lim_{t\to\infty}\mathbb{E}\left[\log\frac{\|\theta_{t}-\theta^{*}\|}{\|\theta_{0}\|}\right].$
(8)

Here, the expectation is taken over the random samplings of the SGD algorithm.
In general, $\Lambda$ does not vanish. The following theorem shows that SGD is
probabilistically stable at a point if and only if its Lyapunov exponent is
negative.

######  Theorem 2.

Assuming that $\Lambda\neq 0$, the linearized dynamics of SGD is
probabilistically stable at $\theta^{*}$ for any $\theta_{0}$ if and only if
$\Lambda<0$.

The proof of this theorem invokes the Furstenberg-Kesten theorem (Furstenberg
& Kesten, 1960). It is straightforward to show that the Lyapunov exponent
exists and to find an upper and lower bound. For a finite-size dataset, we can
define $r_{\max}$ to be larger than the absolute value of the eigenvalues of
$I-\lambda\hat{H}(x)$ for all $x$. Similarly, we can define $r_{\min}>0$ to be
smaller than all the absolute values of all the eigenvalues of
$I-\lambda\hat{H}(x)$ for all $x$. Therefore, it is easy to check that

$\log r_{\min}<\Lambda<\log r_{\max}.$ (9)

However, it is difficult to give a better estimation of the exponent. In fact,
it is a well-known open problem in the field of dynamical systems to find an
analytical expression of the Lyapunov exponent (Crisanti et al., 2012;
Pollicott, 2010; Jurga & Morris, 2019).

Now, we give two quantitative estimates about when the Lyapunov exponent will
be negative. This discussion also implies a sufficient but weak condition for
a general type of multidimensional dynamics to converge in probability. Let
$h^{*}(x)$ be the largest eigenvalue of $\hat{H}(x)$ and assume that
$1-h^{*}(x)>0$ for all $x$. Then, the following condition implies that
$w\to_{p}0$:

$\mathbb{E}_{x}[\log|1-\lambda h^{*}(x)|]<0$ (10)

which mimics the condition we found for rank-$1$ systems. An alternative
estimate can be made by assuming that $\hat{H}(x)$ commute with
$\hat{H}(x^{\prime})$ for all $x$ and $x^{\prime}$. If $\hat{H}$ has rank $d$,
This reduces the problem to $d$ separated rank-1 dynamics, and Theorem 1 gives
the exact solutions in each subspace. Numerical evidence shows that the
commutation approximation quite accurately predicts the onset of low-rank
behavior the actual rank (see Section 4 and Appendix A.4).

Another relevant question is whether this theorem is trivial for SGD at a high
dimension in the sense that it could be the case that $\Lambda$ could be
identically zero independent of the dataset. One can show that the Lyapunov
exponent is generally nonzero for all datasets that satisfy a mild condition.
Let $\mathbb{E}[\hat{H}]$ be full rank. By definition,

$\Lambda=\lim_{t\to\infty}\frac{1}{t}\mathbb{E}\left[\log\theta_{0}^{T}\left(\prod_{j}^{t}(I-\lambda\hat{H}_{i_{j}})\right)\left(\prod_{j}^{t}(I-\lambda\hat{H}_{i_{j}})\right)^{T}\theta_{0}\right]=-\frac{2\lambda\theta_{0}^{T}H\theta_{0}}{\|\theta_{0}\|^{2}}+O(\lambda^{2}).$
(11)

Therefore, as long as $\lambda$ is sufficiently small, the sign of the
Lyapunov exponent is opposite to the sign of the eigenvalues of $H$. This
proves something quite general for SGD at an interpolation minimum: with a
small learning rate, the model converges to the minimum exponentially fast, in
agreement with common analysis in the optimization literature. See Figure
1-right for numerical computation of Lyapunov exponents of a matrix
factorization problem and the corresponding phases.

##  4 Phases of Learning

With this notation of stability, we can study the actual effect of minibatch
noise on a neural network-like landscape. The theory has interesting
implications for the stability of interpolation minima in deep learning, which
we explore in Appendix C.1. In the main text, we focus on the stability of SGD
on saddle points.

A commonly studied minimal model of the landscape of neural networks is a deep
linear net (or deep matrix factorization) (Kawaguchi, 2016; Lu & Kawaguchi,
2017; Ziyin et al., 2022a; Wang & Ziyin, 2022). For these problems, we
understand that all local minima are identical copies of each other, and so
all local minima have the same generalization capability (Kawaguchi, 2016; Ge
et al., 2016). A deep linear net’s special and interesting solutions are the
saddle points, which are low-rank solutions and often achieve similar training
loss with dramatically different generalization performances.333The influence
of these singularities on generalization has been established by the singular
learning theory (Watanabe & Opper, 2010). More importantly, these saddles
points also appear in nonlinear models with similar geometric properties, and
they could be a rather general feature of the deep learning landscape (Brea et
al., 2019). It is thus important to understand how the noise of SGD affects
the stability of a low-rank saddle here. Let the loss function be
$\mathbb{E}_{x}[(\sum_{i}u_{i}\sigma(w_{i}x)-y)^{2}/2]$, where
$\sigma(x)=c_{0}x+O(x^{2})$ is any nonlinearity that is locally linear at
$x=0$. We let $c_{0}=1$ and focus on cases where both $x$ and $y$ are one-
dimensional. Locally around $u,\ w\approx 0$, the model $u^{T}w$ is either
rank-$1$ or rank-$0$. The rank-$0$ point where $u_{i}=w_{i}=0$ for all $i$ is
a saddle point as long as $\mathbb{E}[xy]\neq 0$. In this section, we show
that the stability of this saddle point features complex and dramatic phase
transition-like behaviors as we change the learning rate of SGD.

Consider the linearized dynamics around the saddle at $w_{i}=u_{i}=0$. The
expanded loss function takes the following form:

$\ell(u,w)=-{xy}\sum^{d}_{i}u_{i}w_{i}+const.$ (12)

For learning to happen, SGD needs to escape from the saddle point. Let us
consider a simple data distribution where $xy=1$ and $xy=a$ with equal
probability. When $a>-1$, correct learning happens when ${\rm sign}(w)={\rm
sign}(u)$. We thus focus on the case of $a>-1$. The case of $a<-1$ is
symmetric to this case up to a rescaling. This example is already presented in
Figure 1. The two most important observations are: (1) SGD can indeed converge
to low-rank saddle points; however, this happens only when the gradient noise
is sufficiently strong and when the learning rate is large (but not too
large); (2) the region for convergence to saddles (region III) is exclusive
with the region for convergence in mean square (Ia), and thus one can only
understand the saddle-seeking behavior of SGD within the proposed
probabilistic framework. We prove the following proposition. It becomes
evident that the low-rank solution is reached when both $w_{t}+u_{t}\to_{p}0$
and $w_{t}-u_{t}\to_{p}0$.

######  Proposition 3.

For any $w_{0},u_{0}\in\mathbb{R}/\\{0\\}$. $w_{t}-u_{t}\to_{p}0$ if and only
if $\mathbb{E}_{x}[\log|1-\lambda xy|]<0$. $w_{t}+u_{t}$ converges to $0$ in
probability if and only if $\mathbb{E}_{x}[\log|1+\lambda xy|]<0$.

Figure 2: Phase diagrams of SGD stability. The definitions of the phases are
the same as those in Figure 1. We sample a dataset of size $N$ such that
$x\sim\mathcal{N}(0,1)$ and noise $\epsilon\sim\mathcal{N}(0,4)$, and generate
a noisy label $y=\mu x+(1-\mu)\epsilon$. Left: the $\lambda-\mu$ phase diagram
for $S=1$ and $N=\infty$. Right: The $\lambda-S$ phase diagram for $\mu=0.06$
and $N=\infty$.

Figure 3: Convergence to low-rank solutions in nonlinear neural networks. At
every training step, we sample input $x\sim\mathcal{N}(0,I_{200})$ and noise
$\epsilon\sim\mathcal{N}(0,\sqrt{2}I_{200})$, and generate a noisy label
$y=\mu x+(1-\mu)\epsilon$, where $1-\mu$ controls the level of the noise.
Left: Linear network. Right: tanh network. The white dashed line shows the
theoretical prediction of the appearance of low-rank structure computed by
numerically integrating the condition in Proposition 3.

Figure 4: Density ($1-sparsity$) of the convolutional layers in a ResNet18,
when there is static noise (mislabeling) in the training data. Here, we show
the number of sparse parameters in the two of the largest convolutional
layers, each containing roughly one million parameters in total. The figures
respectively show layer1.1.conv2 (upper left), layer2.1.conv2 (upper right),
layer3.1.conv2 (lower left), and layer4.1.conv2 (lower right).

Figure 5: Rank (left) and test accuracy (right) of the ResNet18 trained in a
data set with static noise. The transition of rank has a clear boundary. The
model has a full rank but random-guess level performance for large noise and
small learning rates. Here, noise refers to the probability that data point is
mislabeled.

The theory shows that the phase diagram of SGD strongly depends on the data
distribution, and it is interesting to explore and compare a few different
settings. Now, we consider a size-$N$ Gaussian dataset. Let
$x_{i}\sim\mathcal{N}(0,1)$ and noise $\epsilon_{i}\sim\mathcal{N}(0,4)$, and
generate a noisy label $y_{i}=\mu x_{i}+(1-\mu)\epsilon_{i}$. See the phase
diagram for this dataset in Figure 2 for an infinite $N$. The phase diagrams
in Figure 7 show the phase diagram for a finite $N$. We see that the phase
diagram has a very rich structure at a finite size. We make three rather
surprising observations about the phase diagrams: (1) as $N\to\infty$, the
phase diagram becomes smoother and smoother and each phase takes a connected
region (cf. finite size experiments in Appendix A.2); (2) phase II seems to
disappear as $N$ becomes large; (3) the lower part of the phase diagram seems
universal, taking the same shape for all samplings of the datasets and across
different sizes of the dataset. This suggests that the convergence to low-rank
structures can be a universal aspect of SGD dynamics, which corroborates the
widely observed phenomenon of collapse in deep learning (Papyan et al., 2020;
Wang & Ziyin, 2022; Tian, 2022). The theory also shows that if we fix the
learning rate and noise level, increasing the batch size makes it more and
more difficult to converge to the low-rank solution (see Figure 8, for
example). This is expected because the larger the batch size, the smaller the
effective noise in the gradient.

Many recent works have suggested how neural networks could be biased toward
low-rank solutions. Theoretically, Galanti et al. (2023) showed that with weak
weight decay, SGD is biased towards low-rank solutions. Ziyin et al. (2022a)
showed that GD converges to a low-rank solution with weight decay. Therefore,
weight decay already induces a low-rank bias in learning, and it is unknown if
SGD alone has any bias toward low-rank solutions. Andriushchenko et al. (2022)
showed empirical hints of a preference for low-rank solutions when training
without SGD. However, it remains to be clarified when or why SGD has such a
preference on its own. To the best of our knowledge, our theory is the first
to precisely characterize the low-rank bias of SGD in a deep learning setting.
Compared with the stability diagram of linear regression, this result implies
that a large learning rate can both help and hinder optimization.

##### Phase Diagram for Neural Networks

There is a strong sense of universality in the lower-left part of the phase
diagram in Figure 1 since they all take a similar shape independent of the
size or sampling of the data points. We now verify its existence in actual
neural networks.

We start with a controlled experiment where, at every training step, we sample
input $x\sim\mathcal{N}(0,I_{200})$ and noise
$\epsilon\sim\mathcal{N}(0,4I_{200})$, and generate a noisy label $y=\mu
x+(1-\mu)\epsilon$. Note that $1-\mu$ controls the level of the noise.
Training proceeds with SGD on the MSE loss. We train a two-layer model with
the architecture: $200\to 200\to 200$. See Figure 3 for the theoretical phase
diagram, which is estimated under the commutation approximation. Under SGD,
the model escapes from the saddle with a finite variance to the right of the
dashed line and has an infinite variance to its left. In the region
$\lambda\in(0,0.2)$, this loss of the $L_{2}$ stability condition coincides
with the condition for the convergence to the saddle. The experiment shows
that the theoretical boundary agrees well with the numerical
results.444Commonly used variants of SGD such as the Adam optimizer (Kingma &
Ba, 2014) also have a qualitatively similar phase diagram. See Appendix A This
suggests that the effects we studied are rather universal, not just a special
feature of SGD.

Lastly, we train independently initialized ResNets on CIFAR-10 with SGD. The
training proceeds with SGD without momentum at a fixed learning rate and batch
size $S=32$ (unless specified otherwise) for $10^{5}$ iterations. Our
implementation of Resnet18 contains $11$M parameters and achieves $94\%$ test
accuracy under the standard training protocol, consistent with the established
values. To probe the effect of noise, we artificially inject a dynamical label
noise during every training step, where, at every step, a correct label is
flipped to a random label with probability $noise$, and we note that the phase
diagrams are similar regardless of whether the noise is dynamical or static
(where the mislabelling is fixed). See Figure 5 for the phase diagram of
static label noise. Interestingly, the best generalization performance is
achieved close to the phase boundary when the noise is strong. This is direct
evidence that SGD noise has a strong regularization effect on the trained
model. We see that the results agree with the theoretical expectation and the
analytical model’s phase diagram. We also study the sparsity of the ResNets in
different layers in Figure 5, and we observe that the phase diagrams are all
qualitatively similar. Also, see Appendix A for the experiment with a varying
batch size.

Figure 6: How SGD selects a solution. Left: The landscape of a two-layer
network with the swish activation function (Ramachandran et al., 2017).
Middle, Right: the generalization performance of the model for different
learning rates. Middle: Initialized at solution B, SGD first jumps to C and
then diverges. Right: Initialized at A, SGD also jumps to C and diverges. In
both cases, the behavior of SGD agrees with the prediction of the
probabilistic stability instead of the $L_{2}$ stability. Instead of jumping
between local minima, SGD, at a large learning rate, transitions from minima
to saddles.

##  5 Which Solution Does SGD Prefer?

We now investigate one of the most fundamental problems in deep learning
through the lens of probabilistic stability: how SGD selects a solution for a
neural network. In this section, we study a two-layer network with a single
hidden neuron with the swish activation function: $f(w,u,x)=u\times{\rm
swish}(wx)$, where ${\rm swish}(x)=x\times{\rm sigmoid}(x)$. Swish is a
differentiable ReLU variant discovered by meta-learning techniques and
consistently outperforms ReLU in various tasks. We generate $100$ data points
$(x,y)$ as $y=0.1{\rm swish}(x)+0.9\epsilon$, where both $x$ and $\epsilon$
are sampled from normal distributions. See Figure 6 for an illustration of the
training loss landscape. There are two local minima: solution A at roughly
$(-0.7,-0.2)$ and solution B at $(1.1,-0.3)$. Here, the solution with better
generalization is A because it captures the correct correlation between $x$
and $y$ when $x$ is small. Solution A is also the sharper one; its largest
Hessian eigenvalue is roughly $h_{a}=7.7$. Solution B is the worse solution;
it is also the flatter one, with the largest Hessian value being $h_{b}=3.0$.
There is also a saddle point C at $(0,0)$, which performs significantly better
than B and slightly worse than A in generalization.

If we initialize the model at A, $L_{2}$ stability theory would predict that
as we increase the learning rate, the model moves from the sharper solution A
to the flatter minimum B when SGD loses $L_{2}$ stability in A; the model
would then lose total stability once SGD becomes $L_{2}$-unstable at B. As
shown by the red arrows in Figure 6. In contrast, probabilistic stability
predicts that SGD will move from A to C as C becomes attractive and then lose
stability, as the black arrows indicate. See the right panel of the figure for
the comparison with the experiment for the model’s generalization performance.
The dashed lines show the predictions of the $L_{2}$ stability and
probabilistic theories, respectively. We see that the probabilistic theory
predicts both the error and the place of transition right, whereas $L_{2}$
stability neither predicts the right transition nor the correct level of
performance.

If we initialize at B, the flatter minimum, $L_{2}$ stability theory would
predict that the model will only have one jump from B to divergence as we
increase the learning rate. Thus, from $L_{2}$ stability, SGD would have
roughly the performance of B until it diverges, and having a large learning
rate will not help increase the performance. In sharp contrast, the
probabilistic stability predicts that the model will have two jumps: it stays
at B for a small $\lambda$ and jumps to C as it becomes attractive at an
intermediate learning rate. The model will ultimately diverge if C loses
stability. Thus, our theory predicts that the model will first have a bad
performance, then show a better performance at an intermediate learning rate,
and finally diverge. See the middle panel of Figure 6. We see that the
prediction of the probabilistic stability agrees with the experiment and
correctly explains why SGD leads to better performance.

##  6 Discussion

In this work, we have demonstrated that the convergence in probability
condition serves as an essential notion for understanding the stability of
SGD, leading to highly nontrivial and practically relevant phase diagrams at a
finite learning rate. Crucially, these effects are only present for SGD and
not for GD, demonstrating that the algorithmic regularization due to SGD is
qualitatively different from that of GD. We also clarified its intimate
connection to Lyapunov exponents, which are fundamental metrics of stability
in the study of dynamical systems. The proposed stability agrees with the
norm-based notion of stability at a small learning rate and large batch size.
At a large learning rate and a small batch size, we have shown that the
proposed notion of stability captures the actual behavior of SGD much better
and successfully explains a series of experiment phenomena that have been
quite puzzling. Among the many implications that we discussed, perhaps the
most fundamental one is a novel understanding of the implicit bias of SGD.
When viewed from a dynamical stability point of view, the implicit bias of
stochastic gradient descent is thus fundamentally different from the implicit
bias of gradient descent. The new perspective we provide is also different
from the popular perspective that SGD influences the learning outcome directly
by making the model converge to flat minima. In fact, in our construction, the
flatter minimum does not have better generalization properties, nor does SGD
favor it over the sharper one. Instead, SGD performs a selection between the
case of converging to saddles and local minima and that of directly and
significantly affecting its performance. Thus, we believe that probabilistic
stability and Lyapunov exponents are important future directions to
investigate the stability of SGD in a deep-learning scenario.

## References

  * Andriushchenko et al. (2022) Maksym Andriushchenko, Aditya Varre, Loucas Pillaud-Vivien, and Nicolas Flammarion.  Sgd with large step sizes learns sparse features.  _arXiv preprint arXiv:2210.05337_ , 2022.
  * Brea et al. (2019) Johanni Brea, Berfin Simsek, Bernd Illing, and Wulfram Gerstner.  Weight-space symmetry in deep networks gives rise to permutation saddles, connected by equal-loss valleys across the loss landscape.  _arXiv preprint arXiv:1907.02911_ , 2019.
  * Crisanti et al. (2012) Andrea Crisanti, Giovanni Paladin, and Angelo Vulpiani.  _Products of random matrices: in Statistical Physics_ , volume 104\.  Springer Science & Business Media, 2012.
  * Dinh et al. (2017) L. Dinh, R. Pascanu, S. Bengio, and Y. Bengio.  Sharp Minima Can Generalize For Deep Nets.  _ArXiv e-prints_ , March 2017.
  * Eckmann & Ruelle (1985) J-P Eckmann and David Ruelle.  Ergodic theory of chaos and strange attractors.  _The theory of chaotic attractors_ , pp. 273–312, 1985.
  * Entezari et al. (2021) Rahim Entezari, Hanie Sedghi, Olga Saukh, and Behnam Neyshabur.  The role of permutation invariance in linear mode connectivity of neural networks.  _arXiv preprint arXiv:2110.06296_ , 2021.
  * Fukumizu & Amari (2000) Kenji Fukumizu and Shun-ichi Amari.  Local minima and plateaus in hierarchical structures of multilayer perceptrons.  _Neural networks_ , 13(3):317–327, 2000.
  * Furstenberg & Kesten (1960) Harry Furstenberg and Harry Kesten.  Products of random matrices.  _The Annals of Mathematical Statistics_ , 31(2):457–469, 1960.
  * Galanti et al. (2023) Tomer Galanti, Zachary S. Siegel, Aparna Gupte, and Tomaso Poggio.  Sgd and weight decay provably induce a low-rank bias in neural networks, 2023.
  * Ge et al. (2016) Rong Ge, Jason D Lee, and Tengyu Ma.  Matrix completion has no spurious local minimum.  _Advances in neural information processing systems_ , 29, 2016.
  * Gower et al. (2019) Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtárik.  Sgd: General analysis and improved rates.  In _International conference on machine learning_ , pp. 5200–5209. PMLR, 2019.
  * Hou et al. (2019) Tianqi Hou, KY Michael Wong, and Haiping Huang.  Minimal model of permutation symmetry in unsupervised learning.  _Journal of Physics A: Mathematical and Theoretical_ , 52(41):414001, 2019.
  * Jurga & Morris (2019) Natalia Jurga and Ian Morris.  Effective estimates on the top lyapunov exponents for random matrix products.  _Nonlinearity_ , 32(11):4117, 2019.
  * Kawaguchi (2016) Kenji Kawaguchi.  Deep learning without poor local minima.  _Advances in Neural Information Processing Systems_ , 29:586–594, 2016.
  * Khas’ minskii (1967) Rafail Zalmanovich Khas’ minskii.  Necessary and sufficient conditions for the asymptotic stability of linear stochastic systems.  _Theory of Probability & Its Applications_, 12(1):144–147, 1967.
  * Kingma & Ba (2014) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _CoRR_ , abs/1412.6980, 2014.  URL http://dblp.uni-trier.de/db/journals/corr/corr1412.html#KingmaB14.
  * Kolb et al. (2023) Chris Kolb, Christian L Müller, Bernd Bischl, and David Rügamer.  Smoothing the edges: A general framework for smooth optimization in sparse regularization using hadamard overparametrization.  _arXiv preprint arXiv:2307.03571_ , 2023.
  * Li et al. (2021) Zhiyuan Li, Sadhika Malladi, and Sanjeev Arora.  On the validity of modeling sgd with stochastic differential equations (sdes).  _Advances in Neural Information Processing Systems_ , 34:12712–12725, 2021.
  * Liu et al. (2021) Kangqiao Liu, Liu Ziyin, and Masahito Ueda.  Noise and fluctuation of finite learning rate stochastic gradient descent, 2021.
  * Liu et al. (2020) Yanli Liu, Yuan Gao, and Wotao Yin.  An improved analysis of stochastic gradient descent with momentum.  _Advances in Neural Information Processing Systems_ , 33:18261–18271, 2020.
  * Lu & Kawaguchi (2017) Haihao Lu and Kenji Kawaguchi.  Depth creates no bad local minima.  _arXiv preprint arXiv:1702.08580_ , 2017.
  * Lyapunov (1992) Aleksandr Mikhailovich Lyapunov.  The general problem of the stability of motion.  _International journal of control_ , 55(3):531–534, 1992.
  * Neyshabur et al. (2014) Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  In search of the real inductive bias: On the role of implicit regularization in deep learning.  _arXiv preprint arXiv:1412.6614_ , 2014.
  * Papyan et al. (2020) Vardan Papyan, XY Han, and David L Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  _Proceedings of the National Academy of Sciences_ , 117(40):24652–24663, 2020.
  * Pollicott (2010) Mark Pollicott.  Maximal lyapunov exponents for random matrix products.  _Inventiones mathematicae_ , 181(1):209–226, 2010\.
  * Poon & Peyré (2021) Clarice Poon and Gabriel Peyré.  Smooth bilevel programming for sparse regularization.  _Advances in Neural Information Processing Systems_ , 34:1543–1555, 2021.
  * Poon & Peyré (2022) Clarice Poon and Gabriel Peyré.  Smooth over-parameterized solvers for non-smooth structured optimization.  _arXiv preprint arXiv:2205.01385_ , 2022.
  * Ramachandran et al. (2017) Prajit Ramachandran, Barret Zoph, and Quoc V. Le.  Searching for activation functions, 2017.
  * Simsek et al. (2021) Berfin Simsek, François Ged, Arthur Jacot, Francesco Spadaro, Clément Hongler, Wulfram Gerstner, and Johanni Brea.  Geometry of the loss landscape in overparameterized neural networks: Symmetries and invariances.  In _International Conference on Machine Learning_ , pp. 9722–9732. PMLR, 2021.
  * Teel et al. (2014) Andrew R Teel, João P Hespanha, and Anantharaman Subbaraman.  A converse lyapunov theorem and robustness for asymptotic stability in probability.  _IEEE Transactions on Automatic Control_ , 59(9):2426–2441, 2014.
  * Tian (2022) Yuandong Tian.  Deep contrastive learning is provably (almost) principal component analysis.  _arXiv preprint arXiv:2201.12680_ , 2022.
  * Vaswani et al. (2019) Sharan Vaswani, Francis Bach, and Mark Schmidt.  Fast and faster convergence of sgd for over-parameterized models and an accelerated perceptron.  In _The 22nd international conference on artificial intelligence and statistics_ , pp. 1195–1204. PMLR, 2019.
  * Wang & Ziyin (2022) Zihao Wang and Liu Ziyin.  Posterior collapse of a linear latent variable model.  _arXiv preprint arXiv:2205.04009_ , 2022.
  * Watanabe & Opper (2010) Sumio Watanabe and Manfred Opper.  Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.  _Journal of machine learning research_ , 11(12), 2010.
  * Wu et al. (2018) Lei Wu, Chao Ma, et al.  How sgd selects the global minima in over-parameterized learning: A dynamical stability perspective.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Zhu et al. (2018) Zhanxing Zhu, Jingfeng Wu, Bing Yu, Lei Wu, and Jinwen Ma.  The anisotropic noise in stochastic gradient descent: Its behavior of escaping from sharp minima and regularization effects.  _arXiv preprint arXiv:1803.00195_ , 2018.
  * Ziyin (2023) Liu Ziyin.  Symmetry leads to structured constraint of learning, 2023.
  * Ziyin & Wang (2023) Liu Ziyin and Zihao Wang.  spred: Solving L1 Penalty with SGD.  In _International Conference on Machine Learning_ , 2023.
  * Ziyin et al. (2022a) Liu Ziyin, Botao Li, and Xiangming Meng.  Exact solutions of a deep linear network.  _arXiv preprint arXiv:2202.04777_ , 2022a.
  * Ziyin et al. (2022b) Liu Ziyin, Kangqiao Liu, Takashi Mori, and Masahito Ueda.  Strength of minibatch noise in SGD.  In _International Conference on Learning Representations_ , 2022b.  URL https://openreview.net/forum?id=uorVGbWV5sw.
  * Ziyin et al. (2023a) Liu Ziyin, Hongchao Li, and Masahito Ueda.  Law of balance and stationary distribution of stochastic gradient descent.  _arXiv preprint arXiv:2308.06671_ , 2023a.
  * Ziyin et al. (2023b) Liu Ziyin, Ekdeep Singh Lubana, Masahito Ueda, and Hidenori Tanaka.  What shapes the loss landscape of self supervised learning?  In _The Eleventh International Conference on Learning Representations_ , 2023b.  URL https://openreview.net/forum?id=3zSn48RUO8M.

##  Appendix A Additional Numerical Results

###  A.1 Experimental Detail for Figure 1-right

The experiment is performed for a two-dimensional system whose dynamics is
specified in (1). The expectation of the Hessian $\mathbb{E}[\hat{H}]$ is
chosen to be $\text{diag}(0.1,-0.1)$, while the noise is generated via a
normal $2\times 2$ random matrix $M_{\text{noise}}$. The noisy Hessian is
obtained as

$\hat{H}=\mathbb{E}[\hat{H}]+M_{\text{noise}}+M_{\text{noise}}^{T},$ (13)

and one can verify that such $\hat{H}$ is symmetric and consistent with our
choice of $\mathbb{E}[\hat{H}]$. The initial state is sampled from a unit
circle. The dynamics stops at time step $t$, and the Lyapunov exponent is
calculated as $\frac{1}{t}\log||\theta_{t}||$, if one of the three following
conditions is satisfied: $||\theta||$ reaches the upper cutoff of $10^{100}$;
$||\theta||$ reaches the lower cutoff of $10^{-140}$; the preset maximal
number of steps of $5000$ is reached. For each learning rate, the Lyapunov is
obtained as the average of the results collected in $800$ independent runs.

###  A.2 Phases of Finite-Size Datasets

See Figure 7.

Figure 7: Phase diagrams of SGD stability for finite-size dataset. The data
sampling is the same as in Figure 2. From upper left to lower right: $N=3,\
4,\ 8,\ 10,\ 24,\ 100$. As the dataset size tends to infinity, the phase
diagram converges to that in Figure 2. The lower parts of all the phase
diagrams look similar, suggesting a universal structure.

###  A.3 Effect of Changing Batchsize on Resnet18

See Figure 8.

Figure 8: The rank of the penultimate-layer representation of Resnet18 trained
with different levels of batch sizes. In agreement with the phase diagram, the
model escapes from the low-rank saddle as one increases the batch size.

###  A.4 Critical learning rate

Here, we compare the empirical rank of the solution with the commutation
approximated critical learning rates obtained in (7). See Figure 9. The
experiment is run on a two-layer fully connected linear network: $50\to 50\to
50$, which is equivalent to a matrix factorization problem. The model is
initialized with the standard Kaiming init. The dataset we consider is one
with a sparse but full-rank signal.

Let $\odot$ denote the Hadamard product. The input data is generated as
$x=m\odot X$, where $X\sim\mathcal{N}(0,I_{50})$ and $m$ is a random mask
where a random element is set to be $1$, and the rest is zero. The labels $Y$
is generated as $Y=\mu x+(1-\mu)(m\odot\epsilon)$, where
$\epsilon\sim\mathcal{N}(0,2\text{diag}(0.01,0.05,....,2.01))$ is the noise.

Figure 9: Rank vs learning rate in a vanilla matrix factorization problem for,
from upper left to lower right, $\mu=0.05,\ 0.15,\ 0.25,\ 0.35$. The
theoretical curve is from the commutation approximation where each subspace of
the model collapses at the critical learning rate
$\lambda=-2\frac{\mathbb{E}[h(x)]}{\mathbb{E}[h^{2}(x)]}$.

##  Appendix B Experiment with Adam

Figure 10: Rank of the converged solution for two-layer linear (upper left),
tanh (upper right), relu (lower left) and swish (lower right) models.

###  B.1 Experiment with Adam

We note that the phenomena we studied is not just a special feature of the
SGD, but, empirically, seems to be a universal feature of first-order
optimization methods that rely on minibatch sampling. Here, we repeat the
experiment in Figure 2. We train with the same data and training procedure,
except that we replace SGD with Adam (Kingma & Ba, 2014), the most popular
first-order optimization method in deep learning. Figure 10 shows that
similarly to SGD, Adam also converges to the low-rank saddles in similar
regions of learning rate and $\mu$.

##  Appendix C Additional Theoretical Concerns

###  C.1 A Simple example of Interpolation Minimum

To illustrate different concepts of stability of the SGD algorithm, we examine
a simple one-dimensional linear regression problem. The training loss function
for this problem is defined as
$L(w)=\frac{1}{N}\sum_{i}^{N}(wx_{i}-y_{i})^{2}$.

For GD, the dynamics diverges when the learning rate is larger than twice the
inverse of the largest Hessian eigenvalue. To see this, let
$H=\mathbb{E}_{x}[\hat{H}(w^{*},x)]$ denote the Hessian of $L$ and $h$ its
largest eigenvalue. Using GD leads to $||w_{t+1}||=||w_{0}(I-\lambda
H)^{t}||\propto|1-\lambda h|^{t}$. Divergence happens when $|1-\lambda h|>1$.
The range of viable learning rates is thus:

$\lambda_{\rm GD}\leq 2/h=2/\mathbb{E}_{x}[x^{2}].$ (14)

Naively, one would expect that a similar condition approximates the stability
condition for the case when mini-batch sampling is used to estimate the
gradient.

For SGD, the variance stability condition is the same as the condition that
the second moment of SGD decreases after every time step, starting from an
arbitrary initialization (see Section C.1.1):

$\lambda_{\rm
DS}\leq\frac{2S^{2}\mathbb{E}[x^{2}]}{\mathbb{E}[x^{4}]+(S-1)^{2}\mathbb{E}[x^{2}]^{2}}.$
(15)

Also related is the stability condition proposed by Ziyin et al. (2022b), who
showed that starting from a stationary distribution, $w$ stays stationary
under the condition $\lambda_{\rm SS}<\frac{2}{h}\frac{1}{1+1/S}$, which we
call the stationary linear stability condition (SS). When we have batch size
$1$, the stability condition halves: $\lambda<1/h$. For all stability
conditions, we denote the maximum stable learning rate with an asterisk as
$\lambda^{*}$.

For the probabilistic stability, let us consider the interpolation regime,
where all data points $(x,y)\in\mathbb{R}^{2}$ lie on a straight line. In this
situation, the loss function has a unique global minimum of
$w^{*}=y_{i}/x_{i}$ for any $i$. Applying Theorem 1, one can immediately prove
the following proposition.

######  Proposition 4.

Let $\lambda$ be such that $\mathbb{E}_{x}[\log|1-\lambda x^{2}|]\neq 0$.
Then, for any $w_{0}$, $w_{t}\to_{p}w^{*}$ if and only if
$\mathbb{E}_{x}[\log|1-\lambda x^{2}|]<0$.

It is worth remarking that this condition is distinctively different from the
case in which the gradient noise is a parameter-independent random vector. For
example, Liu et al. (2021) showed that if the gradient noise is a parameter-
independent Gaussian, SGD diverges in distribution if $\lambda>2/h$. This
suggests that the fact that the noise of SGD is $w$-dependent is crucial for
its probabilistic stability.

One of the implications of the probabilistic stability is that for
$\lambda=1/x_{i}^{2}$, the SGD dynamics is always stable. Therefore, the
largest stable learning rate is roughly given by:

$\lambda_{\max}=1/x^{2}_{\min}.$ (16)

However, for these special choices of learning rates, the moment stability is
not always guaranteed. As mentioned earlier, convergence in mean occurs when
$\lambda\leq\lambda_{DS}^{*}$. However, this condition does not hold when
$\lambda=1/x^{2}_{\min}$ and $x_{\min}<\mathbb{E}[x_{i}]/2$, which is often
the case for standard datasets. This result shows that the maximal learning
rate that ensures stable training can be much larger than the maximal learning
rate required for convergence in mean (cf. (15)). For a fixed value of
$\mathbb{E}[x^{2}]$, $x_{\min}^{2}$ can be arbitrarily small, which means that
the maximal stable learning rate can be arbitrarily large. Another consequence
of this result is that the stability of SGD depends strongly on individual
data points and not just on summary statistics of the whole dataset.

This result highlights the importance of the Lyapunov exponent
$\Lambda=\mathbb{E}_{x}[\log|1-\lambda x^{2}|]$ and its sign in understanding
the convergence of $w_{t}$ to the global minimum. When $\Lambda$ is negative,
the convergence to the global minimum occurs. If $\Lambda$ is positive, SGD
becomes unstable. We can determine when $\Lambda$ is negative for a training
set of finite size by examining the following equation:

$\Lambda=\frac{1}{N}\sum_{i}\log|1-\lambda x_{i}^{2}|,$ (17)

which is negative when $\lambda$ is close to $1/x_{i}^{2}$ for some
$i\in{1,\dots,N}$. What is the range of $\lambda$ values that satisfy this
condition? Suppose that $\lambda$ is in the vicinity of some $1/x_{i}^{2}$:
$\lambda=\delta\lambda+1/x_{i}^{2}$, and the instability is caused by a single
outlier data point $x_{\rm out}\gg 1$. Then, $\Lambda$ is determined by the
competing contributions from the outlier, which destabilizes training, and
$x_{i}^{2}$, which stabilizes training, and the resulting condition is
approximately $|1-\lambda x_{i}^{2}|<{1}/{|\lambda x_{\rm out}^{2}|}$. Because
$\lambda\approx 1/x_{i}^{2}$, this condition leads to:

$|\delta\lambda|<{x_{i}^{2}}/{x_{\rm out}^{2}}.$ (18)

This is a small quantity. However, if we change the learning rate to the
stability region associated with another data point $x_{j}$ as soon as we exit
the stability region of $x_{i}$, we still maintain stability. Therefore, the
global stability region depends on the density of data points near $x_{i}$.
Assuming that there are $N$ data points near $x_{i}$ with a variance of
$\sigma^{2}$, the average distance between $x_{i}$ and its neighbors is
approximately $\sigma^{2}/N$. As long as $\sigma^{2}/N<{x_{i}^{2}}/{x_{\rm
out}^{2}}$, SGD will remain stable in a large neighborhood. In practical
terms, this means that when the number of data points is large, SGD is highly
resilient to outliers in the data as shown in Figure 11. We see that the
region of convergence in probability is dramatic, featuring stripes of
convergent regions that correspond to $1/x_{i}^{2}$ for each data point and
divergent regions where $\Lambda>0$. While simple, this example has a
fundamental implication: there are problems that cannot be learned by SGD at a
small learning rate but can be learned by SGD at a finite learning rate. This
implies the insufficiency of the commonly used stochastic differential
equation theories of SGD (Li et al., 2021).

(a) $N=2$

(b) $N=6$

(c) $N=500$

(d) $N=500$

Figure 11: Stability of SGD against a single outlier data in a dataset of size
$N$. Yellow denotes where SGD converges in probability, and dark blue denotes
divergence. We control the norm of the first data point ($x_{1}^{2}$) while
sampling the rest data from a standard normal distribution. (a-c) Stability of
SGD for different sizes of the dataset; (d) zoom-in of (c) at a small learning
rate. The grey dashed curves show $\lambda_{GD}^{*}$, and the green dashed
curve shows $\lambda_{GD}^{*}/N$. The intermediate finite learning rates are
robust against outliers in the data, whereas the smallest learning rates are
strongly sensitive to outliers in the data.

An important implication is the robustness of SGD to outliers in comparison to
gradient descent. As Figure 11 shows, the bulk region of probabilistic
stability stays roughly unchanged as the outlier data point becomes larger and
larger; in contrast, both $\lambda_{GD}^{*}$ and $\lambda_{DS}^{*}$ decreases
quickly to zero. In the bulk region of the learning rates, SGD is thus
probabilistically stable but not stable in the moments. Meanwhile, in sharp
contrast to this bulk robustness is the sensitivity of the smallest branch of
learning rates of SGD to the outliers. Assuming that there is an outlier data
point with a very large norm $c\gg N$, the largest $\lambda_{\rm GD}$ scales
as $\lambda_{\max}\sim Nc^{-1}$. In contrast, for SGD, the smallest branch of
the probabilistically stable learning rate scales as $c^{-1}$, independent of
the dataset size. This means that if we only consider the smallest learning
rate, SGD is much less stable than GD, and one needs to use a much smaller
learning rate to ensure convergence. For $\lambda_{\rm DS}$ a detailed
analysis in Section C.1.1 shows that $\lambda_{\rm DS}^{*}=(Nc)^{-1}$. Thus,
the threshold of convergence in mean square is yet one order of magnitude
smaller than that of probabilistic convergence. In the limit $N\to\infty$, SGD
cannot converge in variance but can still converge in probability.

####  C.1.1 $L_{2}$ Stability Conditions and the Derivation of Eq. (15)

For a general batch size $S$, the dynamics of SGD reads

$\displaystyle w_{t+1}$ $\displaystyle=w_{t}-\lambda
w_{t}\frac{1}{S}\sum_{i=1}^{S}x_{i}^{2}$ (19)
$\displaystyle=w_{t}\left(1-\lambda\frac{1}{S}\sum_{i=1}^{S}x_{i}^{2}\right).$
(20)

The second moment of $w_{t+1}$ is

$\displaystyle\mathbb{E}_{x}[w_{t+1}^{2}|w_{t}]$
$\displaystyle=w_{t}^{2}\mathbb{E}_{x}\left(1-\lambda\frac{1}{S}\sum_{i}^{S}x_{i}^{2}\right)^{2}$
(21)
$\displaystyle=w_{t}^{2}\left(1-2\lambda\mathbb{E}[x^{2}]+\frac{\lambda^{2}}{S^{2}}\sum_{i,j}^{S}\mathbb{E}[x_{i}^{2}x_{j}^{2}]\right)$
(22)
$\displaystyle=w_{t}^{2}\left(1-2\lambda\mathbb{E}[x^{2}]+\frac{\lambda^{2}}{S^{2}}\mathbb{E}[x^{4}]+\frac{\lambda^{2}(S-1)^{2}}{S^{2}}\mathbb{E}[x^{2}]^{2}\right).$
(23)

Note that this equation applies to any $w_{t}\in\mathbb{R}$. Therefore, the
second moment of $w_{t}$ is convergent if

$1-2\lambda\mathbb{E}[x^{2}]+\frac{\lambda^{2}}{S^{2}}\mathbb{E}[x^{4}]+\frac{\lambda^{2}(S-1)^{2}}{S^{2}}\mathbb{E}[x^{2}]^{2}<1,$
(24)

which solves to

$\lambda<\frac{2S^{2}\mathbb{E}[x^{2}]}{\mathbb{E}[x^{4}]+(S-1)^{2}\mathbb{E}[x^{2}]^{2}}.$
(25)

This condition applies to any data distribution. One immediate observation is
that it only depends on the second and fourth moments of the data distribution
and that both moments need to be finite for convergence at a non-zero learning
rate. It is quite instructive to solve this condition under a few special
conditions.

##### Gaussian Data Distribution

The condition (25) takes a precise form when the data is Gaussian. Using the
fact that for a Gaussian variable $x$ with variance $\sigma^{2}$,
$\mathbb{E}[x^{4}]=3\sigma^{4}$, the condition simplifies to

$\lambda<\frac{2}{\mathbb{E}[x^{2}]}\frac{S^{2}}{3S+(S-1)^{2}}.$ (26)

This is the same as Eq. (15).

##### Bernoulli Dataset

Another instructive case to consider is the case when there are only two data
points in the data: $x_{1}$ and $x_{2}$. The moments are

$\begin{cases}\mathbb{E}[x^{2}]=\frac{1}{2}(x_{1}^{2}+x_{2}^{2}),\\\
\mathbb{E}[x^{4}]=\frac{1}{2}(x_{1}^{4}+x_{2}^{4}).\end{cases}$ (27)

When one of the data points, say $x_{1}$, is large, the condition becomes

$\lambda<\frac{2S^{2}x_{1}^{2}}{2x_{1}^{4}+(S-1)^{2}x_{1}^{4}}=\frac{2S^{2}}{x_{1}^{2}}\frac{1}{2+(S-1)^{2}}.$
(28)

##### Extreme Outlier

We can also consider the general case of a finite dataset with a large
outlier, $x_{\rm max}$. The condition is similar to the Bernoulli case. We
have

$\begin{cases}\mathbb{E}[x^{2}]\approx\frac{1}{N}x_{\rm max}^{2},\\\
\mathbb{E}[x^{4}]\approx\frac{1}{N}x_{\rm max}^{4}.\end{cases}$ (29)

The condition reduces to

$\lambda<\frac{2S^{2}}{Nx_{\rm max}^{2}}\frac{1}{1+\frac{(S-1)^{2}}{N}}.$ (30)

This can be seen as the generalization of the Bernoulli condition. When $S=1$,
this condition becomes

$\lambda<\frac{2}{Nx_{\rm max}^{2}}.$ (31)

There are many other interesting limits of this condition we can consider from
the perspective of extreme value theory. However, this is beyond the scope of
this work and we leave it as an interesting future work.

###  C.2 Proofs

####  C.2.1 Proof of Theorem 1

###### Proof.

Consider Eq. (1):

$\theta_{t+1}=\theta_{t}-\lambda\hat{H}(x)(\theta_{t}-\theta^{*}).$ (32)

Defining $w_{t}=\theta_{t}-\theta^{*}$, this equation can be written as

$w_{t+1}=w_{t}-\lambda\hat{H}(x)w_{t},$ (33)

which is mathematically equivalent to the case when $\theta^{*}=0$. Therefore,
without loss of generality, we write the dynamics in the form of Eq. (33) in
this proof and the rest of the proofs.

Now, when $\hat{H}\propto nn^{T}$ is rank-$1$, we can multiply $n^{T}$ from
the left on both sides of the dynamics to obtain

$n^{T}w_{t+1}=n^{T}w_{t}-\lambda h(x)n^{T}w_{t}.$ (34)

The dynamics thus becomes one-dimensional in the direction of $n^{T}$.

Let $h_{t}$ denote the eigenvalue of the Hessian of the randomly sampled batch
at time step $t$. The dynamics in Eq. (3) implies the following dynamics

$\|n^{T}w_{t+1}\|/\|n^{T}w_{t}\|=|1-\lambda h_{t}|,$ (35)

which implies

$\|n^{T}w_{t+1}\|/\|n^{T}w_{0}\|=\prod_{\tau=1}^{t}|1-\lambda h_{\tau}|.$ (36)

We can define auxiliary variables
$z_{t}:=\log(\|n^{T}w_{t+1}\|/\|n^{T}w_{0}\|)-m$ and
$m:=\mathbb{E}[\log(\|n^{T}w_{t+1}\|/\|n^{T}w_{0}\|)]=t\mathbb{E}_{x}[\log|1-\lambda
h_{t}|]$. Let $\epsilon>0$. We have that

$\displaystyle\mathbb{P}(\|n^{T}w_{t}\|<\epsilon)$
$\displaystyle=\mathbb{P}(\|n^{T}w_{0}\|e^{z_{t}+m}<\epsilon)$ (37)
$\displaystyle=\mathbb{P}\left(\frac{1}{{t}}z_{t}<\frac{1}{{t}}(\log\epsilon/\|n^{T}w_{0}\|-m)\right)$
(38)
$\displaystyle=\mathbb{P}\left(\frac{z_{t}}{t}<-\mathbb{E}_{x}\log|1-\lambda
h_{t}|+o(1)\right).$ (39)

By the law of large numbers, the left-hand side of the inequality converges to
$0$, whereas the right-hand side converges to a constant. Thus, we have, for
all $\epsilon>0$,

$\lim_{t\to\infty}\mathbb{P}(\|n^{T}w_{t}\|<\epsilon)=\begin{cases}1&\text{if
$m<0$};\\\ 0&\text{if $m>1$}.\end{cases}$ (40)

This completes the proof. ∎

####  C.2.2 Proof of Proposition 2

###### Proof.

Part 2 of the proposition follows immediately from Proposition 5, which we
prove below. Here, we prove part 1.

It suffices to consider a dataset with two data points for which
$h(x_{1})=1/\lambda$ and $h(x_{2})=c_{0}$, where each data point is sampled
with equal probability. Let $c_{0}$ be such that

$|1-\lambda c_{0}|^{p}>\frac{1}{2}.$ (41)

Now, we claim that this dynamics converges to zero in probability. To see
this, note that

$\|z_{t+1}\|=\begin{cases}\|z_{t}\||1-\lambda/\lambda|=0&\text{with
probability $0.5$};\\\ \|z_{t}\||1-\lambda c_{0}|&\text{with probability
$0.5$}.\end{cases}$ (42)

Therefore, at time step $t$, $\mathbb{P}(z_{t}=0)\geq 1-2^{-t}$, which
converges to $0$. This means that $z_{t}$ converges in probability to $0$.

Meanwhile, the $p$-norm is

$\displaystyle\mathbb{E}[\|z_{t+1}\|^{p}]$
$\displaystyle=\frac{1}{2}\mathbb{E}[\|z_{t}\|^{p}]|1-\lambda c_{0}|^{p}$ (43)
$\displaystyle\propto\frac{1}{2^{t}}|1-\lambda c_{0}|^{pt}\to 0.$ (44)

The convergence to zero follows from the construction that $|1-\lambda
c_{0}|^{p}>\frac{1}{2}$. This completes the proof. ∎

######  Proposition 5.

(No convergence in $L_{p}$.) For every strict saddle point $\theta^{*}$, there
exists an initialization $\theta_{0}$ such that for any
$\lambda\in\mathbb{R}_{+}$ and distance function $f(\cdot,\theta^{*})$,
$\theta^{*}$ is unstable in $f$.

###### Proof.

This problem is easy to prove when $\theta$ is one-dimensional. For a high-
dimensional $\theta$, the dynamics of SGD is

$\theta_{t+1}=(I-\lambda\hat{H}_{t})\theta_{t}.$ (45)

Note that the expected value of $\theta_{t}$ is the same as the gradient
descent iterations:

$\mathbb{E}[\theta_{t+1}]=(I-\lambda\mathbb{E}[\hat{H}])\mathbb{E}[\theta_{t}]=(I-\lambda\mathbb{E}[\hat{H}])^{t}\theta_{0},$
(46)

which diverges if $\theta_{0}$ is in one of the escape directions of
$\mathbb{E}[\hat{H}]$, which exist by the definition of strict saddle points.

Taking the $f-$distance of both sides and taking expectation, we obtain

$\displaystyle\mathbb{E}[f(\theta_{t},\theta^{*})]$ $\displaystyle\geq
f(\mathbb{E}[\theta_{t}],\theta^{*})$ (47)
$\displaystyle=f\left((I-\lambda\mathbb{E}[\hat{H}])^{t}\theta_{0},\theta^{*}\right)\nto
0.$ (48)

The first line follows from the fact that the distance function is convex by
definition, and so one can apply Jensen’s inequality.

Therefore, as long as $\theta_{0}$ overlaps with the concave directions of
$\mathbb{E}[\hat{H}]$, the argument of $f$ diverges, which implies that the
distance function converges to a nonzero value. The expected value of
$\theta_{t}$ is just the gradient descent trajectory, which diverges for any
strict saddle point.

By definition, $\mathbb{E}[\hat{H}]$ contains at least one negative
eigenvalue, and so the directions that do not overlap with this direction are
strict linear subspaces with dimension lower than the the total available
dimensions. This is a space with Lesbegue measure zero. The proof is complete.
∎

####  C.2.3 Proof of Theorem 2

Let us first state the Furstenberg-Kesten theorem.

######  Theorem 3.

(Furstenberg-Kesten theorem) Let $X_{1},\ X_{2},\ X_{3},...$ be independent
random square matrices drawn from a metrically transitive time-independent
stochastic process and $\mathbb{E}[\log_{+}\|X^{1}\|<\infty]$,
then555$\log_{+}(x)=\max(\log x,0)$.

$\lim_{n\to\infty}\frac{1}{n}\log\|X_{1}X_{2}...X_{n}\|=\lim_{n\to\infty}\mathbb{E}\left[\frac{1}{n}\log\|X_{1}X_{2}...X_{n}\|\right]$
(49)

with probability $1$, where $\|\cdot\|$ denotes any matrix norm.

Namely, the Lyapunov exponent of every trajectory converges to the expected
value almost surely. Essentially, this is a law of large numbers for the
Lyapunov exponent.

Now, we present the proof of Theorem 2.

###### Proof.

First of all, we define $m_{t}=\log\|\theta_{t}-\theta^{*}\|$ and
$z_{t}=m_{t}-\mathbb{E}[m_{t}]$. By definition, we have

$\displaystyle\mathbb{P}(g_{t}<\epsilon)$
$\displaystyle=\mathbb{P}(e^{z_{t}+m_{t}}<\epsilon)$ (50)
$\displaystyle=\mathbb{P}\left(\frac{1}{{t}}(z_{t}+\mathbb{E}[m_{t}])<\frac{1}{{t}}\log\epsilon\right)$
(51)
$\displaystyle=\mathbb{P}\left(\frac{1}{{t}}(z_{t}+\mathbb{E}[m_{t}])<o(1)\right).$
(52)

We can lower bound this probability by

$\mathbb{P}\left(\frac{1}{{t}}(z_{t}+\mathbb{E}[m_{t}])<o(1)\right)\geq\mathbb{P}\left(\frac{1}{{t}}\max_{\theta_{0}}(z_{t}+\mathbb{E}[m_{t}])<o(1)\right).$
(53)

By the definition of SGD, we have

$\displaystyle\frac{1}{t}\max_{\theta_{0}}\left(z_{t}(\theta_{0})+\mathbb{E}[m_{t}]\right)$
$\displaystyle=\frac{1}{t}\max_{\theta_{0}}\log\left\|\prod_{i}^{t}(I-\lambda\hat{H}_{i})(\theta_{t}-\theta_{0})\right\|.$
(54)

By the Furstenberg-Kesten theorem (Furstenberg & Kesten, 1960), this quantity
converges to the constant
$\Lambda=\lim_{t\to\infty}\mathbb{E}[m_{t}]/t\in\mathbb{R}$ almost surely.
Namely, $z_{t}/t$ converges to $0$ for almost every SGD trajectory.

Thus, for every $\epsilon$, if $\Lambda<0$, Eq. (52) can be bounded as

$\lim_{t\to\infty}\mathbb{P}(g_{t}<\epsilon)=\mathbb{P}(\Lambda<0)=1.$ (55)

Because $\Lambda$ is a constant, we have that if $\Lambda<0$, all trajectories
from all initialization converge to $0$. This finishes the first part of the
proof. For the second part, simply let $z_{t}$ be the trajectory starting from
the trajectory that achieves the maximum Lyapunov exponent. Again, this
dynamics escapes with probability $1$ by the Furstenberg-Kesten theorem. The
proof is complete. ∎

####  C.2.4 Proof of Proposition 3

###### Proof.

We consider the dynamics of SGD around a saddle:

$\ell=-\chi\sum_{i}u_{i}w_{i},$ (56)

where we have combined $xy$ into a single variable $\chi$. The dynamics of SGD
is

$\begin{cases}w_{i,t+1}=w_{i,t}+\lambda\chi u_{i,t};\\\
u_{i,t+1}=u_{i,t}+\lambda\chi w_{i,t}.\end{cases}$ (57)

Namely, we obtain a set of coupled stochastic difference equations. Since the
dynamics is the same for all values of the index $i$, we omit $i$ from now on.
This dynamics can be decoupled if we consider two transformed parameters:
$h_{t}=w_{t}+u_{t}$ and $m_{t}=w_{t}-u_{t}$. The dynamics for these two
variables is given by

$\begin{cases}h_{t+1}=h_{t}+\lambda\chi h_{t};\\\ m_{t+1}=m_{t}-\lambda\chi
m_{t}.\end{cases}$ (58)

We have thus obtained two decoupled linear dynamics that take the same form as
that in Theorem 1. Therefore, as immediate corollaries, we know that $h$
converges to $0$ if and only if $\mathbb{E}[\log|1+\lambda\chi|]<0$, and $m$
converges to $0$ if and only if $\mathbb{E}[\log|1-\lambda\chi|]<0$.

When both $h$ and $m$ converge to zero in probability, we have that both $w$
and $u$ converge to zero in probability. For the data distribution under
consideration, we have

$\mathbb{E}[\log|1+\lambda\chi|]=\frac{1}{2}\log\left|(1+\lambda)(1+\lambda
a)\right|$ (59)

and

$\mathbb{E}[\log|1-\lambda\chi|]=\frac{1}{2}\log\left|({1-\lambda})({1-\lambda
a})\right|.$ (60)

There are four cases: (1) both conditions are satisfied; (2) one of the two is
satisfied; (3) neither is satisfied. These correspond to four different phases
of SGD around this saddle. ∎


# On the Antiferromagnetic-Ferromagnetic Phase Transition
in Pinwheel Artificial Spin Ice

Anders Strømberg<EMAIL_ADDRESS>Department of Electronic Systems,
Norwegian University of Science and Technology, Trondheim, Norway.  Einar
Digernes  Department of Electronic Systems, Norwegian University of Science
and Technology, Trondheim, Norway.  Rajesh Vilas Chopdekar  Advanced Light
Source, Lawrence Berkeley National Laboratory, Berkeley, CA, USA  Jostein
Grepstad  Department of Electronic Systems, Norwegian University of Science
and Technology, Trondheim, Norway.  Erik Folven  Department of Electronic
Systems, Norwegian University of Science and Technology, Trondheim, Norway.

(September 3, 2024)

###### Abstract

Nanopatterned magnetic thin films offer a platform for exploration of tailored
magnetic properties such as emergent long-range order. A prominent example is
artificial spin ice (ASI), where an arrangement of nanoscale magnetic
elements, acting as macrospins, interact via their dipolar fields. In this
study, we discuss the transition from antiferromagnetic (AF) to ferromagnetic
(FM) long-range order in a square lattice ASI, as the magnetic elements are
gradually rotated through $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ to a
“pinwheel” configuration. The AF–FM transition is observed experimentally
using synchrotron radiation x-ray spectromicroscopy and occurs for a critical
rotation angle of the nanomagnets, contingent upon the dipolar coupling
determined by their separation in the lattice. Large-scale magnetic dipole
simulations show that the point-dipole approximation fails to capture the
correct AF–FM transition angle. However, excellent agreement with experimental
data is obtained using a dumbbell-dipole model which better reflects the
actual dipolar fields of the magnets. This model explains the coupling-
dependence of the transition angle, another feature not captured by the point-
dipole model. Our findings resolve a discrepancy between measurement and
theory in previous work on “pinwheel” ASIs. Control of the AF–FM transition
and this revised model open for improved design of magnetic order in
nanostructured systems.

††preprint: APS/123-QED

Artificial spin ice (ASI) systems[1] composed of single-domain nanomagnets,
acting as macrospins, exhibit a wide range of exotic behaviors[2, 3, 4, 5].
Recent work suggests that these systems are key to emergent device
technologies such as physical reservoir computing[6, 7]. The interaction
between the macrospins in an ASI is mediated by the stray fields of the
individual magnets, creating a complex energy landscape sensitive to the exact
shape of the elements and the geometric arrangement of the ensemble[3, 8].
With precise nanofabrication, this sensitivity enables effective tuning of
magnetic order. Accurate modeling of the magnetic interactions in these
dipolar-coupled metamaterials is essential to the analysis of emergent
behavior, such as long-range magnetic order.

Previous research has established the predominance of long-range ferromagnetic
(FM) order in a $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ “pinwheel” ASI[9],
in contrast to the antiferromagnetically (AF) ordered ground state of the
square ASI [1]. The FM-ordered pinwheel ASI is a metamaterial with several
exotic emergent properties such as dynamic chirality[10], unique domain-wall
topologies[11], and Heisenberg pseudo-exchange[12].

In this Letter, we investigate the AF to FM ordering transition we observe
when going from a square to a pinwheel ASI. Specifically, we examine
modifications of the square ASI by gradually rotating each nanomagnet on the
square lattice from $\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$ to
$\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$, see fig. 1. We observe the
long-range magnetic order of these metamaterials using x-ray magnetic circular
dichroism (XMCD) spectromicroscopy at the Advanced Light Source (ALS). A
predominance of FM order is seen for the ASIs with a $\theta$ close to the
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ pinwheel configuration. We find a
critical angle for the AF–FM transition dependent on the separation of the
nanomagnetic elements, at odds with previous reports relying on a simple
point-dipole model[9, 13]. For a more accurate representation of the stray
fields, we invoke a dumbbell-dipole model, which is found to accurately
capture the dependence of the AF–FM transition angle on the nanomagnet
separation.

Figure 1:  Geometric configuration of square-to-pinwheel structures. (a)
Rotation of elements from the square
$\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$ ASI (top left) to the
typical $\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$ pinwheel ASI
(bottom right). (b) The same series as in (a), but with the lattice compressed
to compensate for decoupling. (c) Schematic showing the geometric parameters
of the systems.  Figure 2: Relaxed states of a square
($\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$) and pinwheel
($\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$) ASI. (a) and (b) XMCD-
PEEM images, with top left corner insets of a SEM image showing the rotation
and shape of the nanomagnets. (c) and (d) Analyzed magnetization maps. Gray
areas indicate the two possible phases of antiferromagnetically ordered
regions. Colored regions are ferromagnetically ordered with a net
magnetization as depicted in the legend. On the right, the relation between
nanomagnet states and cell magnetization is illustrated. The length of the
cell arrows correspond to the magnitude of the cell magnetization.

For this study, ensembles of $25\times 25$ nanomagnets with lateral dimensions
$l\times w=\qtyproduct{220x80}{\nano}$ and thickness
$t=$3\text{\,}\mathrm{nm}$$ were patterned using electron beam lithography
(EBL) and lift-off. Each sample featured arrays with a range of different
element rotations, $\theta$, corresponding to the schematics in fig. 1. The
magnetic material was deposited by e-beam evaporation of permalloy (Ni80Fe20),
followed by a $2\text{\,}\mathrm{nm}$ Al layer serving as an oxidation
barrier. The thickness of the magnetic film was chosen so that the ensembles
could be easily thermalized in situ.

Figure 3: Measured fraction of antiferromagnetic Voronoi cells for the
experimental ensembles. Five different separations $d$, kept constant across
the rotation series from
$0\text{\,}\mathrm{\SIUnitSymbolDegree}45\text{\,}\mathrm{\SIUnitSymbolDegree}$,
were investiaged. XMCD-PEEM images (top panel) show the PEEM images and
graphical representation of the relaxed states for the
$d=$185\text{\,}\mathrm{nm}$$ series.

Rotation of the elements in this ASI system (see fig. 1(a)) alters the
effective coupling between the nanomagnets due to the concentrated stray
fields at their ends. Increasing the rotation angle $\theta$ on a fixed
lattice will effectively decouple the elements. In this study, we reduce the
lattice constant with increasing rotation angle $\theta$, in order to retain a
strong interisland coupling, see fig. 1(b). We keep the parameter $d$, shown
in fig. 1(c), fixed for all element rotations.

Modeling the AF–FM transition of a full array presents a computational
challenge because of the large differences in scale between the internal
magnetization of a single nanomagnet to the macroscopic scale of the
nanomagnet array. In order to establish an accurate analytical stray-field
model, the stray field of an isolated single nanomagnet is calculated using a
micromagnetic simulator[14, 15, 16]. However, this micromagnetic approach is
not feasible for the interactions of the full ensemble, due to the high
computational cost, and an analytical representation of the stray field for a
single nanomagnet is required. We use the micromagnetic framework
$\text{mumax}^{3}$[14, 15] for the micromagnetic calculations and the Ubermag
package[16] for comparison with analytically expressed fields. The relaxed
magnetization in the absence of an external field was computed using
micromagnetic simulation, and the stray field was calculated from the
resulting magnetization texture. Each nanomagnet was modeled as a 220x80
rectangle with rounded corners and a thickness of $3\text{\,}\mathrm{nm}$. The
saturation magnetization was set at
$M_{\text{S}}=$860\text{\,}\mathrm{kA}\text{\,}{\mathrm{m}}^{-1}$$ and the
exchange stiffness at
$A_{\text{ex}}=$1.3\text{\times}{10}^{-11}\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-1}$$[17].
A simulation mesh of $512\times 512$ cells was adopted, where each cell had a
side length of $2\text{\,}\mathrm{nm}$, well below the exchange length
$l_{\text{ex}}=\sqrt{2A_{\text{ex}}/\mu_{0}M_{\text{S}}^{2}}\approx$5\text{\,}\mathrm{nm}$$.

X-ray magnetic circular dichroism photoemission electron microscopy (XMCD-
PEEM) was carried out on the PEEM3 beamline at ALS and was used to image the
magnetic state of the ensembles. Prior to imaging at room temperature, the
samples were brought to $T=$385\text{\,}\mathrm{K}$$ to thermalize the system,
so as to avoid quenched, metastable, high-energy states. Complete
thermalization was confirmed from the XMCD-PEEM images, showing no magnetic
contrast, thus indicating full superparamagnetic behavior on the time scale of
imaging. After thermalization, the sample was gradually cooled to room
temperature, resulting in a relaxed state governed by the intrinsic dipolar
interactions

Magnetic contrast images of the
$\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and
$\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$ ASIs are shown in fig. 2(a)
and (b) with the typical AF and FM order, respectively. The magnetization of
each nanomagnet is identified from the XMCD-PEEM images, with the help of
machine learning, to create a graphical representation of the ensemble state.
In fig. 2(c) and (d), we divide the ASI into so-called Voronoi cells (VCs),
where each cell is assigned the net magnetization of four nearest neigbor
nanomagnets. The $\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$ ASI in fig.
2(c) shows extended regions with zero net magnetic moment, AF-ordered domains
divided by domain walls of finite magnetic moment. In contrast, the
$\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$ ASI in fig. 2(d) displays
numerous smaller FM-ordered regions with a net magnetic moment, mimicking the
domains of a conventional ferromagnetic material. As these square and pinwheel
ASI ensembles show distinct AF and FM order, respectively, this finding
suggests successful relaxation of the arrays and sufficient dipolar magnetic
coupling, to study their long-range order.

The main results of the measurement series, where the element rotation
$\theta$ is gradually increased while $d$ is kept constant, are shown in fig.
3. (A complete overview of the XMCD-PEEM images discussed in this study is
shown in the Supplemental Material[18].) The AF-ordered VCs dominate for all
rotation angles below $30\text{\,}\mathrm{\SIUnitSymbolDegree}$ irrespective
of separation $d$, except for the decoupled $d>$200\text{\,}\mathrm{nm}$$
system. The transition to FM ordering starts at around
$\theta=$39\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for the ASI with
$d=$195\text{\,}\mathrm{nm}$$. The strongest coupled system, with
$d=$170\text{\,}\mathrm{nm}$$, exhibits the transition at
$\theta=$43\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The transition angle,
$\theta_{\text{t}}$, is seen to decrease monotonically with reduced magnetic
coupling. We note that the observed transition angles differ significantly
from the $\theta_{\text{t}}\approx$35\text{\,}\mathrm{\SIUnitSymbolDegree}$$
reported in previous work[9, 13] where a simple point-dipole model was used.

Figure 4: Ensemble state energies and nanomagnet stray-fields. Two example
states are chosen when calculating the energies, a $100\%$ AF-ordered state
and a polarized state, i.e., a $100\%$ FM-ordered state. (a) energies for the
polarized and AF states, as calculated by the point-dipole model. (b) energies
for the polarized, AF, and experimental states, as calculated by the dumbbell-
dipole model. (c), (d), and (e) show the calculated stray fields for a single
220x80x3 nanomagnet (white rectangle), for the point-dipole model, the
micromagnetic simulation, and the dumbbell-dipole model, respectively. The
field strength color bar is truncated at $5\text{\,}\mathrm{mT}$.

An aberration from the expected complete AF ordering is found for
$\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$$. For this element rotation,
we find a multidomain AF ordering with FM ordered domain walls, which may be
explained by reduced coupling across the nanomagnet corners.

In the AF-ordered square ASI, the next nearest parallel oriented nanomagnets
feature antiparallel spin alignment. The change in long-range order from AF to
FM with $\theta$ increasing from
$$0\text{\,}\mathrm{\SIUnitSymbolDegree}$\rightarrow$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$,
implies that for $\theta>\theta_{\text{t}}$ no neighbors align antiparallel.
Specifically, the antiparallel magnetization alignment of the next nearest
neighbors (a defining feature of AF order in the square ASI) will flip when
their stray-field interaction energy favors parallel spin alignment. With this
assumption, we can estimate the transition angle of an ensemble
$\theta_{\text{t}}$ based on an analytical expression for the stray field.

With the ASI nanomagnets represented by point dipoles, the stray field is
described by,

$\vec{H}_{\text{point
dipole}}(\vec{r})=\frac{1}{4\pi}\left[\frac{3(\vec{m}_{\text{dip}}\cdot\vec{r})\vec{r}}{r^{5}}-\frac{\vec{m}_{\text{dip}}}{r^{3}}\right].$
(1)

where $\vec{m}_{\text{dip}}$ is the magnetic moment of a single nanomagnet and
$\vec{r}$ is the distance vector from the point dipole. By calculating the
angle at which the point-dipole interaction energy for the parallel neighbors
changes sign, favoring parallel alignment, we determine a transition at
$\theta_{\text{t}}=\arctan(1/\sqrt{2})\approx$35.3\text{\,}\mathrm{\SIUnitSymbolDegree}$$,
commonly referred to as the ‘magic angle’, $\theta_{\text{magic}}$[19, 20].

When modeling the full ensembles in the point-dipole approximation, the AF–FM
transition takes place at
$\theta\approx$35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$$, coinciding with
$\theta_{\text{magic}}$, independent of the interisland coupling distance $d$.
At this angle, the nature of the interaction energy for parallel neighbors
changes. The minimum energy state shifts from “two aligned–two anti-aligned”
to “all aligned” macrospins, as a result of the asymmetry between the $x$ and
$y$ components of the point-dipole field (eq. 1). This angle was cast as the
critical angle for the transition between AF and FM order in previous work
relying on Monte-Carlo simulations within the point-dipole model[9, 13].
However, the point-dipole approximation misses key features of the actual
stray-field coupling as it fails to capture the observed variation in
transition angle $\theta_{\text{t}}$ with nanomagnet separation.

The stray-fields of the point-dipole model and the micromagnetic simulation
are shown in fig. 4(c) and (d), respectively. By comparison with the
micromagnetic stray field, we find that a dumbbell-dipole model[21] captures
the characteristics of the nanomagnet stray field better than the point-dipole
model (see fig. 4(e)).

The dumbbell model emulates two magnetic charges, $\pm q_{\text{m}}$,
separated by a distance $d_{\text{dip}}$ on the long-axis of the magnet. The
magnitude of the total magnetic moment of a nanomagnet,
$\lvert{\vec{m}_{\text{dip}}}\rvert$, is given by $M_{\text{S}}\cdot V$, where
$M_{\text{S}}$ represents the saturation magnetization of the magnetic
material comprising the nanomagnet and $V$ is the nanomagnet volume. The
effective stray field of the single magnet at a position $\vec{r}$ relative to
its center in this model is given by

$\displaystyle\vec{H}_{\text{dumbbell}}(\vec{r})$
$\displaystyle=\vec{H}_{\text{monopole}}(\vec{r}-\vec{a})-\vec{H}_{\text{monopole}}(\vec{r}+\vec{a}),$
$\displaystyle=\frac{q_{m}}{4\pi}\left[\frac{\vec{r}-\vec{a}}{\left|\vec{r}-\vec{a}\right|^{3}}-\frac{\vec{r}+\vec{a}}{\left|\vec{r}+\vec{a}\right|^{3}}\right],$
(2)

where $\vec{a}=d_{\text{dip}}\cdot\hat{d}_{\text{dip}}/2$ is the distance
vector separating the center of the magnet from its positive end. For the
systems considered in this letter, we find that a naive, unoptimized dumbbell
model of monopoles with $d_{\text{dip}}=l$ and $q_{\text{m}}\cdot
d_{\text{dip}}=\lvert{\vec{m}_{\text{dip}}}\rvert$ is sufficient to match the
experimental data.

The magnetic long-range order is evaluated by considering the total energies
for the investigated ASI rotation angles
$\theta=$0\text{\,}\mathrm{\SIUnitSymbolDegree}$\rightarrow$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$
within the scope of the dumbbell-dipole model. We find AF and FM long-range
order for angles close to $0\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$, respectively. Taking the transition
angle $\theta_{\text{t}}$ as the pinwheel rotation angle where the total
energy of a FM state drops below the energy of the AF state, the energies of
the ensemble for the point-dipole and dumbbell-dipole stray-field models,
respectively, are shown in fig. 4a-b. The state energy is calculated by
summation over all interactions for each of the experimentally observed states
(circular data points), as well as for the fully polarized and AF-ordered
states. The estimated transition angles from simulations and experimental
results are compared in table 1, with the simulated transition angles based on
the FM and AF energy calculations. The measured angles are those where the AF
fraction is interpolated at $50\%$.

$d$ (nm) | $E_{\text{AF}}=E_{\text{pol}}$ | 50% AF
---|---|---
| point dipole | dumbbell | measured
170 | $35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $43.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $43.0\text{\,}\mathrm{\SIUnitSymbolDegree}$
175 | $35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $43.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $42.1\text{\,}\mathrm{\SIUnitSymbolDegree}$
185 | $35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $42.7\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $40.7\text{\,}\mathrm{\SIUnitSymbolDegree}$
195 | $35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $42.0\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $39.1\text{\,}\mathrm{\SIUnitSymbolDegree}$
205 | $35.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $41.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | N/A
Table 1: Transition angles $\theta_{\text{t}}$, where the calculated energies
of the antiferromagnetic state ($E_{\text{AF}}$) is equal to the energy of the
polarized state ($E_{\text{pol}}$), for point-dipole and dumbbell-dipole
models at different $d$ values. Also included: experimental transition angle
taken to be the angle where the measured fraction of AF Voronoi cells is
$50\%$ (by interpolation).

From the stray-field maps in fig. 4(c-e), it is clear that the dumbbell model
better captures the characteristic anisotropy and magnitude of the calculated
stray fields of the micromagnetically modeled nanomagnet, in particular in the
near-field region depicted here. Moreover, dumbbell-dipole model accurately
captures the experimental observations with a transition angle
$\theta_{\text{t}}$ dependent on the intermagnet coupling. The AF–FM
transition is observed for angles ranging from
$39\text{\,}\mathrm{\SIUnitSymbolDegree}43\text{\,}\mathrm{\SIUnitSymbolDegree}$,
while simulations using this model indicate an AF–FM energy crossover for
angles from $41.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ to
$43.3\text{\,}\mathrm{\SIUnitSymbolDegree}$. This close correspondence
underscores the advantage of using the dumbbell-dipole model for an accurate
representation of the physical behavior of the system.

A notable exception to the AF order, observed at small angles $\theta$, is
found for the $d=$205\text{\,}\mathrm{nm}$$-series in fig. 3, where no
experimental transition is found. The ensembles with this nanomagnet
separation are all nearly polarized in the same direction. This finding
suggests that dipolar coupling in this particular ensemble is too small to
drive long-range order. The observed polarization may be due to a small stray
field present during annealing inside the PEEM-3 microscope.

In conclusion, this study shows a transition from distinct AF order to
pronounced FM order in a square-to-pinwheel artificial spin ice system. We
find that the critical transition angle is dependent on the interisland
coupling, a behavior not captured by the conventional point-dipole
approximation. By introducing a dumbbell-dipole model, we achieve excellent
agreement with experimental observations and micromagnetic simulations of the
near-field interactions. We maintain that this model offers a suitable
framework for simulation of a variety of nanomagnet shapes and textures. Our
approach improves modeling of magnetic ordering in ASIs, which is key to
fundamental research as well as technological applications.

###### Acknowledgements.

This research used resources of the Advanced Light Source, which is a DOE
Office of Science User Facilities under contract no. DE-AC02-05CH11231. The
Research Council of Norway is acknowledged for the support to the Norwegian
Micro- and Nanofabrication Facility, NorFab, project no. 295864. The samples
were fabricated at NTNU NanoLab. This work was funded by the Norwegian
Research Council through the IKTPLUSS project SOCRATES (Grant no. 270961) and
by the EU FET-Open RIA project SpinENGINE (Grant no. 861618).

## References

  * Wang  _et al._ [2006] R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, Nature 439, 303 (2006).
  * Mengotti  _et al._ [2011] E. Mengotti, L. J. Heyderman, A. F. Rodríguez, F. Nolting, R. V. Hügli, and H.-B. Braun, Nature Physics 7, 68 (2011).
  * Skjærvø  _et al._ [2020] S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Nature Reviews Physics 2, 13 (2020).
  * May  _et al._ [2021] A. May, M. Saccone, A. van den Berg, J. Askey, M. Hunt, and S. Ladak, Nature Communications 12, 3217 (2021).
  * Digernes  _et al._ [2021] E. Digernes, A. Strømberg, C. A. F. Vaz, A. Kleibert, J. K. Grepstad, and E. Folven, Applied Physics Letters 118, 202404 (2021).
  * Jensen  _et al._ [2018] J. H. Jensen, E. Folven, and G. Tufte, in  _ALIFE 2018: The 2018 Conference on Artificial Life_ (MIT Press, Tokyo, Japan, 2018) pp. 15–22.
  * Gartside  _et al._ [2022] J. C. Gartside, K. D. Stenning, A. Vanstone, H. H. Holder, D. M. Arroo, T. Dion, F. Caravelli, H. Kurebayashi, and W. R. Branford, Nature Nanotechnology 17, 460 (2022).
  * Digernes  _et al._ [2020] E. Digernes, S. D. Slöetjes, A. Strømberg, A. D. Bang, F. K. Olsen, E. Arenholz, R. V. Chopdekar, J. K. Grepstad, and E. Folven, Physical Review Research 2, 013222 (2020).
  * Macêdo  _et al._ [2018] R. Macêdo, G. M. Macauley, F. S. Nascimento, and R. L. Stamps, Physical Review B 98, 014437 (2018).
  * Gliga  _et al._ [2017] S. Gliga, G. Hrkac, C. Donnelly, J. Büchi, A. Kleibert, J. Cui, A. Farhan, E. Kirk, R. V. Chopdekar, Y. Masaki, N. S. Bingham, A. Scholl, R. L. Stamps, and L. J. Heyderman, Nature Materials 16, 1106 (2017).
  * Li  _et al._ [2019] Y. Li, G. W. Paterson, G. M. Macauley, F. S. Nascimento, C. Ferguson, S. A. Morley, M. C. Rosamond, E. H. Linfield, D. A. MacLaren, R. Macêdo, C. H. Marrows, S. McVitie, and R. L. Stamps, ACS Nano 13, 2213 (2019).
  * Paterson  _et al._ [2019] G. W. Paterson, G. M. Macauley, Y. Li, R. Macêdo, C. Ferguson, S. A. Morley, M. C. Rosamond, E. H. Linfield, C. H. Marrows, R. L. Stamps, and S. McVitie, Physical Review B 100, 174410 (2019).
  * Macauley  _et al._ [2020] G. M. Macauley, G. W. Paterson, Y. Li, R. Macêdo, S. McVitie, and R. L. Stamps, Physical Review B 101, 144403 (2020).
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, AIP Advances 4, 107133 (2014).
  * Leliaert  _et al._ [2018] J. Leliaert, M. Dvornik, J. Mulkers, J. D. Clercq, M. V. Milošević, and B. V. Waeyenberge, Journal of Physics D: Applied Physics 51, 123002 (2018).
  * Beg  _et al._ [2022] M. Beg, M. Lang, and H. Fangohr, IEEE Transactions on Magnetics 58, 1 (2022).
  * Maicas  _et al._ [2008] M. Maicas, R. Ranchal, C. Aroca, P. Sánchez, and E. López, The European Physical Journal B 62, 267 (2008).
  * [18] See Supplemental Material below for all XMCD-PEEM images with their graphical representation.
  * Erickson  _et al._ [1993] S. J. Erickson, R. W. Prost, and M. E. Timins, Radiology 188, 23 (1993).
  * De Sousa  _et al._ [2004] R. De Sousa, J. D. Delgado, and S. Das Sarma, Physical Review A 70, 052304 (2004).
  * Castelnovo  _et al._ [2008] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).

Supplemental Material

Figure S1: XMCD-PEEM images of the ASI systems in their relaxed state. All
separations $d$ and rotation angles $\theta$ included in this paper are
displayed. Figure S2: Analyzed magnetization maps of the ASI systems in their
relaxed state. All separations $d$ and rotation angles $\theta$ included in
this paper are displayed. Gray domains indicate antiferromagnetically ordered
regions, where each shade of gray indicates one of two possible
antiferromagnetic phases. Colored regions are ferromagnetically ordered with a
net magnetization as in the main paper.